source: issm/oecreview/Archive/12321-12677/ISSM-12373-12374.diff@ 12679

Index: /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/model/plot/plot_manager.m
===================================================================
--- /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/model/plot/plot_manager.m       (revision 12373)
+++ /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/model/plot/plot_manager.m       (revision 12374)
@@ -167,6 +167,12 @@
        return;
 end
 
+%Figure out if this is a Profile plot
+if exist(options,'profile')
+       plot_profile(md,data,options,nlines,ncols,i);
+       return;
+end
+
 %process data and model
 [x y z elements is2d isplanet]=processmesh(md,data,options);
 [data2 datatype]=processdata(md,data,options);
Index: /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/model/plot/plot_profile.m
===================================================================
--- /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/model/plot/plot_profile.m       (revision 0)
+++ /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/model/plot/plot_profile.m       (revision 12374)
@@ -0,0 +1,47 @@
+function plot_profile(md,data,options,nlines,ncols,i)
+%PLOT_SECTION - plot a given field on a profile
+%
+%   Usage:
+%      plot_profile(md,data,options,nlines,ncols,i)
+%
+%   See also: PLOTMODEL
+
+%process model
+[x_m y_m z_m elements_m is2d isplanet]=processmesh(md,[],options);
+if ~is2d, error('only 3d model supported'); end
+
+%Get number of curves and generate random colors
+numcurves=size(data,2);
+colorm=getfieldvalue(options,'colormap','lines');
+color=eval([ colorm '(numcurves);']);
+options=removefield(options,'colormap',0); %back to default colormap
+
+%Loop over number of curves
+for i=1:numcurves,
+
+       %Process data
+       [datai datatype]=processdata(md,data(:,i),options);
+
+       %resolution
+       if exist(options,'resolution'),
+               resolution=getfieldvalue(options,'resolution');
+       else %Default resolution
+               resolution=[100];
+               disp(['plot_profile warning: no resolution specified, use default resolution: [horizontal_resolution vertical_resolution]=[' num2str(resolution)  ']']);
+       end
+
+       %Compute profile value
+       [z,data_interp]=ProfileValues(md,data,xprof,yprof,resolution)
+
+       %plot profile
+       subplot(nlines,ncols,i)
+       A=elements(:,1); B=elements(:,2); C=elements(:,3);  D=elements(:,4); 
+       plot(data_interp,z,'color',color(i,:),'LineWidth',getfieldvalue(options,'linewidth',1));
+       hold on;
+end
+
+%apply options
+options=addfielddefault(options,'title','Profile');
+options=addfielddefault(options,'colorbar',0);
+options=addfielddefault(options,'ylabel','z');
+applyoptions(md,[],options);
Index: /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/model/plot/plotdoc.m
===================================================================
--- /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/model/plot/plotdoc.m    (revision 12373)
+++ /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/model/plot/plotdoc.m    (revision 12374)
@@ -93,6 +93,7 @@
 disp('                       horizontal_resolution must be in meter, and vertical_resolution a number of layers');
 disp('       ''showsection'': show section used by ''sectionvalue'' (string ''on'' or a number of labels)');
 disp('       ''sectionvalue'': give the value of data on a profile given by an Argus file (string ''Argusfile_name.exp'')');
+disp('       ''profile'': give the value of data along a vertical profile ([xlocation ylocation])');
 disp('       ''smooth'': smooth element data (string ''yes'' or integer)');
 disp('       ''title'': same as standard matlab option');
 disp('       ''view'': same as standard matlab option (ex: 2, 3 or [90 180]');
Index: /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/model/ProfileValues.m
===================================================================
--- /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/model/ProfileValues.m   (revision 0)
+++ /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/model/ProfileValues.m   (revision 12374)
@@ -0,0 +1,21 @@
+function [Z,data_interp]=ProfileValues(md,data,xprof,yprof,resolution)
+%PROFILEVALUES - compute the value of a field on a vertical profile
+%
+%   This routine gets the value of a given field of the model on points
+%   given by filname (Argus type file)
+%
+%   Usage:
+%      [z,data]=ProfileValues(md,data,filename,resolution)
+%      [z,data]=ProfileValues(md,data,profile_structure,resolution)
+
+%Get bed and surface for each 2d point, offset to make sure that it is inside the glacier system
+offset=10^-3;
+bed=InterpFromMeshToMesh2d(md.mesh.elements2d,md.mesh.x2d,md.mesh.y2d,project2d(md,md.geometry.bed,1),xprof,yprof)+offset;
+surface=InterpFromMeshToMesh2d(md.mesh.elements2d,md.mesh.x2d,md.mesh.y2d,project2d(md,md.geometry.surface,1),xprof,yprof)-offset;
+
+%Some useful parameters
+layers=ceil(mean(md.geometry.thickness)/res_v);
+Z=bed:resolution:surface;
+X=xprof*ones(size(Z));
+Y=yprof*ones(size(Z));
+data_interp=InterpFromMeshToMesh3d(md.mesh.elements,md.mesh.x,md.mesh.y,md.mesh.z,data,X,Y,Z,NaN);