source: issm/trunk/test/IsmipHomResults/plotting_scripts/plot_exp_a.asv@ 7888

%-----------------------------------------------------------------
% Author:       James Fishbaugh, Laura Perichon & Frank Pattyn
% Date:         06/03/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 
%--------------------------------------------------------

clear all;
close all;

% The names of the folders (also prefix for filenames)
dirs = {'aas1' 'aas2' 'ahu1' 'ahu2' 'bds1' 'cma1' 'cma2' 'dpo1' ...
        'fpa1' 'fpa2' 'fsa1' 'jvj1' 'lpe1' 'mbr1' 'mmr1' 'mtk1' ...
        'oga1' 'oso1' 'rhi1' 'rhi2' 'rhi3' 'rhi4' 'rhi5' 'spr1' ...
        'ssu1' 'tpa1' 'yko1'};
% dirs = {'fpa1'};
% 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 1 0 0 0 1 0 0 ...
                  0 1 0 0 0 0 1 0 ...
                  1 0 1 0 1 0 0 1 ...
                  1 0 1];
% fullStokesMask = [0];
% 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 = 'C:\Program Files\MATLAB\R2006a\work\ismip_all\ismip_mat\';

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

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

% 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) fileNames(j)];

        % Open the current file for reading
        load(cell2mat(fullFilePath));

        % If that file exists
        if (isnan(B(1,1))==0)

            % Store the values 
            x = B(:,1);
            y = B(:,2);
            xvel = B(:,3);
            yvel = B(:,4);
            zvel = B(:,5);
            
            % Make sure the data is normalized
            x = x./max(x);
            
            % Make sure the data is normalized
            y = y./max(y);
            
            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);

            % If the current data is from a full stokes
            if (fullStokesMask(i) == 1)
                [i j fsIndex]
                % 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);
                MaxVelfs(i,j)=max(fsVels{j,fsIndex});
                MeanVelfs(i,j)=mean(fsVels{j,fsIndex});
                
                % 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)))
                    
                    % Increment the full stokes index
                    fsIndex = fsIndex + 1;
                
                end
                
            % Else its non-full stokes
            else
                [i j nfsIndex]
                % 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);
                MaxVelnfs(i,j)=max(nfsVels{j,nfsIndex});
                MeanVelnfs(i,j)=mean(nfsVels{j,nfsIndex});
                
                % Keep track of the number of entries in the array
                if (length(nfsVels{j, nfsIndex}) > 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)))
                
                    % Increment the non-full stokes index
                    nfsIndex = nfsIndex + 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 = 'normal';
% The font size to use for the axis labels
axisLabelSize = 7;
% The font weight to use for the axis labels
axisFontWeight = 'normal';

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

% The font size to use for the title
titleFontSize = 10;
% 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', 'Location','West');

set(gcf, 'paperpositionmode','manual','paperorientation', ...
        'portrait','papertype','A4','paperunits', 'centimeters', ...
        'paperposition',[1 4 21 13]);
saveas(gcf, 'expa.eps', 'psc2');

% Analysis of maximum velocities (F. Pattyn)
% Full Stokes stored in MaxVelfs (else MaxVelnfs)
% cols = L, rows = participants, 0 = no result

for i=1:length(MaxVelfs)
    for j=1:size(MaxVelfs,2)
        if MaxVelfs(i,j)==0
            MaxVelfs(i,j)=NaN;
        end
    end
end
for j=1:size(MaxVelfs,2)
    MeanMaxfs(j)=mean(MaxVelfs(~isnan(MaxVelfs(:,j)),j));
    StdDevfs(j)=std(MaxVelfs(~isnan(MaxVelfs(:,j)),j));
    Numfs(j)=sum(~isnan(MaxVelfs(:,j)));
end
for i=1:length(MaxVelnfs)
    for j=1:size(MaxVelnfs,2)
        if MaxVelnfs(i,j)==0
            MaxVelnfs(i,j)=NaN;
        end
    end
end
for j=1:size(MaxVelnfs,2)
    MeanMaxnfs(j)=mean(MaxVelnfs(~isnan(MaxVelnfs(:,j)),j));
    StdDevnfs(j)=std(MaxVelnfs(~isnan(MaxVelnfs(:,j)),j));
    Numnfs(j)=sum(~isnan(MaxVelnfs(:,j)));
end

% print in format for LaTeX file
fprintf('NFS & ');
for j=1:size(MaxVelnfs,2)
    fprintf('%.2f & %.2f & %d & ', MeanMaxnfs(j), StdDevnfs(j), Numnfs(j));
end
fprintf('\\\\ \n');
fprintf('FS & ');
for j=1:size(MaxVelfs,2)
    fprintf('%.2f & %.2f & %d & ', MeanMaxfs(j), StdDevfs(j), Numfs(j));
end
fprintf('\\\\ \n');

% Analysis of mean velocities (F. Pattyn)
% Full Stokes stored in MaxVelfs (else MaxVelnfs)
% cols = L, rows = participants, 0 = no result

for i=1:length(MaxVelfs)
    for j=1:size(MaxVelfs,2)
        if MaxVelfs(i,j)==0
            MaxVelfs(i,j)=NaN;
        end
    end
end
for j=1:size(MaxVelfs,2)
    MeanMaxfs(j)=mean(MaxVelfs(~isnan(MaxVelfs(:,j)),j));
    StdDevfs(j)=std(MaxVelfs(~isnan(MaxVelfs(:,j)),j));
    Numfs(j)=sum(~isnan(MaxVelfs(:,j)));
end
for i=1:length(MaxVelnfs)
    for j=1:size(MaxVelnfs,2)
        if MaxVelnfs(i,j)==0
            MaxVelnfs(i,j)=NaN;
        end
    end
end
for j=1:size(MaxVelnfs,2)
    MeanMaxnfs(j)=mean(MaxVelnfs(~isnan(MaxVelnfs(:,j)),j));
    StdDevnfs(j)=std(MaxVelnfs(~isnan(MaxVelnfs(:,j)),j));
    Numnfs(j)=sum(~isnan(MaxVelnfs(:,j)));
end

% print in format for LaTeX file
fprintf('NFS & ');
for j=1:size(MaxVelnfs,2)
    fprintf('%.2f & %.2f & %d & ', MeanMaxnfs(j), StdDevnfs(j), Numnfs(j));
end
fprintf('\\\\ \n');
fprintf('FS & ');
for j=1:size(MaxVelfs,2)
    fprintf('%.2f & %.2f & %d & ', MeanMaxfs(j), StdDevfs(j), Numfs(j));
end
fprintf('\\\\ \n');