%-----------------------------------------------------------------
% Author:       James Fishbaugh
% Date:         05/21/2008
% Description:  Reads experiment data and generates a plot
%               of a region of one standard deviation from the
%               mean for the different types of models used.
%----------------------------------------------------------------

%--------------------------------------------------------
% THIS IS THE FILE IO SECTION WHERE THE DATA FILES ARE
% READ AND VALUES STORED 
%--------------------------------------------------------

% The names of the folders (also prefix for filenames)
dirs = {'aas2' 'ahu1' 'ahu2' 'ahu2' 'bds1' 'cma1' 'cma2' 'dpo1' ...
        'fpa1' 'fpa2' 'fsa1' 'jvj1' 'lpe1' 'mbr1' 'mmr1' 'mtk1' ...
        'oga1' 'oso1' 'rhi1' 'rhi2' 'rhi3' 'rhi4' 'rhi5' 'spr1' ...
        'ssu1' 'tpa1' 'yuv1'};

% The entries in this vector correspond to the vector of directory names
% For example:  aas2 is full stokes, ahu1 is not, ahu2 is not, etc
% NOTE: THIS NEEDS TO BE CHECKED FOR ACCURACY
fullStokesMask = [1 0 0 0 0 1 0 0 0 1 0 1 0 0 0 0 1 0 1 0 1 0 0 1 1 0 1];

% The location of the data (this is the path to the directory where all the
% sub dirs lie)
% NOTE:  CHANGE THIS TO THE LOCATION OF YOUR DATA
location = '/Applications/MATLAB_R2007b/work/ismip_all/';

% The experiment we're interested in
exp = 'c';

% The file we're intersted in
fileNames = {[exp '005.txt'], [exp '010.txt'], [exp '020.txt'], ...
             [exp '040.txt'], [exp '080.txt'], [exp '160.txt']};

% The desired vector for normalized x.  Used for interpolation so every
% data set has the same number of points.
newX = linspace(0.0,1.0,250);
newY = 1/4*ones(250,1)';

% Will hold all the full stokes surface velocities
fsVels = cell(length(fileNames), length(fullStokesMask(fullStokesMask == 1)));
% Will hold all the non-full stokes surface velocities
nfsVels = cell(length(fileNames), length(fullStokesMask(fullStokesMask == 0)));

% Keeps track of the index we are currently on (this is necessary because
% some submittors didn't do all experiments so we can't rely on the loop
% counter to keep track of where we are at)
fsIndex = 1;
nfsIndex = 1;

% Keeps track of the number of submittors for each domain length (this will
% help us later on in allocating the correct size of data structures)
fsSizes = zeros(1, length(fileNames));
nfsSizes = zeros(1, length(fileNames));

% Loop over all folders
for i=1:length(dirs)
    
    % Loop over all files
    for j=1:length(fileNames)
    
        % Create the full file path from the base path and current directory
        fullFilePath = [location dirs(i) '/' dirs(i) fileNames(j)];

        % Open the current file for reading
        file = fopen(cell2mat(fullFilePath),'r');

        % If that file exists
        if (file ~= -1)

            % Creates a cell array from the file
            data = textscan(file, '%f %f %f %f %f %f %f %f %f %f');

            % Store the values 
            x = data{1,1};
            y = data{1,2};
            xvel = data{1,3};
            yvel = data{1,4};
            zvel = data{1,5};
            
            % Make sure the data is normalized
            x = x./max(x);
            
            % Make sure the data is normalized
            y = y./max(y);
            
            % Make sure there are no nans in the data
            x = inpaint_nans(x, 4);
            y = inpaint_nans(y, 4);
            xvel = inpaint_nans(xvel, 4);
            yvel = inpaint_nans(yvel, 4);
            zvel = inpaint_nans(zvel, 4);
            
            % Calculates the magnitude of the surface velocity
            surfVel = sqrt(xvel.^2 + yvel.^2 + zvel.^2);
            surfVel = sqrt(xvel.^2 + yvel.^2);

            % If the current data is from a full stokes
            if (fullStokesMask(i) == 1)
                
                % Store the data in the full stokes structure
                fsVels{j, fsIndex} = griddata(x, y, surfVel, newX, newY);
                fsVels{j, fsIndex} = inpaint_nans(fsVels{j,fsIndex}, 4);
                
                % Keep track of the number of entries in the array
                if (length(fsVels{j, fsIndex}) > 0)
                   
                    % Increment the number who performed this experiment
                    fsSizes(j) = fsSizes(j) + 1;
                
                end
                
                % If we are on the last file we can increment the index
                if (j == (length(fileNames)-1))
                    
                    % Increment the full stokes index
                    fsIndex = fsIndex + 1;
                
                end
                
            % Else its non-full stokes
            else
                
                % Store the data in the non-full stokes structure
                nfsVels{j, nfsIndex} = griddata(x, y, surfVel, newX, newY);
                nfsVels{j, nfsIndex} = inpaint_nans(nfsVels{j, nfsIndex}, 4);
                
                % Keep track of the number of entries in the array
                if (length(nfsVels{j, fsIndex}) > 0)
                   
                    % Increment the number who performed this experiment
                    nfsSizes(j) = nfsSizes(j) + 1;
                
                end
                
                % If we are on the last file we can increment the index
                if (j == (length(fileNames)-1))
                
                    % Increment the non-full stokes index
                    nfsIndex = fsIndex + 1;
                
                end
            
            end
            
            % We are done with this file, we can close it
            fclose('all');
            
        end
    
    end
    
end 

%--------------------------------------------------------
% THIS IS THE PROCESSING AND DISPLAY SECTION WHERE THE 
% MEAN AND STANDARD DEVIATION ARE CALCULATED AND PLOTTED
%--------------------------------------------------------

% The titles for the plots
plotTitles = {'5km' '10km' '20km' '40km' '80km' '160km'};
% The label for the x-axis
xLabel = 'Normalized x';
% The label for the y-axis
yLabel = 'Velocity (m/yr)';

% The font to use for the axis labels
axisLabelFont = 'times new roman';
% The font size to use for the axis labels
axisLabelSize = 14;
% The font weight to use for the axis labels
axisFontWeight = 'normal';

% The font size to use for tick mark labels
axisTickLabelSize = 11;
% The font weight to use for tick mark labels
axisTickFontWeight = 'normal';

% The font size to use for the title
titleFontSize = 12;
% The font weight to use for the title
titleFontWeight = 'bold';

% Loop over each domain length
for i=1:length(fileNames)
   
    % Allocate a matrix for full stokes surface velocities
    fsCurDomainVels = zeros(250, fsSizes(i));
    % Again we need a seperate index
    fsIndex = 1;
    
    % Allocate a matrix for non-full stokes surface velocities
    nfsCurDomainVels = zeros(250, fsSizes(i));
    % Again we need a seperate index
    nfsIndex = 1;
    
    % Loop over full stokes
    for j=1:length(fullStokesMask(fullStokesMask == 1))
        
        % If there are data at this entry
        if (length(fsVels{i,j}) ~= 0)
            
            % Add it to the matrix
            fsCurDomainVels(:,fsIndex) = fsVels{i,j};
            
            % Increment the index
            fsIndex = fsIndex + 1;
            
        end
        
    end
    
    % Loop over non-full stokes
    for j=1:length(fullStokesMask(fullStokesMask == 0))
        
        % If there are data at this entry
        if (length(nfsVels{i,j}) ~= 0)
            
            % Add it to the matrix
            nfsCurDomainVels(:,nfsIndex) = nfsVels{i,j};
            
            % Increment the index
            nfsIndex = nfsIndex + 1;
        end
        
    end
    
    % Calcuate the mean for each domain length (full stokes)
    fsMeanVel = mean(fsCurDomainVels, 2)';
    % Calculate the std deviation for each domain length (full stokes)
    fsStdVel = std(fsCurDomainVels, 0, 2)';
    
    % Calcuate the mean for each domain length (non-full stokes)
    nfsMeanVel = mean(nfsCurDomainVels, 2)';
    % Calculate the std deviation for each domain length (non-full stokes)
    nfsStdVel = std(nfsCurDomainVels, 0, 2)';
    
    % Tell MATLAB which subplot we are currently on
    subplot(2, 3, i);    
    
    % The fill color (rgb) for full stokes
    fsFillColor = [0.254 0.625 0.660];
    % The fill color (rgb) for non-full stokes
    nfsFillColor = [0.422 0.688 0.367];

    % Plot full stokes
    [fh1, msg1] = jbfill(newX, fsMeanVel+fsStdVel, fsMeanVel-fsStdVel, ...
                       fsFillColor, fsFillColor, 1, 0.5);
    
    % Plot non-full stokes
    [fh2, msg2] = jbfill(newX, nfsMeanVel+nfsStdVel, nfsMeanVel-nfsStdVel, ...
                       nfsFillColor, nfsFillColor, 1, 0.35);               
           
    % Turn the grid on
    grid on;
        
    % Set title properties
    th = title(plotTitles(i));
    set(th, 'FontSize', titleFontSize);
    set(th, 'FontWeight', titleFontWeight);
    
    % Set x label properties
    xlabel(xLabel);
    xh = get(gca, 'xlabel');
    set(xh, 'FontName', axisLabelFont);
    set(xh, 'FontSize', axisLabelSize);
    set(xh, 'FontWeight', axisFontWeight);

    % Set y label properties
    ylabel('Velocity (m/yr)');
    yh = get(gca, 'ylabel');
    set(yh, 'FontName', axisLabelFont);
    set(yh, 'FontSize', axisLabelSize);
    set(yh, 'FontWeight', axisFontWeight);

    % Set tick mark properties
    set(gca, 'FontSize', axisTickLabelSize);
    set(gca, 'FontWeight', axisTickFontWeight);
    
    % This turns the box around the axis on
    set(gca, 'Box', 'on' );
    
    % These lines plot the mean
    hold on;
    plot(newX, fsMeanVel, 'k');
    set(findobj(gca,'Type','line','Color',[0 0 0]),'Color',fsFillColor,'LineWidth',2);
    plot(newX, nfsMeanVel, 'b');
    set(findobj(gca,'Type','line','Color',[0 0 1]),'Color',nfsFillColor,'LineWidth',2);
    
end

% Add the legend to the final subplot (in upper left corner)
legend('FS', 'NFS', 'FS Mean', 'NFS Mean', 2);