Index: /issm/trunk/test/IsmipHomResults/README
===================================================================
--- /issm/trunk/test/IsmipHomResults/README	(revision 7888)
+++ /issm/trunk/test/IsmipHomResults/README	(revision 7888)
@@ -0,0 +1,4 @@
+These are the results from all participating model of the ISMIP-HOM benchmarks
+They have been provided by Frank Pattyn (pers. comm.)
+
+Some files have been modified from their original version
Index: /issm/trunk/test/IsmipHomResults/plotting_scripts/inpaint_nans.m
===================================================================
--- /issm/trunk/test/IsmipHomResults/plotting_scripts/inpaint_nans.m	(revision 7888)
+++ /issm/trunk/test/IsmipHomResults/plotting_scripts/inpaint_nans.m	(revision 7888)
@@ -0,0 +1,533 @@
+function B=inpaint_nans(A,method)
+% INPAINT_NANS: in-paints over nans in an array
+% usage: B=INPAINT_NANS(A)          % default method
+% usage: B=INPAINT_NANS(A,method)   % specify method used
+%
+% Solves approximation to one of several pdes to
+% interpolate and extrapolate holes in an array
+%
+% arguments (input):
+%   A - nxm array with some NaNs to be filled in
+%
+%   method - (OPTIONAL) scalar numeric flag - specifies
+%       which approach (or physical metaphor to use
+%       for the interpolation.) All methods are capable
+%       of extrapolation, some are better than others.
+%       There are also speed differences, as well as
+%       accuracy differences for smooth surfaces.
+%
+%       methods {0,1,2} use a simple plate metaphor.
+%       method  3 uses a better plate equation,
+%                 but may be much slower and uses
+%                 more memory.
+%       method  4 uses a spring metaphor.
+%       method  5 is an 8 neighbor average, with no
+%                 rationale behind it compared to the
+%                 other methods. I do not recommend
+%                 its use.
+%
+%       method == 0 --> (DEFAULT) see method 1, but
+%         this method does not build as large of a
+%         linear system in the case of only a few
+%         NaNs in a large array.
+%         Extrapolation behavior is linear.
+%         
+%       method == 1 --> simple approach, applies del^2
+%         over the entire array, then drops those parts
+%         of the array which do not have any contact with
+%         NaNs. Uses a least squares approach, but it
+%         does not modify known values.
+%         In the case of small arrays, this method is
+%         quite fast as it does very little extra work.
+%         Extrapolation behavior is linear.
+%         
+%       method == 2 --> uses del^2, but solving a direct
+%         linear system of equations for nan elements.
+%         This method will be the fastest possible for
+%         large systems since it uses the sparsest
+%         possible system of equations. Not a least
+%         squares approach, so it may be least robust
+%         to noise on the boundaries of any holes.
+%         This method will also be least able to
+%         interpolate accurately for smooth surfaces.
+%         Extrapolation behavior is linear.
+%         
+%       method == 3 --+ See method 0, but uses del^4 for
+%         the interpolating operator. This may result
+%         in more accurate interpolations, at some cost
+%         in speed.
+%         
+%       method == 4 --+ Uses a spring metaphor. Assumes
+%         springs (with a nominal length of zero)
+%         connect each node with every neighbor
+%         (horizontally, vertically and diagonally)
+%         Since each node tries to be like its neighbors,
+%         extrapolation is as a constant function where
+%         this is consistent with the neighboring nodes.
+%
+%       method == 5 --+ See method 2, but use an average
+%         of the 8 nearest neighbors to any element.
+%         This method is NOT recommended for use.
+%
+%
+% arguments (output):
+%   B - nxm array with NaNs replaced
+%
+%
+% Example:
+%  [x,y] = meshgrid(0:.01:1);
+%  z0 = exp(x+y);
+%  znan = z0;
+%  znan(20:50,40:70) = NaN;
+%  znan(30:90,5:10) = NaN;
+%  znan(70:75,40:90) = NaN;
+%
+%  z = inpaint_nans(znan);
+%
+%
+% See also: griddata, interp1
+%
+% Author: John D'Errico
+% e-mail address: woodchips@rochester.rr.com
+% Release: 2
+% Release date: 4/15/06
+
+
+% I always need to know which elements are NaN,
+% and what size the array is for any method
+[n,m]=size(A);
+A=A(:);
+nm=n*m;
+k=isnan(A(:));
+
+% list the nodes which are known, and which will
+% be interpolated
+nan_list=find(k);
+known_list=find(~k);
+
+% how many nans overall
+nan_count=length(nan_list);
+
+% convert NaN indices to (r,c) form
+% nan_list==find(k) are the unrolled (linear) indices
+% (row,column) form
+[nr,nc]=ind2sub([n,m],nan_list);
+
+% both forms of index in one array:
+% column 1 == unrolled index
+% column 2 == row index
+% column 3 == column index
+nan_list=[nan_list,nr,nc];
+
+% supply default method
+if (nargin<2) || isempty(method)
+  method = 0;
+elseif ~ismember(method,0:5)
+  error 'If supplied, method must be one of: {0,1,2,3,4,5}.'
+end
+
+% for different methods
+switch method
+ case 0
+  % The same as method == 1, except only work on those
+  % elements which are NaN, or at least touch a NaN.
+  
+  % horizontal and vertical neighbors only
+  talks_to = [-1 0;0 -1;1 0;0 1];
+  neighbors_list=identify_neighbors(n,m,nan_list,talks_to);
+  
+  % list of all nodes we have identified
+  all_list=[nan_list;neighbors_list];
+  
+  % generate sparse array with second partials on row
+  % variable for each element in either list, but only
+  % for those nodes which have a row index > 1 or < n
+  L = find((all_list(:,2) > 1) & (all_list(:,2) < n)); 
+  nl=length(L);
+  if nl>0
+    fda=sparse(repmat(all_list(L,1),1,3), ...
+      repmat(all_list(L,1),1,3)+repmat([-1 0 1],nl,1), ...
+      repmat([1 -2 1],nl,1),nm,nm);
+  else
+    fda=spalloc(n*m,n*m,size(all_list,1)*5);
+  end
+  
+  % 2nd partials on column index
+  L = find((all_list(:,3) > 1) & (all_list(:,3) < m)); 
+  nl=length(L);
+  if nl>0
+    fda=fda+sparse(repmat(all_list(L,1),1,3), ...
+      repmat(all_list(L,1),1,3)+repmat([-n 0 n],nl,1), ...
+      repmat([1 -2 1],nl,1),nm,nm);
+  end
+  
+  % eliminate knowns
+  rhs=-fda(:,known_list)*A(known_list);
+  k=find(any(fda(:,nan_list(:,1)),2));
+  
+  % and solve...
+  B=A;
+  B(nan_list(:,1))=fda(k,nan_list(:,1))\rhs(k);
+  
+ case 1
+  % least squares approach with del^2. Build system
+  % for every array element as an unknown, and then
+  % eliminate those which are knowns.
+
+  % Build sparse matrix approximating del^2 for
+  % every element in A.
+  % Compute finite difference for second partials
+  % on row variable first
+  [i,j]=ndgrid(2:(n-1),1:m);
+  ind=i(:)+(j(:)-1)*n;
+  np=(n-2)*m;
+  fda=sparse(repmat(ind,1,3),[ind-1,ind,ind+1], ...
+      repmat([1 -2 1],np,1),n*m,n*m);
+  
+  % now second partials on column variable
+  [i,j]=ndgrid(1:n,2:(m-1));
+  ind=i(:)+(j(:)-1)*n;
+  np=n*(m-2);
+  fda=fda+sparse(repmat(ind,1,3),[ind-n,ind,ind+n], ...
+      repmat([1 -2 1],np,1),nm,nm);
+  
+  % eliminate knowns
+  rhs=-fda(:,known_list)*A(known_list);
+  k=find(any(fda(:,nan_list),2));
+  
+  % and solve...
+  B=A;
+  B(nan_list(:,1))=fda(k,nan_list(:,1))\rhs(k);
+  
+ case 2
+  % Direct solve for del^2 BVP across holes
+
+  % generate sparse array with second partials on row
+  % variable for each nan element, only for those nodes
+  % which have a row index > 1 or < n
+  L = find((nan_list(:,2) > 1) & (nan_list(:,2) < n)); 
+  nl=length(L);
+  if nl>0
+    fda=sparse(repmat(nan_list(L,1),1,3), ...
+      repmat(nan_list(L,1),1,3)+repmat([-1 0 1],nl,1), ...
+      repmat([1 -2 1],nl,1),n*m,n*m);
+  else
+    fda=spalloc(n*m,n*m,size(nan_list,1)*5);
+  end
+  
+  % 2nd partials on column index
+  L = find((nan_list(:,3) > 1) & (nan_list(:,3) < m)); 
+  nl=length(L);
+  if nl>0
+    fda=fda+sparse(repmat(nan_list(L,1),1,3), ...
+      repmat(nan_list(L,1),1,3)+repmat([-n 0 n],nl,1), ...
+      repmat([1 -2 1],nl,1),n*m,n*m);
+  end
+  
+  % fix boundary conditions at extreme corners
+  % of the array in case there were nans there
+  if ismember(1,nan_list(:,1))
+    fda(1,[1 2 n+1])=[-2 1 1];
+  end
+  if ismember(n,nan_list(:,1))
+    fda(n,[n, n-1,n+n])=[-2 1 1];
+  end
+  if ismember(nm-n+1,nan_list(:,1))
+    fda(nm-n+1,[nm-n+1,nm-n+2,nm-n])=[-2 1 1];
+  end
+  if ismember(nm,nan_list(:,1))
+    fda(nm,[nm,nm-1,nm-n])=[-2 1 1];
+  end
+  
+  % eliminate knowns
+  rhs=-fda(:,known_list)*A(known_list);
+  
+  % and solve...
+  B=A;
+  k=nan_list(:,1);
+  B(k)=fda(k,k)\rhs(k);
+  
+ case 3
+  % The same as method == 0, except uses del^4 as the
+  % interpolating operator.
+  
+  % del^4 template of neighbors
+  talks_to = [-2 0;-1 -1;-1 0;-1 1;0 -2;0 -1; ...
+      0 1;0 2;1 -1;1 0;1 1;2 0];
+  neighbors_list=identify_neighbors(n,m,nan_list,talks_to);
+  
+  % list of all nodes we have identified
+  all_list=[nan_list;neighbors_list];
+  
+  % generate sparse array with del^4, but only
+  % for those nodes which have a row & column index
+  % >= 3 or <= n-2
+  L = find( (all_list(:,2) >= 3) & ...
+            (all_list(:,2) <= (n-2)) & ...
+            (all_list(:,3) >= 3) & ...
+            (all_list(:,3) <= (m-2)));
+  nl=length(L);
+  if nl>0
+    % do the entire template at once
+    fda=sparse(repmat(all_list(L,1),1,13), ...
+        repmat(all_list(L,1),1,13) + ...
+        repmat([-2*n,-n-1,-n,-n+1,-2,-1,0,1,2,n-1,n,n+1,2*n],nl,1), ...
+        repmat([1 2 -8 2 1 -8 20 -8 1 2 -8 2 1],nl,1),nm,nm);
+  else
+    fda=spalloc(n*m,n*m,size(all_list,1)*5);
+  end
+  
+  % on the boundaries, reduce the order around the edges
+  L = find((((all_list(:,2) == 2) | ...
+             (all_list(:,2) == (n-1))) & ...
+            (all_list(:,3) >= 2) & ...
+            (all_list(:,3) <= (m-1))) | ...
+           (((all_list(:,3) == 2) | ...
+             (all_list(:,3) == (m-1))) & ...
+            (all_list(:,2) >= 2) & ...
+            (all_list(:,2) <= (n-1))));
+  nl=length(L);
+  if nl>0
+    fda=fda+sparse(repmat(all_list(L,1),1,5), ...
+      repmat(all_list(L,1),1,5) + ...
+        repmat([-n,-1,0,+1,n],nl,1), ...
+      repmat([1 1 -4 1 1],nl,1),nm,nm);
+  end
+  
+  L = find( ((all_list(:,2) == 1) | ...
+             (all_list(:,2) == n)) & ...
+            (all_list(:,3) >= 2) & ...
+            (all_list(:,3) <= (m-1)));
+  nl=length(L);
+  if nl>0
+    fda=fda+sparse(repmat(all_list(L,1),1,3), ...
+      repmat(all_list(L,1),1,3) + ...
+        repmat([-n,0,n],nl,1), ...
+      repmat([1 -2 1],nl,1),nm,nm);
+  end
+  
+  L = find( ((all_list(:,3) == 1) | ...
+             (all_list(:,3) == m)) & ...
+            (all_list(:,2) >= 2) & ...
+            (all_list(:,2) <= (n-1)));
+  nl=length(L);
+  if nl>0
+    fda=fda+sparse(repmat(all_list(L,1),1,3), ...
+      repmat(all_list(L,1),1,3) + ...
+        repmat([-1,0,1],nl,1), ...
+      repmat([1 -2 1],nl,1),nm,nm);
+  end
+  
+  % eliminate knowns
+  rhs=-fda(:,known_list)*A(known_list);
+  k=find(any(fda(:,nan_list(:,1)),2));
+  
+  % and solve...
+  B=A;
+  B(nan_list(:,1))=fda(k,nan_list(:,1))\rhs(k);
+  
+ case 4
+  % Spring analogy
+  % interpolating operator.
+  
+  % list of all springs between a node and a horizontal
+  % or vertical neighbor
+  hv_list=[-1 -1 0;1 1 0;-n 0 -1;n 0 1];
+  hv_springs=[];
+  for i=1:4
+    hvs=nan_list+repmat(hv_list(i,:),nan_count,1);
+    k=(hvs(:,2)>=1) & (hvs(:,2)<=n) & (hvs(:,3)>=1) & (hvs(:,3)<=m);
+    hv_springs=[hv_springs;[nan_list(k,1),hvs(k,1)]];
+  end
+
+  % delete replicate springs
+  hv_springs=unique(sort(hv_springs,2),'rows');
+  
+  % build sparse matrix of connections, springs
+  % connecting diagonal neighbors are weaker than
+  % the horizontal and vertical springs
+  nhv=size(hv_springs,1);
+  springs=sparse(repmat((1:nhv)',1,2),hv_springs, ...
+     repmat([1 -1],nhv,1),nhv,nm);
+  
+  % eliminate knowns
+  rhs=-springs(:,known_list)*A(known_list);
+  
+  % and solve...
+  B=A;
+  B(nan_list(:,1))=springs(:,nan_list(:,1))\rhs;
+  
+ case 5
+  % Average of 8 nearest neighbors
+  
+  % generate sparse array to average 8 nearest neighbors
+  % for each nan element, be careful around edges
+  fda=spalloc(n*m,n*m,size(nan_list,1)*9);
+  
+  % -1,-1
+  L = find((nan_list(:,2) > 1) & (nan_list(:,3) > 1)); 
+  nl=length(L);
+  if nl>0
+    fda=fda+sparse(repmat(nan_list(L,1),1,2), ...
+      repmat(nan_list(L,1),1,2)+repmat([-n-1, 0],nl,1), ...
+      repmat([1 -1],nl,1),n*m,n*m);
+  end
+  
+  % 0,-1
+  L = find(nan_list(:,3) > 1);
+  nl=length(L);
+  if nl>0
+    fda=fda+sparse(repmat(nan_list(L,1),1,2), ...
+      repmat(nan_list(L,1),1,2)+repmat([-n, 0],nl,1), ...
+      repmat([1 -1],nl,1),n*m,n*m);
+  end
+
+  % +1,-1
+  L = find((nan_list(:,2) < n) & (nan_list(:,3) > 1));
+  nl=length(L);
+  if nl>0
+    fda=fda+sparse(repmat(nan_list(L,1),1,2), ...
+      repmat(nan_list(L,1),1,2)+repmat([-n+1, 0],nl,1), ...
+      repmat([1 -1],nl,1),n*m,n*m);
+  end
+
+  % -1,0
+  L = find(nan_list(:,2) > 1);
+  nl=length(L);
+  if nl>0
+    fda=fda+sparse(repmat(nan_list(L,1),1,2), ...
+      repmat(nan_list(L,1),1,2)+repmat([-1, 0],nl,1), ...
+      repmat([1 -1],nl,1),n*m,n*m);
+  end
+
+  % +1,0
+  L = find(nan_list(:,2) < n);
+  nl=length(L);
+  if nl>0
+    fda=fda+sparse(repmat(nan_list(L,1),1,2), ...
+      repmat(nan_list(L,1),1,2)+repmat([1, 0],nl,1), ...
+      repmat([1 -1],nl,1),n*m,n*m);
+  end
+
+  % -1,+1
+  L = find((nan_list(:,2) > 1) & (nan_list(:,3) < m)); 
+  nl=length(L);
+  if nl>0
+    fda=fda+sparse(repmat(nan_list(L,1),1,2), ...
+      repmat(nan_list(L,1),1,2)+repmat([n-1, 0],nl,1), ...
+      repmat([1 -1],nl,1),n*m,n*m);
+  end
+  
+  % 0,+1
+  L = find(nan_list(:,3) < m);
+  nl=length(L);
+  if nl>0
+    fda=fda+sparse(repmat(nan_list(L,1),1,2), ...
+      repmat(nan_list(L,1),1,2)+repmat([n, 0],nl,1), ...
+      repmat([1 -1],nl,1),n*m,n*m);
+  end
+
+  % +1,+1
+  L = find((nan_list(:,2) < n) & (nan_list(:,3) < m));
+  nl=length(L);
+  if nl>0
+    fda=fda+sparse(repmat(nan_list(L,1),1,2), ...
+      repmat(nan_list(L,1),1,2)+repmat([n+1, 0],nl,1), ...
+      repmat([1 -1],nl,1),n*m,n*m);
+  end
+  
+  % eliminate knowns
+  rhs=-fda(:,known_list)*A(known_list);
+  
+  % and solve...
+  B=A;
+  k=nan_list(:,1);
+  B(k)=fda(k,k)\rhs(k);
+  
+end
+
+% all done, make sure that B is the same shape as
+% A was when we came in.
+B=reshape(B,n,m);
+
+
+% ====================================================
+%      end of main function
+% ====================================================
+% ====================================================
+%      begin subfunctions
+% ====================================================
+function neighbors_list=identify_neighbors(n,m,nan_list,talks_to)
+% identify_neighbors: identifies all the neighbors of
+%   those nodes in nan_list, not including the nans
+%   themselves
+%
+% arguments (input):
+%  n,m - scalar - [n,m]=size(A), where A is the
+%      array to be interpolated
+%  nan_list - array - list of every nan element in A
+%      nan_list(i,1) == linear index of i'th nan element
+%      nan_list(i,2) == row index of i'th nan element
+%      nan_list(i,3) == column index of i'th nan element
+%  talks_to - px2 array - defines which nodes communicate
+%      with each other, i.e., which nodes are neighbors.
+%
+%      talks_to(i,1) - defines the offset in the row
+%                      dimension of a neighbor
+%      talks_to(i,2) - defines the offset in the column
+%                      dimension of a neighbor
+%      
+%      For example, talks_to = [-1 0;0 -1;1 0;0 1]
+%      means that each node talks only to its immediate
+%      neighbors horizontally and vertically.
+% 
+% arguments(output):
+%  neighbors_list - array - list of all neighbors of
+%      all the nodes in nan_list
+
+if ~isempty(nan_list)
+  % use the definition of a neighbor in talks_to
+  nan_count=size(nan_list,1);
+  talk_count=size(talks_to,1);
+  
+  nn=zeros(nan_count*talk_count,2);
+  j=[1,nan_count];
+  for i=1:talk_count
+    nn(j(1):j(2),:)=nan_list(:,2:3) + ...
+        repmat(talks_to(i,:),nan_count,1);
+    j=j+nan_count;
+  end
+  
+  % drop those nodes which fall outside the bounds of the
+  % original array
+  L = (nn(:,1)<1)|(nn(:,1)>n)|(nn(:,2)<1)|(nn(:,2)>m); 
+  nn(L,:)=[];
+  
+  % form the same format 3 column array as nan_list
+  neighbors_list=[sub2ind([n,m],nn(:,1),nn(:,2)),nn];
+  
+  % delete replicates in the neighbors list
+  neighbors_list=unique(neighbors_list,'rows');
+  
+  % and delete those which are also in the list of NaNs.
+  neighbors_list=setdiff(neighbors_list,nan_list,'rows');
+  
+else
+  neighbors_list=[];
+end
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
Index: /issm/trunk/test/IsmipHomResults/plotting_scripts/jbfill.m
===================================================================
--- /issm/trunk/test/IsmipHomResults/plotting_scripts/jbfill.m	(revision 7888)
+++ /issm/trunk/test/IsmipHomResults/plotting_scripts/jbfill.m	(revision 7888)
@@ -0,0 +1,44 @@
+function[fillhandle,msg]=jbfill(xpoints,upper,lower,color,edge,add,transparency)
+%USAGE: [fillhandle,msg]=jbfill(xpoints,upper,lower,color,edge,add,transparency)
+%This function will fill a region with a color between the two vectors provided
+%using the Matlab fill command.
+%
+%fillhandle is the returned handle to the filled region in the plot.
+%xpoints= The horizontal data points (ie frequencies). Note length(Upper)
+%         must equal Length(lower)and must equal length(xpoints)!
+%upper = the upper curve values (data can be less than lower)
+%lower = the lower curve values (data can be more than upper)
+%color = the color of the filled area 
+%edge  = the color around the edge of the filled area
+%add   = a flag to add to the current plot or make a new one.
+%transparency is a value ranging from 1 for opaque to 0 for invisible for
+%the filled color only.
+%
+%John A. Bockstege November 2006;
+%Example:
+%     a=rand(1,20);%Vector of random data
+%     b=a+2*rand(1,20);%2nd vector of data points;
+%     x=1:20;%horizontal vector
+%     [ph,msg]=jbfill(x,a,b,rand(1,3),rand(1,3),0,rand(1,1))
+%     grid on
+%     legend('Datr')
+if nargin<7;transparency=.5;end %default is to have a transparency of .5
+if nargin<6;add=1;end     %default is to add to current plot
+if nargin<5;edge='k';end  %dfault edge color is black
+if nargin<4;color='b';end %default color is blue
+
+if length(upper)==length(lower) && length(lower)==length(xpoints)
+    msg='';
+    filled=[upper,fliplr(lower)];
+    xpoints=[xpoints,fliplr(xpoints)];
+    if add
+        hold on
+    end
+    fillhandle=fill(xpoints,filled,color);%plot the data
+    set(fillhandle,'EdgeColor',edge,'FaceAlpha',transparency,'EdgeAlpha',transparency);%set edge color
+    if add
+        hold off
+    end
+else
+    msg='Error: Must use the same number of points in each vector';
+end
Index: /issm/trunk/test/IsmipHomResults/plotting_scripts/normtest.m
===================================================================
--- /issm/trunk/test/IsmipHomResults/plotting_scripts/normtest.m	(revision 7888)
+++ /issm/trunk/test/IsmipHomResults/plotting_scripts/normtest.m	(revision 7888)
@@ -0,0 +1,17 @@
+clear all;
+close all;
+
+tot=20;
+
+dx=zeros(tot);
+norm=zeros(tot);
+dx(1)=1.;
+for i=1:tot
+    if i>1
+        dx(i)=dx(i-1)/2.;
+    end
+    x=1:dx(i):10;
+    norm(i)=sum(sqrt(sin(x).^2))/length(x);
+end
+
+semilogx(dx,norm,'-o');
Index: /issm/trunk/test/IsmipHomResults/plotting_scripts/plot_exp_a.asv
===================================================================
--- /issm/trunk/test/IsmipHomResults/plotting_scripts/plot_exp_a.asv	(revision 7888)
+++ /issm/trunk/test/IsmipHomResults/plotting_scripts/plot_exp_a.asv	(revision 7888)
@@ -0,0 +1,389 @@
+%-----------------------------------------------------------------
+% 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');
+
Index: /issm/trunk/test/IsmipHomResults/plotting_scripts/plot_exp_a.m
===================================================================
--- /issm/trunk/test/IsmipHomResults/plotting_scripts/plot_exp_a.m	(revision 7888)
+++ /issm/trunk/test/IsmipHomResults/plotting_scripts/plot_exp_a.m	(revision 7888)
@@ -0,0 +1,389 @@
+%-----------------------------------------------------------------
+% 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 = {'dpo1' 'lpe1' 'rhi4' 'rhi5'};
+% 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 1 0 0 1 0 ...
+                  1 0 1 0 1 0 0 1 ...
+                  1 0 1];
+% fullStokesMask = [0 0 0 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:\Documents and Settings\frankp\My Documents\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 a^{-1})';
+
+% 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 a^{-1})');
+    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 surface velocities (F. Pattyn)
+% Full Stokes stored in MaxVelfs (else MaxVelnfs)
+% cols = L, rows = participants, 0 = no result
+
+for i=1:length(MeanVelfs)
+    for j=1:size(MeanVelfs,2)
+        if MeanVelfs(i,j)==0
+            MeanVelfs(i,j)=NaN;
+        end
+    end
+end
+for j=1:size(MaxVelfs,2)
+    MeanMaxfs(j)=mean(MeanVelfs(~isnan(MeanVelfs(:,j)),j));
+    StdDevfs(j)=std(MeanVelfs(~isnan(MeanVelfs(:,j)),j));
+    Numfs(j)=sum(~isnan(MeanVelfs(:,j)));
+end
+for i=1:length(MeanVelnfs)
+    for j=1:size(MeanVelnfs,2)
+        if MeanVelnfs(i,j)==0
+            MeanVelnfs(i,j)=NaN;
+        end
+    end
+end
+for j=1:size(MeanVelnfs,2)
+    MeanMaxnfs(j)=mean(MeanVelnfs(~isnan(MeanVelnfs(:,j)),j));
+    StdDevnfs(j)=std(MeanVelnfs(~isnan(MeanVelnfs(:,j)),j));
+    Numnfs(j)=sum(~isnan(MeanVelnfs(:,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');
+
Index: /issm/trunk/test/IsmipHomResults/plotting_scripts/plot_exp_b.asv
===================================================================
--- /issm/trunk/test/IsmipHomResults/plotting_scripts/plot_exp_b.asv	(revision 7888)
+++ /issm/trunk/test/IsmipHomResults/plotting_scripts/plot_exp_b.asv	(revision 7888)
@@ -0,0 +1,294 @@
+%-----------------------------------------------------------------
+% 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' 'ahu1' 'ahu2' 'bds1' 'cma1' 'cma2' 'dpo1' ...
+        'fpa1' 'fpa2' 'fsa1' 'jvj1' 'lpe1' 'mbr1' 'mmr1' 'mtk1' ...
+        'oga1' 'rhi1' 'rhi2' 'rhi3' 'rhi4' 'rhi5' 'spr1' ...
+        'ssu1' 'tpa1' 'yko1'};
+dirs = {'fsa1'};
+% 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 1 0 0 ...
+                  0 1 0 0 0 0 1 0 ...
+                  1 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 = 'b';
+
+% 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);
+
+% 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);
+            xvel = B(:,2);
+            zvel = B(:,3);
+            
+            % Make sure the data is normalized
+            x = x./max(x);
+             
+            x = inpaint_nans(x, 4);
+            xvel = inpaint_nans(xvel, 4);
+            zvel = inpaint_nans(zvel, 4);
+            
+            % Calculates the magnitude of the surface velocity
+            
+            surfVel = sqrt(xvel.^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} = interp1(x, surfVel, newX);
+                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)))
+                    
+                    % 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} = interp1(x, surfVel, newX);
+                nfsVels{j, nfsIndex} = inpaint_nans(nfsVels{j, nfsIndex}, 4);
+                
+                % 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 = 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);
+
+set(gcf, 'paperpositionmode','manual','paperorientation', ...
+        'landscape','papertype','A4','paperunits', 'centimeters', ...
+        'paperposition',[1 4 25 15]);
+saveas(gcf, 'expd.eps', 'psc2');
Index: /issm/trunk/test/IsmipHomResults/plotting_scripts/plot_exp_b.m
===================================================================
--- /issm/trunk/test/IsmipHomResults/plotting_scripts/plot_exp_b.m	(revision 7888)
+++ /issm/trunk/test/IsmipHomResults/plotting_scripts/plot_exp_b.m	(revision 7888)
@@ -0,0 +1,381 @@
+%-----------------------------------------------------------------
+% 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 = {'dpo1' 'lpe1' 'rhi4' 'rhi5'};
+% 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 1 0 0 1 0 ...
+                  1 0 1 0 1 0 0 1 ...
+                  1 0 1];
+% fullStokesMask = [0 0 0 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 = 'b';
+
+% 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);
+
+% 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);
+            xvel = B(:,2);
+            zvel = B(:,3);
+            
+            % Make sure the data is normalized
+            x = x./max(x);
+             
+            x = inpaint_nans(x, 4);
+            xvel = inpaint_nans(xvel, 4);
+            zvel = inpaint_nans(zvel, 4);
+            
+            % Calculates the magnitude of the surface velocity
+            
+            surfVel = sqrt(xvel.^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} = interp1(x, surfVel, newX);
+                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} = interp1(x, surfVel, newX);
+                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 a^{-1})';
+
+% 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 a^{-1})');
+    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, 'expb.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 surface velocities (F. Pattyn)
+% Full Stokes stored in MaxVelfs (else MaxVelnfs)
+% cols = L, rows = participants, 0 = no result
+
+for i=1:length(MeanVelfs)
+    for j=1:size(MeanVelfs,2)
+        if MeanVelfs(i,j)==0
+            MeanVelfs(i,j)=NaN;
+        end
+    end
+end
+for j=1:size(MaxVelfs,2)
+    MeanMaxfs(j)=mean(MeanVelfs(~isnan(MeanVelfs(:,j)),j));
+    StdDevfs(j)=std(MeanVelfs(~isnan(MeanVelfs(:,j)),j));
+    Numfs(j)=sum(~isnan(MeanVelfs(:,j)));
+end
+for i=1:length(MeanVelnfs)
+    for j=1:size(MeanVelnfs,2)
+        if MeanVelnfs(i,j)==0
+            MeanVelnfs(i,j)=NaN;
+        end
+    end
+end
+for j=1:size(MeanVelnfs,2)
+    MeanMaxnfs(j)=mean(MeanVelnfs(~isnan(MeanVelnfs(:,j)),j));
+    StdDevnfs(j)=std(MeanVelnfs(~isnan(MeanVelnfs(:,j)),j));
+    Numnfs(j)=sum(~isnan(MeanVelnfs(:,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');
+
Index: /issm/trunk/test/IsmipHomResults/plotting_scripts/plot_exp_c.m
===================================================================
--- /issm/trunk/test/IsmipHomResults/plotting_scripts/plot_exp_c.m	(revision 7888)
+++ /issm/trunk/test/IsmipHomResults/plotting_scripts/plot_exp_c.m	(revision 7888)
@@ -0,0 +1,389 @@
+%-----------------------------------------------------------------
+% 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 = {'dpo1' 'lpe1' 'rhi4' 'rhi5'};
+% 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 1 0 0 1 0 ...
+                  1 0 1 0 1 0 0 1 ...
+                  1 0 1];
+% fullStokesMask = [0 0 0 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:\Documents and Settings\frankp\My Documents\work\ismip_all\ismip_mat\';
+
+% The experiment we're interested in
+exp = 'c';
+
+% 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 a^{-1})';
+
+% 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 a^{-1})');
+    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, 'expc.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 surface velocities (F. Pattyn)
+% Full Stokes stored in MaxVelfs (else MaxVelnfs)
+% cols = L, rows = participants, 0 = no result
+
+for i=1:length(MeanVelfs)
+    for j=1:size(MeanVelfs,2)
+        if MeanVelfs(i,j)==0
+            MeanVelfs(i,j)=NaN;
+        end
+    end
+end
+for j=1:size(MaxVelfs,2)
+    MeanMaxfs(j)=mean(MeanVelfs(~isnan(MeanVelfs(:,j)),j));
+    StdDevfs(j)=std(MeanVelfs(~isnan(MeanVelfs(:,j)),j));
+    Numfs(j)=sum(~isnan(MeanVelfs(:,j)));
+end
+for i=1:length(MeanVelnfs)
+    for j=1:size(MeanVelnfs,2)
+        if MeanVelnfs(i,j)==0
+            MeanVelnfs(i,j)=NaN;
+        end
+    end
+end
+for j=1:size(MeanVelnfs,2)
+    MeanMaxnfs(j)=mean(MeanVelnfs(~isnan(MeanVelnfs(:,j)),j));
+    StdDevnfs(j)=std(MeanVelnfs(~isnan(MeanVelnfs(:,j)),j));
+    Numnfs(j)=sum(~isnan(MeanVelnfs(:,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');
+
Index: /issm/trunk/test/IsmipHomResults/plotting_scripts/plot_exp_d.asv
===================================================================
--- /issm/trunk/test/IsmipHomResults/plotting_scripts/plot_exp_d.asv	(revision 7888)
+++ /issm/trunk/test/IsmipHomResults/plotting_scripts/plot_exp_d.asv	(revision 7888)
@@ -0,0 +1,381 @@
+%-----------------------------------------------------------------
+% 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' 'lpe1' 'mbr1' 'mmr1' 'mtk1' ...
+        'oga1' 'oso1' 'rhi1' 'rhi2' 'rhi3' 'rhi4' 'rhi5' 'spr1' ...
+        'ssu1' 'tpa1' 'yko1'};
+dirs = {'dpo1' 'lpe1' 'rhi4' 'rhi5'};
+% 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 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 = 'd';
+
+% 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);
+
+% 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);
+            xvel = B(:,2);
+            zvel = B(:,3);
+            
+            % Make sure the data is normalized
+            x = x./max(x);
+             
+            x = inpaint_nans(x, 4);
+            xvel = inpaint_nans(xvel, 4);
+            zvel = inpaint_nans(zvel, 4);
+            
+            % Calculates the magnitude of the surface velocity
+            
+            surfVel = sqrt(xvel.^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} = interp1(x, surfVel, newX);
+                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} = interp1(x, surfVel, newX);
+                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 a^{-1})';
+
+% 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 a^{-1})');
+    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, 'expd.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 surface velocities (F. Pattyn)
+% Full Stokes stored in MaxVelfs (else MaxVelnfs)
+% cols = L, rows = participants, 0 = no result
+
+for i=1:length(MeanVelfs)
+    for j=1:size(MeanVelfs,2)
+        if MeanVelfs(i,j)==0
+            MeanVelfs(i,j)=NaN;
+        end
+    end
+end
+for j=1:size(MaxVelfs,2)
+    MeanMaxfs(j)=mean(MeanVelfs(~isnan(MeanVelfs(:,j)),j));
+    StdDevfs(j)=std(MeanVelfs(~isnan(MeanVelfs(:,j)),j));
+    Numfs(j)=sum(~isnan(MeanVelfs(:,j)));
+end
+for i=1:length(MeanVelnfs)
+    for j=1:size(MeanVelnfs,2)
+        if MeanVelnfs(i,j)==0
+            MeanVelnfs(i,j)=NaN;
+        end
+    end
+end
+for j=1:size(MeanVelnfs,2)
+    MeanMaxnfs(j)=mean(MeanVelnfs(~isnan(MeanVelnfs(:,j)),j));
+    StdDevnfs(j)=std(MeanVelnfs(~isnan(MeanVelnfs(:,j)),j));
+    Numnfs(j)=sum(~isnan(MeanVelnfs(:,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');
+
Index: /issm/trunk/test/IsmipHomResults/plotting_scripts/plot_exp_d.m
===================================================================
--- /issm/trunk/test/IsmipHomResults/plotting_scripts/plot_exp_d.m	(revision 7888)
+++ /issm/trunk/test/IsmipHomResults/plotting_scripts/plot_exp_d.m	(revision 7888)
@@ -0,0 +1,381 @@
+%-----------------------------------------------------------------
+% 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' 'lpe1' 'mbr1' 'mmr1' 'mtk1' ...
+        'oga1' 'oso1' 'rhi1' 'rhi2' 'rhi3' 'rhi4' 'rhi5' 'spr1' ...
+        'ssu1' 'tpa1' 'yko1'};
+% dirs = {'dpo1' 'lpe1' 'rhi4' 'rhi5' 'rhi1' 'oga1'};
+% 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 1 0 ...
+                  1 0 1 0 1 0 0 1 ...
+                  1 0 1];
+% fullStokesMask = [0 0 0 0 1 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 = 'C:\Program Files\MATLAB\R2006a\work\ismip_all\ismip_mat\';
+
+% The experiment we're interested in
+exp = 'd';
+
+% 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);
+
+% 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);
+            xvel = B(:,2);
+            zvel = B(:,3);
+            
+            % Make sure the data is normalized
+            x = x./max(x);
+             
+            x = inpaint_nans(x, 4);
+            xvel = inpaint_nans(xvel, 4);
+            zvel = inpaint_nans(zvel, 4);
+            
+            % Calculates the magnitude of the surface velocity
+            
+            surfVel = sqrt(xvel.^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} = interp1(x, surfVel, newX);
+                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} = interp1(x, surfVel, newX);
+                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 a^{-1})';
+
+% 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 a^{-1})');
+    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, 'expd.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 surface velocities (F. Pattyn)
+% Full Stokes stored in MaxVelfs (else MaxVelnfs)
+% cols = L, rows = participants, 0 = no result
+
+for i=1:length(MeanVelfs)
+    for j=1:size(MeanVelfs,2)
+        if MeanVelfs(i,j)==0
+            MeanVelfs(i,j)=NaN;
+        end
+    end
+end
+for j=1:size(MaxVelfs,2)
+    MeanMaxfs(j)=mean(MeanVelfs(~isnan(MeanVelfs(:,j)),j));
+    StdDevfs(j)=std(MeanVelfs(~isnan(MeanVelfs(:,j)),j));
+    Numfs(j)=sum(~isnan(MeanVelfs(:,j)));
+end
+for i=1:length(MeanVelnfs)
+    for j=1:size(MeanVelnfs,2)
+        if MeanVelnfs(i,j)==0
+            MeanVelnfs(i,j)=NaN;
+        end
+    end
+end
+for j=1:size(MeanVelnfs,2)
+    MeanMaxnfs(j)=mean(MeanVelnfs(~isnan(MeanVelnfs(:,j)),j));
+    StdDevnfs(j)=std(MeanVelnfs(~isnan(MeanVelnfs(:,j)),j));
+    Numnfs(j)=sum(~isnan(MeanVelnfs(:,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');
+
Index: /issm/trunk/test/IsmipHomResults/plotting_scripts/plot_exp_e.asv
===================================================================
--- /issm/trunk/test/IsmipHomResults/plotting_scripts/plot_exp_e.asv	(revision 7888)
+++ /issm/trunk/test/IsmipHomResults/plotting_scripts/plot_exp_e.asv	(revision 7888)
@@ -0,0 +1,380 @@
+%-----------------------------------------------------------------
+% 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' 'spr1' 'dpo1' ...
+        'fpa1' 'fpa2' 'fsa1' 'lpe1' 'mbr1' 'mmr1' 'mtk1' ...
+        'oga1' 'oso1' 'rhi1' 'rhi2' 'rhi3' 'rhi4' 'rhi5' ...
+        'ssu1' 'tpa1'};
+dirs = {'jvj1'};
+% 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 1 ...
+                  0 1 0 1 0 0 1 0 ...
+                  1 0 1 0 1 0 0 ...
+                  1 0];
+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 = 'e';
+
+% The file we're intersted in
+fileNames = {[exp '000'], [exp '001']};
+
+% 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);
+
+% 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);
+            xvel = B(:,2);
+            zvel = B(:,3);
+            
+            % Make sure the data is normalized
+            x = x./max(x);
+             
+            x = inpaint_nans(x, 4);
+            xvel = inpaint_nans(xvel, 4);
+            zvel = inpaint_nans(zvel, 4);
+            
+            % Calculates the magnitude of the surface velocity
+            
+            surfVel = sqrt(xvel.^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} = interp1(x, surfVel, newX);
+                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} = interp1(x, surfVel, newX);
+                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 = {'No-Slip Bed' 'Slip Bed'};
+% The label for the x-axis
+xLabel = 'Normalized x';
+% The label for the y-axis
+yLabel = 'Velocity (m a^{-1})';
+
+% 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, 1, 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 a^{-1})');
+    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);
+
+
+set(gcf, 'paperpositionmode','manual','paperorientation', ...
+        'portrait','papertype','A4','paperunits', 'centimeters', ...
+        'paperposition',[1 4 11 17]);
+saveas(gcf, 'expe.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 surface velocities (F. Pattyn)
+% Full Stokes stored in MaxVelfs (else MaxVelnfs)
+% cols = L, rows = participants, 0 = no result
+
+for i=1:length(MeanVelfs)
+    for j=1:size(MeanVelfs,2)
+        if MeanVelfs(i,j)==0
+            MeanVelfs(i,j)=NaN;
+        end
+    end
+end
+for j=1:size(MaxVelfs,2)
+    MeanMaxfs(j)=mean(MeanVelfs(~isnan(MeanVelfs(:,j)),j));
+    StdDevfs(j)=std(MeanVelfs(~isnan(MeanVelfs(:,j)),j));
+    Numfs(j)=sum(~isnan(MeanVelfs(:,j)));
+end
+for i=1:length(MeanVelnfs)
+    for j=1:size(MeanVelnfs,2)
+        if MeanVelnfs(i,j)==0
+            MeanVelnfs(i,j)=NaN;
+        end
+    end
+end
+for j=1:size(MeanVelnfs,2)
+    MeanMaxnfs(j)=mean(MeanVelnfs(~isnan(MeanVelnfs(:,j)),j));
+    StdDevnfs(j)=std(MeanVelnfs(~isnan(MeanVelnfs(:,j)),j));
+    Numnfs(j)=sum(~isnan(MeanVelnfs(:,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');
Index: /issm/trunk/test/IsmipHomResults/plotting_scripts/plot_exp_e.m
===================================================================
--- /issm/trunk/test/IsmipHomResults/plotting_scripts/plot_exp_e.m	(revision 7888)
+++ /issm/trunk/test/IsmipHomResults/plotting_scripts/plot_exp_e.m	(revision 7888)
@@ -0,0 +1,380 @@
+%-----------------------------------------------------------------
+% 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' 'spr1' 'dpo1' ...
+        'fpa1' 'fpa2' 'fsa1' 'lpe1' 'mbr1' 'mmr1' 'mtk1' ...
+        'oga1' 'oso1' 'rhi1' 'rhi2' 'rhi3' 'rhi4' 'rhi5' ...
+        'ssu1' 'tpa1'};
+% dirs = {'jvj1'};
+% 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 1 ...
+                  0 1 0 0 0 1 0 ...
+                  1 0 1 0 1 0 0 ...
+                  1 0];
+% 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 = 'e';
+
+% The file we're intersted in
+fileNames = {[exp '000'], [exp '001']};
+
+% 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);
+
+% 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);
+            xvel = B(:,2);
+            zvel = B(:,3);
+            
+            % Make sure the data is normalized
+            x = x./max(x);
+             
+            x = inpaint_nans(x, 4);
+            xvel = inpaint_nans(xvel, 4);
+            zvel = inpaint_nans(zvel, 4);
+            
+            % Calculates the magnitude of the surface velocity
+            
+            surfVel = sqrt(xvel.^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} = interp1(x, surfVel, newX);
+                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} = interp1(x, surfVel, newX);
+                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 = {'No-Slip Bed' 'Slip Bed'};
+% The label for the x-axis
+xLabel = 'Normalized x';
+% The label for the y-axis
+yLabel = 'Velocity (m a^{-1})';
+
+% 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, 1, 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 a^{-1})');
+    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);
+
+
+set(gcf, 'paperpositionmode','manual','paperorientation', ...
+        'portrait','papertype','A4','paperunits', 'centimeters', ...
+        'paperposition',[1 4 11 17]);
+saveas(gcf, 'expe.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 surface velocities (F. Pattyn)
+% Full Stokes stored in MaxVelfs (else MaxVelnfs)
+% cols = L, rows = participants, 0 = no result
+
+for i=1:length(MeanVelfs)
+    for j=1:size(MeanVelfs,2)
+        if MeanVelfs(i,j)==0
+            MeanVelfs(i,j)=NaN;
+        end
+    end
+end
+for j=1:size(MaxVelfs,2)
+    MeanMaxfs(j)=mean(MeanVelfs(~isnan(MeanVelfs(:,j)),j));
+    StdDevfs(j)=std(MeanVelfs(~isnan(MeanVelfs(:,j)),j));
+    Numfs(j)=sum(~isnan(MeanVelfs(:,j)));
+end
+for i=1:length(MeanVelnfs)
+    for j=1:size(MeanVelnfs,2)
+        if MeanVelnfs(i,j)==0
+            MeanVelnfs(i,j)=NaN;
+        end
+    end
+end
+for j=1:size(MeanVelnfs,2)
+    MeanMaxnfs(j)=mean(MeanVelnfs(~isnan(MeanVelnfs(:,j)),j));
+    StdDevnfs(j)=std(MeanVelnfs(~isnan(MeanVelnfs(:,j)),j));
+    Numnfs(j)=sum(~isnan(MeanVelnfs(:,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');
Index: /issm/trunk/test/IsmipHomResults/plotting_scripts/plot_exp_e_txz.asv
===================================================================
--- /issm/trunk/test/IsmipHomResults/plotting_scripts/plot_exp_e_txz.asv	(revision 7888)
+++ /issm/trunk/test/IsmipHomResults/plotting_scripts/plot_exp_e_txz.asv	(revision 7888)
@@ -0,0 +1,298 @@
+%-----------------------------------------------------------------
+% 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'};
+% 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];
+% 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 = 'e';
+
+% The file we're intersted in
+fileNames = {[exp '001'], [exp '000']};
+
+% 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);
+
+% 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);
+            xvel = B(:,2);
+            zvel = B(:,3);
+            txz = B(:,4);
+            
+            % Make sure the data is normalized
+            x = x./max(x);
+             
+            x = inpaint_nans(x, 4);
+            xvel = inpaint_nans(xvel, 4);
+            zvel = inpaint_nans(zvel, 4);
+            txz = inpaint_nans(txz, 4);
+            
+            % Calculates the magnitude of the basal shear
+            
+            surfVel = txz;
+
+            % 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} = interp1(x, surfVel, newX);
+                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)))
+                    
+                    % 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} = interp1(x, surfVel, newX);
+                nfsVels{j, nfsIndex} = inpaint_nans(nfsVels{j, nfsIndex}, 4);
+                
+                % 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 = {'Slip Bed' 'No-Slip Bed'};
+% The label for the x-axis
+xLabel = 'Normalized x';
+% The label for the y-axis
+yLabel = 'Basal shear stress (kPa)';
+
+% 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
+    if i==1
+        subplot(2, 1, 2);
+    else
+        subplot(2, 1, 1);
+    end
+    
+    % 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, fsCurDomainVels, 'k');
+    set(findobj(gca,'Type','line','Color',[0 0 0]),'Color',fsFillColor,'LineWidth',2);
+    plot(newX, nfsCurDomainVels, '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);
+
+
+set(gcf, 'paperpositionmode','manual','paperorientation', ...
+        'portrait','papertype','A4','paperunits', 'centimeters', ...
+        'paperposition',[1 4 11 17]);
+saveas(gcf, 'expe.eps', 'psc2');
Index: /issm/trunk/test/IsmipHomResults/plotting_scripts/plot_exp_e_txz.m
===================================================================
--- /issm/trunk/test/IsmipHomResults/plotting_scripts/plot_exp_e_txz.m	(revision 7888)
+++ /issm/trunk/test/IsmipHomResults/plotting_scripts/plot_exp_e_txz.m	(revision 7888)
@@ -0,0 +1,296 @@
+%-----------------------------------------------------------------
+% 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 = {'aas2' 'ahu1' 'ahu2' 'bds1' 'cma1' 'cma2' 'spr1' 'dpo1' ...
+        'fpa1' 'fpa2' 'fsa1' 'lpe1' 'mbr1' 'mmr1' 'mtk1' ...
+        'oga1' 'oso1' 'rhi1' 'rhi2' 'rhi3' 'rhi4' 'rhi5' ...
+        'ssu1' 'tpa1'};
+% dirs = {'dpo1'};
+% 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 1 0 0 1 ...
+                  0 1 0 0 0 1 0 ...
+                  1 0 1 0 1 0 0 ...
+                  1 0];
+% 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 = 'e';
+
+% The file we're intersted in
+fileNames = {[exp '000'], [exp '001']};
+
+% 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);
+
+% 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);
+            xvel = B(:,2);
+            zvel = B(:,3);
+            txz = B(:,4);
+            
+            % Make sure the data is normalized
+            x = x./max(x);
+             
+            x = inpaint_nans(x, 4);
+            xvel = inpaint_nans(xvel, 4);
+            zvel = inpaint_nans(zvel, 4);
+            txz = inpaint_nans(txz, 4);
+            
+            % Calculates the magnitude of the basal shear
+            
+            surfVel = txz;
+
+            % 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} = interp1(x, surfVel, newX);
+                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)))
+                    
+                    % 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} = interp1(x, surfVel, newX);
+                nfsVels{j, nfsIndex} = inpaint_nans(nfsVels{j, nfsIndex}, 4);
+                
+                % 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 = {'No-Slip Bed' 'Slip Bed'};
+% The label for the x-axis
+xLabel = 'Normalized x';
+% The label for the y-axis
+yLabel = 'Basal shear stress (kPa)';
+
+% 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, 1, 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('Basal shear stress (kPa)');
+    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);
+
+
+set(gcf, 'paperpositionmode','manual','paperorientation', ...
+        'portrait','papertype','A4','paperunits', 'centimeters', ...
+        'paperposition',[1 4 11 17]);
+saveas(gcf, 'expetxz.eps', 'psc2');
Index: /issm/trunk/test/IsmipHomResults/plotting_scripts/plot_exp_f.asv
===================================================================
--- /issm/trunk/test/IsmipHomResults/plotting_scripts/plot_exp_f.asv	(revision 7888)
+++ /issm/trunk/test/IsmipHomResults/plotting_scripts/plot_exp_f.asv	(revision 7888)
@@ -0,0 +1,414 @@
+%-----------------------------------------------------------------
+% 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 = {'mtk1'};
+% 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 = 'f';
+
+% The file we're intersted in
+fileNames = {[exp '000'], [exp '001']};
+
+% The desired vector for normalized x.  Used for interpolation so every
+% data set has the same number of points.
+newX = linspace(-50.0,50.0,250);
+newY = 0*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(yvel.^2 + zvel.^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 
+
+% Analytical solution by Gudmundsson
+
+% Create the full file path from the base path and current directory
+fullFilePath = [location 'ghg1' fileNames(1)];
+
+% Open the current file for reading
+load(cell2mat(fullFilePath));
+
+x = B(:,1);
+y = B(:,2);
+xvel = B(:,3);
+yvel = B(:,4);
+zvel = B(:,5);
+            
+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);
+
+anVels{1, 1} = griddata(x, y, surfVel, newX, newY);
+
+% Calculates the magnitude of the surface velocity
+            
+surfVel = sqrt(yvel.^2 + zvel.^2);
+
+%--------------------------------------------------------
+% THIS IS THE PROCESSING AND DISPLAY SECTION WHERE THE 
+% MEAN AND STANDARD DEVIATION ARE CALCULATED AND PLOTTED
+%--------------------------------------------------------
+
+% The titles for the plots
+plotTitles = {'No-Slip Bed' 'Slip Bed'};
+% The label for the x-axis
+xLabel = 'Distance for center (km)';
+% The label for the y-axis
+yLabel = 'Velocity (m a^{-1})';
+
+% 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, 1, 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 a^{-1})');
+    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);
+
+
+set(gcf, 'paperpositionmode','manual','paperorientation', ...
+        'portrait','papertype','A4','paperunits', 'centimeters', ...
+        'paperposition',[1 4 11 17]);
+saveas(gcf, 'expf.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 surface velocities (F. Pattyn)
+% Full Stokes stored in MaxVelfs (else MaxVelnfs)
+% cols = L, rows = participants, 0 = no result
+
+for i=1:length(MeanVelfs)
+    for j=1:size(MeanVelfs,2)
+        if MeanVelfs(i,j)==0
+            MeanVelfs(i,j)=NaN;
+        end
+    end
+end
+for j=1:size(MaxVelfs,2)
+    MeanMaxfs(j)=mean(MeanVelfs(~isnan(MeanVelfs(:,j)),j));
+    StdDevfs(j)=std(MeanVelfs(~isnan(MeanVelfs(:,j)),j));
+    Numfs(j)=sum(~isnan(MeanVelfs(:,j)));
+end
+for i=1:length(MeanVelnfs)
+    for j=1:size(MeanVelnfs,2)
+        if MeanVelnfs(i,j)==0
+            MeanVelnfs(i,j)=NaN;
+        end
+    end
+end
+for j=1:size(MeanVelnfs,2)
+    MeanMaxnfs(j)=mean(MeanVelnfs(~isnan(MeanVelnfs(:,j)),j));
+    StdDevnfs(j)=std(MeanVelnfs(~isnan(MeanVelnfs(:,j)),j));
+    Numnfs(j)=sum(~isnan(MeanVelnfs(:,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');
Index: /issm/trunk/test/IsmipHomResults/plotting_scripts/plot_exp_f.m
===================================================================
--- /issm/trunk/test/IsmipHomResults/plotting_scripts/plot_exp_f.m	(revision 7888)
+++ /issm/trunk/test/IsmipHomResults/plotting_scripts/plot_exp_f.m	(revision 7888)
@@ -0,0 +1,418 @@
+%-----------------------------------------------------------------
+% 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 = {'mtk1'};
+% 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 1 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 = 'f';
+
+% The file we're intersted in
+fileNames = {[exp '000'], [exp '001']};
+
+% The desired vector for normalized x.  Used for interpolation so every
+% data set has the same number of points.
+newX = linspace(-50.0,50.0,250);
+newY = 0*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(yvel.^2 + zvel.^2);
+%             surfVel = yvel;
+
+            % 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 
+
+% Analytical solution by Gudmundsson
+
+% Create the full file path from the base path and current directory
+fullFilePath = [location 'ghg1' fileNames(1)];
+
+% Open the current file for reading
+load(cell2mat(fullFilePath));
+
+x = B(:,1);
+y = B(:,2);
+xvel = B(:,3);
+yvel = B(:,4);
+zvel = B(:,5);
+            
+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(yvel.^2 + zvel.^2);
+% surfVel = yvel;
+anVels = griddata(x, y, surfVel, newX, newY);
+
+%--------------------------------------------------------
+% THIS IS THE PROCESSING AND DISPLAY SECTION WHERE THE 
+% MEAN AND STANDARD DEVIATION ARE CALCULATED AND PLOTTED
+%--------------------------------------------------------
+
+% The titles for the plots
+plotTitles = {'No-Slip Bed' 'Slip Bed'};
+% The label for the x-axis
+xLabel = 'Distance from center (km)';
+% The label for the y-axis
+yLabel = 'Velocity (m a^{-1})';
+
+% 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, 1, 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 a^{-1})');
+    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);
+    
+    if (i==1)
+        plot(newX, anVels, 'k', 'LineWidth', 1.5);
+    end
+end
+
+% Add the legend to the final subplot (in upper left corner)
+legend('FS', 'NFS', 'FS Mean', 'NFS Mean', 2);
+
+
+set(gcf, 'paperpositionmode','manual','paperorientation', ...
+        'portrait','papertype','A4','paperunits', 'centimeters', ...
+        'paperposition',[1 4 11 17]);
+saveas(gcf, 'expf.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 surface velocities (F. Pattyn)
+% Full Stokes stored in MaxVelfs (else MaxVelnfs)
+% cols = L, rows = participants, 0 = no result
+
+for i=1:length(MeanVelfs)
+    for j=1:size(MeanVelfs,2)
+        if MeanVelfs(i,j)==0
+            MeanVelfs(i,j)=NaN;
+        end
+    end
+end
+for j=1:size(MaxVelfs,2)
+    MeanMaxfs(j)=mean(MeanVelfs(~isnan(MeanVelfs(:,j)),j));
+    StdDevfs(j)=std(MeanVelfs(~isnan(MeanVelfs(:,j)),j));
+    Numfs(j)=sum(~isnan(MeanVelfs(:,j)));
+end
+for i=1:length(MeanVelnfs)
+    for j=1:size(MeanVelnfs,2)
+        if MeanVelnfs(i,j)==0
+            MeanVelnfs(i,j)=NaN;
+        end
+    end
+end
+for j=1:size(MeanVelnfs,2)
+    MeanMaxnfs(j)=mean(MeanVelnfs(~isnan(MeanVelnfs(:,j)),j));
+    StdDevnfs(j)=std(MeanVelnfs(~isnan(MeanVelnfs(:,j)),j));
+    Numnfs(j)=sum(~isnan(MeanVelnfs(:,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');
Index: /issm/trunk/test/IsmipHomResults/plotting_scripts/plot_exp_fs.m
===================================================================
--- /issm/trunk/test/IsmipHomResults/plotting_scripts/plot_exp_fs.m	(revision 7888)
+++ /issm/trunk/test/IsmipHomResults/plotting_scripts/plot_exp_fs.m	(revision 7888)
@@ -0,0 +1,417 @@
+%-----------------------------------------------------------------
+% 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 = {'mtk1'};
+% 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 1 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 = 'f';
+
+% The file we're intersted in
+fileNames = {[exp '000'], [exp '001']};
+
+% The desired vector for normalized x.  Used for interpolation so every
+% data set has the same number of points.
+newX = linspace(-50.0,50.0,250);
+newY = 0*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);
+            surfVel = xvel;
+
+            % 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 
+
+% Analytical solution by Gudmundsson
+
+% Create the full file path from the base path and current directory
+fullFilePath = [location 'ghg1' fileNames(1)];
+
+% Open the current file for reading
+load(cell2mat(fullFilePath));
+
+x = B(:,1);
+y = B(:,2);
+xvel = B(:,3);
+yvel = B(:,4);
+zvel = B(:,5);
+            
+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 = xvel;
+anVels = griddata(x, y, surfVel, newX, newY);
+
+%--------------------------------------------------------
+% THIS IS THE PROCESSING AND DISPLAY SECTION WHERE THE 
+% MEAN AND STANDARD DEVIATION ARE CALCULATED AND PLOTTED
+%--------------------------------------------------------
+
+% The titles for the plots
+plotTitles = {'No-Slip Bed' 'Slip Bed'};
+% The label for the x-axis
+xLabel = 'Distance from center (km)';
+% The label for the y-axis
+yLabel = 'Surface (m)';
+
+% 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, 1, 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('Surface (m)');
+    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);
+    
+    if (i==1)
+        plot(newX, anVels, 'k', 'LineWidth', 1.5);
+    end
+end
+
+% Add the legend to the final subplot (in upper left corner)
+legend('FS', 'NFS', 'FS Mean', 'NFS Mean', 2);
+
+
+set(gcf, 'paperpositionmode','manual','paperorientation', ...
+        'portrait','papertype','A4','paperunits', 'centimeters', ...
+        'paperposition',[1 4 11 17]);
+saveas(gcf, 'expfs.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 surface velocities (F. Pattyn)
+% Full Stokes stored in MaxVelfs (else MaxVelnfs)
+% cols = L, rows = participants, 0 = no result
+
+for i=1:length(MeanVelfs)
+    for j=1:size(MeanVelfs,2)
+        if MeanVelfs(i,j)==0
+            MeanVelfs(i,j)=NaN;
+        end
+    end
+end
+for j=1:size(MaxVelfs,2)
+    MeanMaxfs(j)=mean(MeanVelfs(~isnan(MeanVelfs(:,j)),j));
+    StdDevfs(j)=std(MeanVelfs(~isnan(MeanVelfs(:,j)),j));
+    Numfs(j)=sum(~isnan(MeanVelfs(:,j)));
+end
+for i=1:length(MeanVelnfs)
+    for j=1:size(MeanVelnfs,2)
+        if MeanVelnfs(i,j)==0
+            MeanVelnfs(i,j)=NaN;
+        end
+    end
+end
+for j=1:size(MeanVelnfs,2)
+    MeanMaxnfs(j)=mean(MeanVelnfs(~isnan(MeanVelnfs(:,j)),j));
+    StdDevnfs(j)=std(MeanVelnfs(~isnan(MeanVelnfs(:,j)),j));
+    Numnfs(j)=sum(~isnan(MeanVelnfs(:,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');
Index: /issm/trunk/test/IsmipHomResults/plotting_scripts/plot_ismip_exp_a.m
===================================================================
--- /issm/trunk/test/IsmipHomResults/plotting_scripts/plot_ismip_exp_a.m	(revision 7888)
+++ /issm/trunk/test/IsmipHomResults/plotting_scripts/plot_ismip_exp_a.m	(revision 7888)
@@ -0,0 +1,294 @@
+%-----------------------------------------------------------------
+% 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 = 'C:\Documents and Settings\frankp\My Documents\work\ismip_all\ismip_mat\';
+
+% The experiment we're interested in
+exp = 'a';
+
+% 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');
+
+            % 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);
+            
+            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 = 'normal';
+% 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);
Index: /issm/trunk/test/IsmipHomResults/plotting_scripts/plot_ismip_exp_b.m
===================================================================
--- /issm/trunk/test/IsmipHomResults/plotting_scripts/plot_ismip_exp_b.m	(revision 7888)
+++ /issm/trunk/test/IsmipHomResults/plotting_scripts/plot_ismip_exp_b.m	(revision 7888)
@@ -0,0 +1,282 @@
+%-----------------------------------------------------------------
+% 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 = 'b';
+
+% 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);
+
+% 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');
+
+            % Store the values into the current spot in the cell arrays
+            x = data{1,1};
+            xvel = data{1,2};
+            yvel = data{1,3};
+            
+            % Make sure the data is normalized
+            x = x./max(x);
+            
+            % Calculates the magnitude of the surface velocity
+            surfVel = sqrt(xvel.^2 + yvel.^2);
+            surfVel = sqrt(xvel.^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} = interp1(x, surfVel, newX);
+                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} = interp1(x, surfVel, newX);
+                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);
Index: /issm/trunk/test/IsmipHomResults/plotting_scripts/plot_ismip_exp_c.m
===================================================================
--- /issm/trunk/test/IsmipHomResults/plotting_scripts/plot_ismip_exp_c.m	(revision 7888)
+++ /issm/trunk/test/IsmipHomResults/plotting_scripts/plot_ismip_exp_c.m	(revision 7888)
@@ -0,0 +1,295 @@
+%-----------------------------------------------------------------
+% 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);
Index: /issm/trunk/test/IsmipHomResults/plotting_scripts/plot_ismip_exp_d.m
===================================================================
--- /issm/trunk/test/IsmipHomResults/plotting_scripts/plot_ismip_exp_d.m	(revision 7888)
+++ /issm/trunk/test/IsmipHomResults/plotting_scripts/plot_ismip_exp_d.m	(revision 7888)
@@ -0,0 +1,286 @@
+%-----------------------------------------------------------------
+% 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 = [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 ON YOUR MACHINE
+location = '/Applications/MATLAB_R2007b/work/ismip_all/';
+
+% The experiment we're interested in
+exp = 'd';
+
+% 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);
+
+% 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');
+
+            % Store the values into the current spot in the cell arrays
+            x = data{1,1};
+            xvel = data{1,2};
+            yvel = data{1,3};
+            
+            % Make sure the data is normalized
+            x = x./max(x);
+            
+            x = inpaint_nans(x, 4);
+            xvel = inpaint_nans(xvel, 4);
+            yvel = inpaint_nans(yvel, 4);
+            
+            % Calculates the magnitude of the surface velocity
+            surfVel = sqrt(xvel.^2 + yvel.^2);
+            surfVel = sqrt(xvel.^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} = interp1(x, surfVel, newX);
+                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} = interp1(x, surfVel, newX);
+                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);
Index: /issm/trunk/test/IsmipHomResults/plotting_scripts/plot_ismip_exp_e
===================================================================
--- /issm/trunk/test/IsmipHomResults/plotting_scripts/plot_ismip_exp_e	(revision 7888)
+++ /issm/trunk/test/IsmipHomResults/plotting_scripts/plot_ismip_exp_e	(revision 7888)
@@ -0,0 +1,281 @@
+%-----------------------------------------------------------------
+% 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 = [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 = '/home/jfishbaugh/Desktop/frank/ismip_all/';
+
+% The experiment we're interested in
+exp = 'e';
+
+% The file we're intersted in
+fileNames = {[exp '000.txt'] [exp '001.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);
+
+% 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');
+
+            % Store the values into the current spot in the cell arrays
+            x = data{1,1};
+            xvel = data{1,2};
+            yvel = data{1,3};
+            
+            x = inpaint_nans(x, 4);
+            xvel = inpaint_nans(xvel, 4);
+            yvel = inpaint_nans(yvel, 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)
+                
+                % Store the data in the full stokes structure
+                fsVels{j, fsIndex} = interp1(x, surfVel, newX);
+                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} = interp1(x, surfVel, newX);
+                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, 1, 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);
Index: /issm/trunk/test/IsmipHomResults/plotting_scripts/plot_ismip_exp_e.m
===================================================================
--- /issm/trunk/test/IsmipHomResults/plotting_scripts/plot_ismip_exp_e.m	(revision 7888)
+++ /issm/trunk/test/IsmipHomResults/plotting_scripts/plot_ismip_exp_e.m	(revision 7888)
@@ -0,0 +1,281 @@
+%-----------------------------------------------------------------
+% 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' '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 = [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 = '/home/jfishbaugh/Desktop/ismip_all/';
+
+% The experiment we're interested in
+exp = 'e';
+
+% The file we're intersted in
+fileNames = {[exp '000.txt'] [exp '001.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);
+
+% 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');
+
+            % Store the values into the current spot in the cell arrays
+            x = data{1,1};
+            
+            % Make sure the data is normalized
+            x = x./max(x);
+            
+            xvel = data{1,2};
+            yvel = data{1,3};
+            
+            x = inpaint_nans(x, 4);
+            xvel = inpaint_nans(xvel, 4);
+            yvel = inpaint_nans(yvel, 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)
+                
+                % Store the data in the full stokes structure
+                fsVels{j, fsIndex} = interp1(x, surfVel, newX);
+                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
+                
+                                
+            % Else its non-full stokes
+            else
+                
+                % Store the data in the non-full stokes structure
+                nfsVels{j, nfsIndex} = interp1(x, surfVel, newX);
+                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
+                
+                      
+            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;
+                % Increment the non-full stokes index
+                nfsIndex = fsIndex + 1;
+                
+            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 = {'No-Slip Bed' 'Slip Bed'};
+% 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, 1, 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);