Index: /issm/trunk-jpl/src/m/consistency/QueueRequirements.m =================================================================== --- /issm/trunk-jpl/src/m/consistency/QueueRequirements.m (revision 13010) +++ /issm/trunk-jpl/src/m/consistency/QueueRequirements.m (revision 13010) @@ -0,0 +1,35 @@ +function QueueRequirements(available_queues,queue_requirements_time,queue_requirements_np,queue,np,time) +%QUEUEREQUIREMENTS - queue requirements in time, number of cpus, by name of queue. +% +% Usage: +% QueueRequirements(available_queues,queue_requirements_time,queue_requirements_np,np,time) + +%Ok, go through requirements for current queue: +index=ismemberi(queue,available_queues); +if ~index, + %ok, either we a generic cluster, with 'none' queue, or we could not find the queue reqruirements + if strcmpi(available_queues{1},'none'), + %reset index to 1, so we can fish the requirements + index=1; + else + string=available_queues{1}; + for i=2:length(available_queues), + string=[string ' ' available_queues{i}]; + end + error(['QueueRequirements error message: availables queues are ' string]); + end +end + +%check on time requirements +rtime=queue_requirements_time(index); +if time<=0, + error('QueueRequirements: time should be a positive number'); +end +if time>rtime, + error(['QueueRequirements: time should be < ' num2str(rtime) ' for queue: ' queue]); +end + +%check on np requirements +if np<=0, + error('QueueRequirements: np should be a positive number'); +end Index: /issm/trunk-jpl/src/m/consistency/ismodelselfconsistent.m =================================================================== --- /issm/trunk-jpl/src/m/consistency/ismodelselfconsistent.m (revision 13010) +++ /issm/trunk-jpl/src/m/consistency/ismodelselfconsistent.m (revision 13010) @@ -0,0 +1,100 @@ +function ismodelselfconsistent(md), +%ISMODELSELFCONSISTENT - check that model forms a closed form solvable problem. +% +% Usage: +% ismodelselfconsistent(md), + +%initialize consistency as true +md.private.isconsistent=true; + +%Get solution and associated analyses +solution=md.private.solution; +[analyses,numanalyses]=AnalysisConfiguration(solution); + +%Go through a model field, check that it is a class, and call checkconsistency +fields=properties('model'); +for i=1:length(fields), + field=fields{i}; + + %Some properties do not need to be checked + if ismember(field,{'results' 'debug' 'radaroverlay'}), + continue; + end + + %Check that current field is an object + if ~isobject(md.(field)) + md=checkmessage(md,['field ''' char(field) ''' is not an object']); + end + + %Check consistency of the object + if verLessThan('matlab', '7.6') + md=checkconsistency(md.(field),md,solution,analyses); + else + md=md.(field).checkconsistency(md,solution,analyses); + end +end + +%error message if mode is not consistent +if md.private.isconsistent==false, + error('Model not consistent, see messages above'); +end +end + +function [analyses,numanalyses]=AnalysisConfiguration(solutiontype), % {{{ +%ANALYSISCONFIGURATION - return type of analyses, number of analyses +% +% Usage: +% [analyses, numanalyses]=AnalysisConfiguration(solutiontype); + + + +switch solutiontype, + + case DiagnosticSolutionEnum, + numanalyses=5; + analyses=[DiagnosticHorizAnalysisEnum;DiagnosticVertAnalysisEnum;DiagnosticHutterAnalysisEnum;SurfaceSlopeAnalysisEnum;BedSlopeAnalysisEnum]; + + case SteadystateSolutionEnum, + numanalyses=7; + analyses=[DiagnosticHorizAnalysisEnum;DiagnosticVertAnalysisEnum;DiagnosticHutterAnalysisEnum;SurfaceSlopeAnalysisEnum;BedSlopeAnalysisEnum;ThermalAnalysisEnum;MeltingAnalysisEnum]; + + case ThermalSolutionEnum, + numanalyses=2; + analyses=[ThermalAnalysisEnum;MeltingAnalysisEnum]; + + case EnthalpySolutionEnum, + numanalyses=1; + analyses=[EnthalpyAnalysisEnum]; + + case PrognosticSolutionEnum, + numanalyses=1; + analyses=[PrognosticAnalysisEnum]; + + case BalancethicknessSolutionEnum, + numanalyses=1; + analyses=[BalancethicknessAnalysisEnum]; + + case SurfaceSlopeSolutionEnum, + numanalyses=1; + analyses=[SurfaceSlopeAnalysisEnum]; + + case BedSlopeSolutionEnum, + numanalyses=1; + analyses=[BedSlopeAnalysisEnum]; + + case TransientSolutionEnum, + numanalyses=9; + analyses=[DiagnosticHorizAnalysisEnum;DiagnosticVertAnalysisEnum;DiagnosticHutterAnalysisEnum;SurfaceSlopeAnalysisEnum;BedSlopeAnalysisEnum;ThermalAnalysisEnum;MeltingAnalysisEnum;EnthalpyAnalysisEnum;PrognosticAnalysisEnum]; + + case FlaimSolutionEnum, + numanalyses=1; + analyses=[FlaimAnalysisEnum]; + + case HydrologySolutionEnum, + numanalyses=3; + analyses=[BedSlopeAnalysisEnum;SurfaceSlopeAnalysisEnum;HydrologyAnalysisEnum]; + + otherwise + error('%s%s%s',' solution type: ',EnumToString(solutiontype),' not supported yet!'); + +end % }}} Index: /issm/trunk-jpl/src/m/consistency/ismodelselfconsistent.py =================================================================== --- /issm/trunk-jpl/src/m/consistency/ismodelselfconsistent.py (revision 13010) +++ /issm/trunk-jpl/src/m/consistency/ismodelselfconsistent.py (revision 13010) @@ -0,0 +1,95 @@ +from AnalysisConfiguration import * +from EnumDefinitions import * + +def AnalysisConfiguration(solutiontype): #{{{ + """ + ANALYSISCONFIGURATION - return type of analyses, number of analyses + + Usage: + [analyses, numanalyses]=AnalysisConfiguration(solutiontype); + """ + + if solutiontype == DiagnosticSolutionEnum: + numanalyses=5 + analyses=[DiagnosticHorizAnalysisEnum,DiagnosticVertAnalysisEnum,DiagnosticHutterAnalysisEnum,SurfaceSlopeAnalysisEnum,BedSlopeAnalysisEnum] + + elif solutiontype == SteadystateSolutionEnum: + numanalyses=7 + analyses=[DiagnosticHorizAnalysisEnum,DiagnosticVertAnalysisEnum,DiagnosticHutterAnalysisEnum,SurfaceSlopeAnalysisEnum,BedSlopeAnalysisEnum,ThermalAnalysisEnum,MeltingAnalysisEnum] + + elif solutiontype == ThermalSolutionEnum: + numanalyses=2 + analyses=[ThermalAnalysisEnum,MeltingAnalysisEnum] + + elif solutiontype == EnthalpySolutionEnum: + numanalyses=1 + analyses=[EnthalpyAnalysisEnum] + + elif solutiontype == PrognosticSolutionEnum: + numanalyses=1 + analyses=[PrognosticAnalysisEnum] + + elif solutiontype == BalancethicknessSolutionEnum: + numanalyses=1 + analyses=[BalancethicknessAnalysisEnum] + + elif solutiontype == SurfaceSlopeSolutionEnum: + numanalyses=1 + analyses=[SurfaceSlopeAnalysisEnum] + + elif solutiontype == BedSlopeSolutionEnum: + numanalyses=1 + analyses=[BedSlopeAnalysisEnum] + + elif solutiontype == TransientSolutionEnum: + numanalyses=9 + analyses=[DiagnosticHorizAnalysisEnum,DiagnosticVertAnalysisEnum,DiagnosticHutterAnalysisEnum,SurfaceSlopeAnalysisEnum,BedSlopeAnalysisEnum,ThermalAnalysisEnum,MeltingAnalysisEnum,EnthalpyAnalysisEnum,PrognosticAnalysisEnum] + + elif solutiontype == FlaimSolutionEnum: + numanalyses=1 + analyses=[FlaimAnalysisEnum] + + elif solutiontype == HydrologySolutionEnum: + numanalyses=3 + analyses=[BedSlopeAnalysisEnum,SurfaceSlopeAnalysisEnum,HydrologyAnalysisEnum] + + else: + raise TypeError("solution type: '%s' not supported yet!" % EnumToString(solutiontype)) + + return analyses,numanalyses +#}}} + +def ismodelselfconsistent(md): + """ + ISMODELSELFCONSISTENT - check that model forms a closed form solvable problem. + + Usage: + ismodelselfconsistent(md), + """ + + #initialize consistency as true + md.private.isconsistent=True + + #Get solution and associated analyses + solution=md.private.solution + analyses,numanalyses=AnalysisConfiguration(solution) + + #Go through a model fields, check that it is a class, and call checkconsistency + fields=vars(md) + for field in fields.iterkeys(): + + #Some properties do not need to be checked + if field in ['results','debug','radaroverlay']: + continue + + #Check that current field is an object + if not hasattr(getattr(md,field),'checkconsistency'): + md.checkmessage("field '%s' is not an object." % field) + + #Check consistency of the object + exec("md.%s.checkconsistency(md,solution,analyses)" % field) + + #error message if mode is not consistent + if not md.private.isconsistent: + raise RuntimeError('Model not consistent, see messages above.') + Index: /issm/trunk-jpl/src/m/contrib/ecco/ecco32issm.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/ecco/ecco32issm.m (revision 13010) +++ /issm/trunk-jpl/src/m/contrib/ecco/ecco32issm.m (revision 13010) @@ -0,0 +1,8 @@ +function nodefield=ecco32issm(field,transition,xecco3,yecco3) + + xecco3linear=xecco3(:); yecco3linear=yecco3(:); %linearize + nodefieldlinear=zeros(length(xecco3linear),1); + nodefieldlinear(transition(:,1))=field(transition(:,2)); + nodefield=xecco3; + nodefield(:)=nodefieldlinear; + %nodefield=nodefield'; %not sure we need that Index: /issm/trunk-jpl/src/m/contrib/ecco/issm2ecco3.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/ecco/issm2ecco3.m (revision 13010) +++ /issm/trunk-jpl/src/m/contrib/ecco/issm2ecco3.m (revision 13010) @@ -0,0 +1,8 @@ +function nodefield=issm2ecco3(field,transition,xecco3,yecco3) + + xecco3linear=xecco3(:); yecco3linear=yecco3(:); %linearize + nodefieldlinear=zeros(length(xecco3linear),1); + nodefieldlinear(transition(:,1))=field(transition(:,2)); + nodefield=xecco3; + nodefield(:)=nodefieldlinear; + %nodefield=nodefield'; %not sure we need that Index: /issm/trunk-jpl/src/m/contrib/gslib/gamv.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/gslib/gamv.m (revision 13010) +++ /issm/trunk-jpl/src/m/contrib/gslib/gamv.m (revision 13010) @@ -0,0 +1,70 @@ +function output = gamv(x,y,data,varargin); +%GAMV - use gslib for Kriging +% +% Usage: +% output = gamv(x,y,data,varargin) + +options=pairoptions(varargin{:}); + +nlag = getfieldvalue(options,'nlag', 20); +dlag = getfieldvalue(options,'dlag', 1000); + +%Write data file +fid=fopen('cluster.dat','w'); +fprintf(fid,'%s\n','Data file'); +fprintf(fid,'%i\n',3); +fprintf(fid,'%s\n','Xlocation'); +fprintf(fid,'%s\n','Ylocation'); +fprintf(fid,'%s\n','Data'); +fprintf(fid,'%g %g %g\n',[x y data]'); +fclose(fid); + +%Write parameter file +fid=fopen('gamv.par','w'); +fprintf(fid,'\t\t\t\t%s\n','Parameters for GAMV'); +fprintf(fid,'\t\t\t\t%s\n','*******************'); +fprintf(fid,'\n'); +fprintf(fid,'%s\n','START OF PARAMETERS:'); +fprintf(fid,'%-30s %s\n','./cluster.dat' ,'\file with data'); +fprintf(fid,'%-30s %s\n','1 2 0' ,'\columns for X, Y, Z coordinates'); +fprintf(fid,'%-30s %s\n','1 3 ' ,'\number of variables, column number'); +fprintf(fid,'%-30s %s\n','-1.0e21 1.0e21' ,'\trimming limits'); +fprintf(fid,'%-30s %s\n','gamv.out' ,'\file for variogram output'); +fprintf(fid,'%-30s %s\n',num2str(nlag,'%i') ,'\number of lags'); +fprintf(fid,'%-30s %s\n',num2str(dlag,'%g') ,'\lag separation distance'); +fprintf(fid,'%-30s %s\n',num2str(dlag/2,'%g') ,'\lag tolerance'); +fprintf(fid,'%-30s %s\n','3' ,'\number of directions'); +fprintf(fid,'%-30s %s\n','0.0 90.0 50.0 0.0 90.0 50.0','\azm, atol, bandh, dip, dtol, bandv'); +fprintf(fid,'%-30s %s\n','0.0 22.5 25.0 0.0 22.5 25.0','\azm, atol, bandh, dip, dtol, bandv'); +fprintf(fid,'%-30s %s\n','90. 22.5 25.0 0.0 22.5 25.0','\azm, atol, bandh, dip, dtol, bandv'); +fprintf(fid,'%-30s %s\n','0' ,'\standardize sill? (0=no, 1=yes)'); +fprintf(fid,'%-30s %s\n','2' ,'\number of variograms'); +fprintf(fid,'%-30s %s\n','1 1 1' ,'\tail var., head vars., variogram type'); +fprintf(fid,'%-30s %s\n','1 1 3' ,'\tail var., head vars., variogram type'); +fclose(fid); + +%Call gamv +system([issmdir() '/externalpackages/gslib/install/gamv gamv.par']); +delete('gamv.par'); + +%Read output +output = struct('Semivariogram',[],'Covariance',[]); +counter1 = 1; +counter2 = 1; +fid=fopen('gamv.out','r'); +while (~feof(fid)), + A=fscanf(fid,'%s',1); + if strcmp(A,'Covariance'); + A=fscanf(fid,'%s',4); %Read tail:Data head:Data direction 2 + output(counter1).Covariance=fscanf(fid,'%i %g %g %i %g %g',[6 nlag+2])'; + counter1=counter1+1; + elseif strcmp(A,'Semivariogram'), + A=fscanf(fid,'%s',4); %Read tail:Data head:Data direction 2 + output(counter2).Semivariogram=fscanf(fid,'%i %g %g %i %g %g',[6 nlag+2])'; + counter2=counter2+1; + else + %do nothing + end +end +fclose(fid); +delete('gamv.out') Index: /issm/trunk-jpl/src/m/contrib/gslib/gslib.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/gslib/gslib.m (revision 13010) +++ /issm/trunk-jpl/src/m/contrib/gslib/gslib.m (revision 13010) @@ -0,0 +1,112 @@ +function output = gslib(x,y,data,varargin); +%GSLIB - use gslib for Kriging +% +% Usage: +% output = gslib(x,y,data,varargin) + +%Output Matrix +xmin = xl(1); +ymin = yl(1); +nx = 101; +ny = 101; +deltax = 5000; +deltay = 5000; + +%Variogram +nugget=10; +sill =164; +range =25763; + +%Kriging options +mindata = 1; +maxdata = 50; +maxsearchradius = 50000; + +%Some intermediaries (Convert to gslib's parameters); +c = (sill-nugget); +a = sqrt(3)*range; + +%Write data file +fid=fopen('cluster.dat','w'); +fprintf(fid,'%s\n','Data file'); +fprintf(fid,'%i\n',3); +fprintf(fid,'%s\n','Xlocation'); +fprintf(fid,'%s\n','Ylocation'); +fprintf(fid,'%s\n','Data'); +fprintf(fid,'%g %g %g\n',[x y data]'); +fclose(fid); + +if 0, %GAMV + %Write parameter file + fid=fopen('gamv.par','w'); + fprintf(fid,'\t\t\t\t%s\n','Parameters for GAMV'); + fprintf(fid,'\t\t\t\t%s\n','*******************'); + fprintf(fid,'\n'); + fprintf(fid,'%s\n','START OF PARAMETERS:'); + fprintf(fid,'%-30s %s\n','./cluster.dat' ,'\file with data'); + fprintf(fid,'%-30s %s\n','1 2 0' ,'\columns for X, Y, Z coordinates'); + fprintf(fid,'%-30s %s\n','1 3 ' ,'\number of variables, column number'); + fprintf(fid,'%-30s %s\n','-1.0e21 1.0e21' ,'\trimming limits'); + fprintf(fid,'%-30s %s\n','gamv.out' ,'\file for variogram output'); + fprintf(fid,'%-30s %s\n','20' ,'\number of lags'); + fprintf(fid,'%-30s %s\n','5.0' ,'\lag separation distance'); + fprintf(fid,'%-30s %s\n','3.0' ,'\lag tolerance'); + fprintf(fid,'%-30s %s\n','3' ,'\number of directions'); + fprintf(fid,'%-30s %s\n','0.0 90.0 50.0 0.0 90.0 50.0','\azm, atol, bandh, dip, dtol, bandv'); + fprintf(fid,'%-30s %s\n','0.0 22.5 25.0 0.0 22.5 25.0','\azm, atol, bandh, dip, dtol, bandv'); + fprintf(fid,'%-30s %s\n','90. 22.5 25.0 0.0 22.5 25.0','\azm, atol, bandh, dip, dtol, bandv'); + fprintf(fid,'%-30s %s\n','0' ,'\standardize sill? (0=no, 1=yes)'); + fprintf(fid,'%-30s %s\n','2' ,'\number of variograms'); + fprintf(fid,'%-30s %s\n','1 1 1' ,'\tail var., head vars., variogram type'); + fprintf(fid,'%-30s %s\n','1 1 3' ,'\tail var., head vars., variogram type'); + fclose(fid); + + %Call gamv + system([issmdir() '/externalpackages/gslib/install/gamv gamv.par']); + +else, %Kriging KB2D + %Write parameter file + fid=fopen('kb2d.par','w'); + fprintf(fid,'\t\t\t\t%s\n','Parameters for KB2D'); + fprintf(fid,'\t\t\t\t%s\n','*******************'); + fprintf(fid,'\n'); + fprintf(fid,'%s\n','START OF PARAMETERS:'); + fprintf(fid,'%-30s %s\n','./cluster.dat' ,'\file with data'); + fprintf(fid,'%-30s %s\n','1 2 3' ,'\columns for X, Y and variable'); + fprintf(fid,'%-30s %s\n','-1.0e21 1.0e21' ,'\trimming limits'); + fprintf(fid,'%-30s %s\n','0' ,'\debugging level: 0,1,2,3'); + fprintf(fid,'%-30s %s\n','kb2d.dbg' ,'\file for debuggging output'); + fprintf(fid,'%-30s %s\n','kb2d.out' ,'\file for kriged output'); + fprintf(fid,'%-30s %s\n',num2str([nx xmin deltax],'%i %10g %6g') ,'\nx, xmn, xsiz'); + fprintf(fid,'%-30s %s\n',num2str([ny ymin deltay],'%i %10g %6g') ,'\nx, xmn, xsiz'); + fprintf(fid,'%-30s %s\n','1 1' ,'\x and y block discretization'); + fprintf(fid,'%-30s %s\n',num2str([mindata maxdata],'%6g') ,'\min and max data for kriging'); + fprintf(fid,'%-30s %s\n',num2str(maxsearchradius,'%6g') ,'\max search radius'); + fprintf(fid,'%-30s %s\n','1 2.302' ,'\0=SK, 1=OK, (mean if SK)'); + fprintf(fid,'%-30s %s\n',['1 ' num2str(nugget)] ,'\nst, nugget effect'); + fprintf(fid,'%-30s %s\n',['3 ' num2str([c 0.0 a a],'%10g')],'\it, c, azm, a_max, a_min'); + fclose(fid); + + tic;system([issmdir() '/externalpackages/gslib/install/kb2d kb2d.par']);toc; + delete('kb2d.par'); + + %Read output + fid=fopen('kb2d.out','r'); + while (~feof(fid)), + A=fscanf(fid,'%s',1); + if strcmp(A,'KB2D'); + A=fscanf(fid,'%s',1); %Read output + params=fscanf(fid,'%i %i %i %i %g %g %g %g %g %g %1',[11 1]); + elseif strcmp(A,' Estimate'), + continue; + elseif strcmp(A,'Estimation'), + A=fscanf(fid,'%s',1); %Read Variance + A=fscanf(fid,'%g %g',[params(1) params(2)*params(3)]); + B=A(1,:); B=reshape(B,[params(3),params(2)])'; + E=A(2,:); E=reshape(E,[params(3),params(2)])'; + else + %do nothing + end + end + fclose(fid); +end Index: /issm/trunk-jpl/src/m/contrib/gslib/pkriging.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/gslib/pkriging.m (revision 13010) +++ /issm/trunk-jpl/src/m/contrib/gslib/pkriging.m (revision 13010) @@ -0,0 +1,55 @@ +function [B E]=pkriging(x,y,observations,x_interp,y_interp,varargin); +%PKRIGING - parallel Kriging +% +% Usage: +% [B E]=pkriging(x,y,observations,x_interp,y_interp,varargin); + +options=pairoptions(varargin{:}); +cluster=getfieldvalue(options,'cluster',generic('np',10)); +options=removefield(options,'cluster',0); +name = ['krig' num2str(feature('GetPid'))]; + +% ========================================= MARSHALL.m ================================================= +disp(['marshalling file ' name '.bin']); +fid=fopen([name '.bin'],'wb'); +if fid==-1, + error(['marshall error message: could not open ' name '.bin file for binary writing']); +end + +%First, write MaximumNumberOfEnum to make sure that the Enums are synchronized +WriteData(fid,'enum',MaximumNumberOfEnums(),'data',true,'format','Boolean'); + +%Write all data +WriteData(fid,'enum',0,'data',x,'format','DoubleMat'); +WriteData(fid,'enum',1,'data',y,'format','DoubleMat'); +WriteData(fid,'enum',2,'data',observations,'format','DoubleMat'); +WriteData(fid,'enum',3,'data',x_interp,'format','DoubleMat'); +WriteData(fid,'enum',4,'data',y_interp,'format','DoubleMat'); +options.marshall(fid,5); +st=fclose(fid); +if st==-1, + error(['marshall error message: could not close file ' name '.bin']); +end +% ========================================= MARSHALL.m ================================================= + +%Launch job on remote cluster +BuildKrigingQueueScript(cluster,name,'',1,0,0); %gather, valgrind, gprof +tic +LaunchQueueJob(cluster,name,name,{[name '.bin'] [name '.queue']}); +toc +choice=input('Is the job successfully completed? (y/n)','s'); +Download(cluster,name,{[name '.outbin']}); +structure=parseresultsfromdisk([name '.outbin'],0); +delete([name '.outlog']); +delete([name '.errlog']); +delete([name '.outbin']); +delete([name '.bin']); +if ~ispc, + delete([name '.tar.gz']); +end + +%Process results +B=structure.AutodiffForward; +B=reshape(B,size(x_interp,2),size(x_interp,1))'; +E=structure.AutodiffIsautodiff; +E=reshape(E,size(x_interp,2),size(x_interp,1))'; Index: /issm/trunk-jpl/src/m/contrib/gslib/varmap.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/gslib/varmap.m (revision 13010) +++ /issm/trunk-jpl/src/m/contrib/gslib/varmap.m (revision 13010) @@ -0,0 +1,55 @@ +function output = varmap(x,y,data,varargin); +%VARMAP - use gslib for Kriging +% +% Usage: +% output = varmap(x,y,data,varargin) + +options=pairoptions(varargin{:}); + +nxlag = getfieldvalue(options,'nxlag', 20); +nylag = getfieldvalue(options,'nylag', 20); +dxlag = getfieldvalue(options,'dxlag', 1000); +dylag = getfieldvalue(options,'dylag', 1000); + +%Write data file +fid=fopen('cluster.dat','w'); +fprintf(fid,'%s\n','Data file'); +fprintf(fid,'%i\n',3); +fprintf(fid,'%s\n','Xlocation'); +fprintf(fid,'%s\n','Ylocation'); +fprintf(fid,'%s\n','Data'); +fprintf(fid,'%g %g %g\n',[x y data]'); +fclose(fid); + +%Write parameter file +fid=fopen('varmap.par','w'); +fprintf(fid,'\t\t\t\t%s\n','Parameters for GAMV'); +fprintf(fid,'\t\t\t\t%s\n','*******************'); +fprintf(fid,'\n'); +fprintf(fid,'%s\n','START OF PARAMETERS:'); +fprintf(fid,'%-30s %s\n','./cluster.dat' ,'\file with data'); +fprintf(fid,'%-30s %s\n','1 3 ' ,'\number of variables, column number'); +fprintf(fid,'%-30s %s\n','-1.0e21 1.0e21' ,'\trimming limits'); +fprintf(fid,'%-30s %s\n','0 ' ,'\1=regular grid, 0=scattered values'); +fprintf(fid,'%-30s %s\n','50 50 1' ,'\if =1: nx, ny, nz'); +fprintf(fid,'%-30s %s\n','1.0 1.0 1.0' ,'\ xsiz, ysiz, zsiz if igrid=1'); +fprintf(fid,'%-30s %s\n','1 2 0' ,'\if =0: columns for x, y and z coordinates'); +fprintf(fid,'%-30s %s\n','varmap.out' ,'\file for variogram output'); +fprintf(fid,'%-30s %s\n',num2str([nxlag nylag 0],'%i '),'\nxlag, nylag, nzlag'); +fprintf(fid,'%-30s %s\n',num2str([dxlag dylag 1],'%g %g %i'),'\dxlag, dylag, dzlag'); +fprintf(fid,'%-30s %s\n','5' ,'\minimum number of pairs'); +fprintf(fid,'%-30s %s\n','0' ,'\standardize sill? (0=no, 1=yes)'); +fprintf(fid,'%-30s %s\n','1' ,'\number of variograms'); +fprintf(fid,'%-30s %s\n','1 1 1' ,'\tail, head, variogram type'); +fclose(fid); + +%Call varmap +system([issmdir() '/externalpackages/gslib/install/varmap varmap.par']); +delete('varmap.par'); + +%Read output +fid=fopen('varmap.out','r'); +A = textscan(fid,'%f %f %f %f %f %f','headerlines',8); +fclose(fid); +delete('varmap.out') +output = reshape(A{1},[2*nxlag+1 2*nylag+1]); Index: /issm/trunk-jpl/src/m/contrib/uci/expremovestraightsegments.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/uci/expremovestraightsegments.m (revision 13010) +++ /issm/trunk-jpl/src/m/contrib/uci/expremovestraightsegments.m (revision 13010) @@ -0,0 +1,31 @@ +function expremovestraightsegments(newfilename,filename,cutoff) +%EXPREMOVESTRAIGHTSEGMENTS: remove straight segments connecting contours. +% +% Usage: expremovestraightsegments('argus.exp',100); +% +% + +a=expread(filename,1); +newcontours=a(1); + +for i=1:length(a), + contour=a(i); + + s=sqrt(contour.x.^2+contour.y.^2); + d=diff(s); + + pos=find(abs(d)>cutoff); + pos=[0;pos;length(contour.x)]; + + for j=1:length(pos)-1, + + newcontour=contour; + newcontour.x=contour.x(pos(j)+1:pos(j+1)); + newcontour.y=contour.y(pos(j)+1:pos(j+1)); + newcontour.nods=length(newcontour.x); + newcontours(end+1)=newcontour; + end +end +newcontours=newcontours(2:end); + +expwrite(newcontours,newfilename); Index: /issm/trunk-jpl/src/m/contrib/uci/expsplit.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/uci/expsplit.m (revision 13010) +++ /issm/trunk-jpl/src/m/contrib/uci/expsplit.m (revision 13010) @@ -0,0 +1,31 @@ +function expsplit(domainoutline) +%EXPSPLIT - split exp file into sub-contours +% +% This routine reads in a domain outline file (Argus format) and plots all the contours +% This will create as many files there are contours in the domain, each file will be postfix with _i +% where i is the contour name. +% +% Usage: +% expsplit(domainoutline) +% +% Example: +% expsplit('Domain.exp'); +% +% See also EXPMASTER, EXPDOC + +%check nargin +if ~nargin | nargin>1 + help expsplit + error('expsplit error message: bad usage'); +end + +[path,root,ext]=fileparts(domainoutline); + +%Read file: +domains=expread(domainoutline,1); + +%split and write contours: +for i=1:length(domains), + subdomain=domains(i); + expwrite(subdomain,[root '_' num2str(i) ext]); +end Index: /issm/trunk-jpl/src/m/contrib/uci/pargenerate.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/uci/pargenerate.m (revision 13010) +++ /issm/trunk-jpl/src/m/contrib/uci/pargenerate.m (revision 13010) @@ -0,0 +1,45 @@ +function pargenerate(filename,xm,ym,data,varargin) +%PARGENERATE - generate parameter file for images ala Rignot +% +% Usage: +% pargenerate(filename,data,xm,ym,options); +% +% Supported options: +% - title: dataset title +% - latitude: standard latitude (degree) +% - meridian: meridian (degree) +% - format: binary format + +%process options +options=pairoptions(varargin{:}); +Title=getfieldvalue(options,'title','N/A'); +latitude=getfieldvalue(options,'latitude','N/A'); +meridian=getfieldvalue(options,'meridian','N/A'); +format=getfieldvalue(options,'format','single'); + +%Get data info +[nlines ncols]=size(data); +xmin=min(xm); +ymax=max(ym); +postx=abs(xm(2)-xm(1)); +posty=abs(ym(2)-ym(1)); + +%Open header file and get machine type +fid=fopen(filename,'wt'); +[filename, permission, machineformat, encoding] = fopen(fid); + +%write header file +fprintf(fid,'%s\n','ISSM gridded dataset parameter file'); +fprintf(fid,'%s%s\n','title: ',Title); +fprintf(fid,'%s \n','DEM_projection: PS'); +fprintf(fid,'%s%s\n','data_format: ',format); +fprintf(fid,'%s%s\n','endian: ',machineformat); +fprintf(fid,'%s%d\n','width: ',ncols); +fprintf(fid,'%s%d\n','nlines: ',nlines); +fprintf(fid,'%s%g%s\n','PS_secant_lat: ',latitude,' decimal degrees'); +fprintf(fid,'%s%g%s\n','PS_meridian_long: ',meridian,' decimal degrees'); +fprintf(fid,'%s%-15.10g\n','PS_corner_north: ',ymax); +fprintf(fid,'%s%-15.10g\n','PS_corner_east: ',xmin); +fprintf(fid,'%s%g%s\n','PS_post_north: ',postx,' m'); +fprintf(fid,'%s%g%s\n','PS_post_east: ',posty,' m'); +fclose(fid); Index: /issm/trunk-jpl/src/m/exp/expbox.m =================================================================== --- /issm/trunk-jpl/src/m/exp/expbox.m (revision 13009) +++ /issm/trunk-jpl/src/m/exp/expbox.m (revision 13010) @@ -1,4 +1,4 @@ function expbox(filename) -%EXPBOX - Create a ARGUS file using to clicks +%EXPBOX - Create an ARGUS file using two clicks % % Two clicks on a plot are used to generate a rectangular box Index: sm/trunk-jpl/src/m/exp/expboxgen.m =================================================================== --- /issm/trunk-jpl/src/m/exp/expboxgen.m (revision 13009) +++ (revision ) @@ -1,99 +1,0 @@ -function expboxgen(x0,y0,nx,ny,parameter_filename,box_filename) -%EXPBOXGEN - creates a domain outline box for a tiff image -% -% This function creates a domain outline box for a .tif image for the -% mosaic tiff. -% x0,y0 and x1,y1 are the cropping coordinates (upper left and lower right -% corners in the larger tiff image). -% paramter_filename is the parameter file name for the mosaic tiff image. -% box_filename is self explanatory. -% -% Usage: -% expboxgen(x0,y0,nx,ny,parameter_filename,box_filename) -% -% See also EXPMASTER, EXPDOC - -%initialize -nrows=-1; -ncols=-1; -X0=-1; -Y0=-1; -posting=-1; - -%first recover necessary information from the paramter_filename. -fid=fopen(parameter_filename,'r'); -count=1; -while(1), -tline=fgetl(fid); -if ~isstr(tline), break; end; - -ind=findstr('ant125m',tline); -if ~isempty(ind), -a=strsplit(tline,' '); -a=a(2); a=char(a); a=a(1:length(a)-1); -posting=str2num(a); -end - - -ind=findstr('no. of rows :',tline); -if ~isempty(ind), -a=strsplit(tline,' '); -nrows=str2num(char(a(length(a)))); -end - -ind=findstr('no. of columns:',tline); -if ~isempty(ind), -a=strsplit(tline,' '); -ncols=str2num(char(a(length(a)))); -end - -ind=findstr('Upper Left X:',tline); -if ~isempty(ind), -a=strsplit(tline,' '); -X0=str2num(char(a(length(a)))); -end - -ind=findstr('Upper Left Y:',tline); -if ~isempty(ind), -a=strsplit(tline,' '); -Y0=str2num(char(a(length(a)))); -end - - -end %while(1), - - -fclose(fid); - -if (X0==-1 | Y0==-1 | nrows==-1 | ncols==-1 | posting==-1), -disp(' '); -disp(['Could not recover all parameters from ' parameter_filename]); -disp('Here are the paramters recovered thus far: '); -disp(['no. of rows: ' num2str(nrows)]); -disp(['no. of columns: ' num2str(ncols)]); -disp(['Upper Left X: ' num2str(X0)]); -disp(['Upper Left Y: ' num2str(Y0)]); -disp(['Posting: ' num2str(posting)]); -end - -disp(' '); -disp(['Recovered the following parameters from ' parameter_filename]); -disp(['no. of rows: ' num2str(nrows)]); -disp(['no. of columns: ' num2str(ncols)]); -disp(['Upper Left X: ' num2str(X0)]); -disp(['Upper Left Y: ' num2str(Y0)]); -disp(['Posting: ' num2str(posting)]); - -%Create X,Y, arrays of coordinates: -X(1)=X0+x0*posting -Y(1)=Y0-y0*posting -X(2)=X(1); -Y(2)=Y(1)-ny*posting; -X(3)=X(1)+nx*posting; -Y(3)=Y(2); -Y(4)=Y(1); -X(4)=X(3); -plot(X,Y,'r*'); - -%Create box exp file using X and Y. Loop it. -expgen(box_filename,X,Y,1); Index: /issm/trunk-jpl/src/m/exp/expcircle.m =================================================================== --- /issm/trunk-jpl/src/m/exp/expcircle.m (revision 13010) +++ /issm/trunk-jpl/src/m/exp/expcircle.m (revision 13010) @@ -0,0 +1,31 @@ +function expcircle(filename,x0,y0,radius,numberofnodes) +%EXPCIRCLE - create a circular contour corresponding to given parameters +% +% Creates a closed argus contour centered on x,y of radius size. +% The contour is made of numberofnodes +% +% Usage: +% expcircle(filename,x0,y0,radius,numberofnodes) +% +% See also EXPMASTER, EXPDOC + +%Calculate the cartesians coordinates of the points +x_list=ones(numberofnodes+1,1); +y_list=ones(numberofnodes+1,1); + +theta=(0:2*pi/numberofnodes:2*pi*(1-1/numberofnodes))'; +theta=[theta;0]; + +x_list=radius*x_list.*cos(theta); +y_list=radius*y_list.*sin(theta); + +%offset x_list and y_list by x0 and y0: +x_list=x_list+x0; +y_list=y_list+y0; + +contour.x=x_list; +contour.y=y_list; +contour.density=1; +contour.name='circle'; + +expwrite(contour,filename); Index: sm/trunk-jpl/src/m/exp/expcreatecircle.m =================================================================== --- /issm/trunk-jpl/src/m/exp/expcreatecircle.m (revision 13009) +++ (revision ) @@ -1,31 +1,0 @@ -function expcreatecircle(filename,x0,y0,radius,numberofnodes) -%EXPCREATECIRCLE - create a circular contour corresponding to given parameters -% -% Creates a closed argus contour centered on x,y of radius size. -% The contour is made of numberofnodes -% -% Usage: -% expcreatecircle(filename,x0,y0,radius,numberofnodes) -% -% See also EXPMASTER, EXPDOC - -%Calculate the cartesians coordinates of the points -x_list=ones(numberofnodes+1,1); -y_list=ones(numberofnodes+1,1); - -theta=(0:2*pi/numberofnodes:2*pi*(1-1/numberofnodes))'; -theta=[theta;0]; - -x_list=radius*x_list.*cos(theta); -y_list=radius*y_list.*sin(theta); - -%offset x_list and y_list by x0 and y0: -x_list=x_list+x0; -y_list=y_list+y0; - -contour.x=x_list; -contour.y=y_list; -contour.density=1; -contour.name='circle'; - -expwrite(contour,filename); Index: /issm/trunk-jpl/src/m/exp/expcreatecontour.m =================================================================== --- /issm/trunk-jpl/src/m/exp/expcreatecontour.m (revision 13009) +++ /issm/trunk-jpl/src/m/exp/expcreatecontour.m (revision 13010) @@ -16,5 +16,5 @@ %Get contour disp('Click on contour points you desire. Type RETURN to end input of points'); -[x,y]=ginputquick; +[x,y]=ginput(); %close contour Index: /issm/trunk-jpl/src/m/exp/expcreateprofile.m =================================================================== --- /issm/trunk-jpl/src/m/exp/expcreateprofile.m (revision 13009) +++ /issm/trunk-jpl/src/m/exp/expcreateprofile.m (revision 13010) @@ -16,5 +16,5 @@ %Get profile disp('Click on profile points you desire. Type RETURN to end input of points'); -[x,y]=ginputquick; +[x,y]=ginput(); %plot contour Index: sm/trunk-jpl/src/m/exp/expdoc.m =================================================================== --- /issm/trunk-jpl/src/m/exp/expdoc.m (revision 13009) +++ (revision ) @@ -1,36 +1,0 @@ -function expdoc() -%EXPDOC - create a doc for EXP routines -% -% Usage: -% expdoc() - -disp(' '); -disp(' Utilities dedicated to ARGUS Files (''.exp'' extension)'); -disp(' expboxgen: creates a domain outline box for a .tif image for the mosaic tiff'); -disp(' usage: expboxgen(x0,y0,nx,ny,parameter_filename,box_filename)'); -disp(' x0,y0 and nx,ny are the cropping coordinates (upper left and lower right corners in the larger tiff image)'); -disp(' paramter_filename is the parameter file name for the mosaic tiff image. box_filename is self explanatory.'); -disp(' expconcatenate: display all the profiles present in an Argus file and merge them accordingly to the tips selected by the user') -disp(' usage: expconcatenate(newfile,oldfile)') -disp(' expcreatecircle: creates a closed Argus contour centered on (x,y), of radius radius_length. The contour is made of N points'); -disp(' usage: expcreatecircle(filename,x0,y0,radius_length,N)'); -disp(' expcreatecontour: creates a closed Argus contour delimited by the ''clicks'' of the user'); -disp(' usage: expcreatecontour(filename)'); -disp(' expcreateprofile: creates an Argus profile delimited by the ''clicks'' of the user'); -disp(' usage: expcreateprofile(filename)'); -disp(' expcut: display the contents of an Argus file and allow the user cut the profiles'); -disp(' usage: expcut(newfilename,oldfilename)'); -disp(' expdisp: display the contents of an Argus file'); -disp(' usage: expdisp(filename,[figure number],[line style])'); -disp(' expgen: creates an Argus file from a contour (x,y) and a flag indicating if the contour must be closed or open'); -disp(' usage: expgen(filename,contours,close_flag)'); -disp(' explink: takes a domain outline made of various segments, and links them together in one domain outline. Use expview to see end result'); -disp(' usage: explink(domainoutline,minthreshold,step)'); -disp(' exptool: allows the user to create, close, merge, remove,... Argus files and save the result in newfile'); -disp(' usage: exptool(newfile,[oldfile1],[oldfile2],[oldfile3],[...]'); -disp(' expread: reads an Argus file and build a structure that holds all the information of the file'); -disp(' usage: expread(file,close_flag)'); -disp(' expselect: display all the profiles of oldfile, the user clicks on the contour he/she wants to remove. Results saved in newfile'); -disp(' usage: expselect(newfile,oldfile)'); -disp(' expwrite: writes an Argus file from a structure given in input???'); -disp(' usage: expwrite(structure,filename)'); Index: /issm/trunk-jpl/src/m/exp/expexcludeoutliers.m =================================================================== --- /issm/trunk-jpl/src/m/exp/expexcludeoutliers.m (revision 13009) +++ /issm/trunk-jpl/src/m/exp/expexcludeoutliers.m (revision 13010) @@ -1,4 +1,5 @@ function excludeoutliers(newcontourname,contourname,domainname) -%EXCLUDEOUTLIERS exclude points of contour that are not within the domain contour. return new contours in a different file. +%EXCLUDEOUTLIERS - exclude points of contour that are not within the domain +%contour. return new contours in a different file. % % Usage: excludeoutliers('NewContour.exp','Contour.exp','DomainOutline.exp'); Index: sm/trunk-jpl/src/m/exp/expgen.m =================================================================== --- /issm/trunk-jpl/src/m/exp/expgen.m (revision 13009) +++ (revision ) @@ -1,47 +1,0 @@ -function expgen(file_name,contours,close_flag); -%EXPGEN - create an Argus file from x and y arrays -% -% Create .exp domain outline type out of x,y coordinates. -% The contour defined by arrays x and y should not be closed. -% Generated domain outline will not be closed, except if close_flag is set to 1. -% -% Usage: -% expgen(file_name,contours,close_flag) -% -% See also EXPMASTER, EXPDOC - -%Check on inputs -if((close_flag~=0) & (close_flag~=1)), -error('close flag must be 0 of 1'); -end -fid=fopen(file_name,'wt'); - -for i=1:length(contours), - if(length(contours(i).x)~=length(contours(i).y)), - error('contours x and y coordinates must be of identical size'); - end - - %get density for this profile. - if isfield(contours,'density'), - density=contours(i).density; - end - - fprintf(fid,'%s\n','## Name:'); - fprintf(fid,'%s\n','## Icon:0'); - fprintf(fid,'%s\n','# Points Count Value'); - if(close_flag==0), - fprintf(fid,'%i %i\n',length(contours(i).x),density); - else - fprintf(fid,'%i %i\n',length(contours(i).x)+1,density); - end - fprintf(fid,'%s\n','# X pos Y pos'); - for j=1:length(contours(i).x), - fprintf(fid,'%f %f\n',contours(i).x(j),contours(i).y(j)); - end - - if(close_flag==1), - fprintf(fid,'%f %f\n',contours(i).x(1),contours(i).y(1)); - end - fprintf(fid,'%s\n',''); -end -fclose(fid); Index: /issm/trunk-jpl/src/m/exp/explink.m =================================================================== --- /issm/trunk-jpl/src/m/exp/explink.m (revision 13009) +++ /issm/trunk-jpl/src/m/exp/explink.m (revision 13010) @@ -27,3 +27,2 @@ end end - Index: sm/trunk-jpl/src/m/exp/expmaster.m =================================================================== --- /issm/trunk-jpl/src/m/exp/expmaster.m (revision 13009) +++ (revision ) @@ -1,2 +1,0 @@ -function expmaster(newfile,varargin) - disp('expmaster has been renamed exptool due to the unpopularity of its name') Index: sm/trunk-jpl/src/m/exp/exporientation.m =================================================================== --- /issm/trunk-jpl/src/m/exp/exporientation.m (revision 13009) +++ (revision ) @@ -1,12 +1,0 @@ -function exporientation(filename) - -a=expread(filename); - -dx=diff(a.x); -dx=[dx;dx(end)]; - -dy=diff(a.y); -dy=[dy;dy(end)]; - -quiver(a.x,a.y,dx,dy); - Index: /issm/trunk-jpl/src/m/exp/expread.m =================================================================== --- /issm/trunk-jpl/src/m/exp/expread.m (revision 13010) +++ /issm/trunk-jpl/src/m/exp/expread.m (revision 13010) @@ -0,0 +1,81 @@ +function Struct=expread(filename); +%EXPREAD - read a file exp and build a Structure +% +% This routine reads a file .exp and build a Structure containing the +% fields x and y corresponding to the coordinates, one for the filename of +% the exp file, for the density, for the nodes, and a field closed to +% indicate if the domain is closed. +% The first argument is the .exp file to be read and the second one (optional) +% indicate if the last point shall be read (1 to read it, 0 not to). +% +% Usage: +% Struct=expread(filename) +% +% Example: +% Struct=expread('domainoutline.exp') +% Struct=expread('domainoutline.exp') +% +% See also EXPDOC, EXPWRITEASVERTICES + +%some checks +if ~exist(filename), + error(['expread error message: file ' filename ' not found!']); +end + +%initialize number of profile +count=0; + +%open file +fid=fopen(filename,'r'); + +%loop over the number of profiles +while (~feof(fid)), + + %update number of profiles + count=count+1; + + %Get file name + A=fscanf(fid,'%s %s',2); + if ~strncmp(A,'##Name:',7), break; end + if length(A)>7, + Struct(count).name=A(8:end); + else + Struct(count).name=''; + end + + %Get Icon + A=fscanf(fid,'%s %s',2); + if ~strncmp(A,'##Icon:',6), break; end + + %Get Info + A=fscanf(fid,'%s %s %s %s',4); + if ~strncmp(A,'#Points',7), break; end + + %Get number of nods and density + A=fscanf(fid,'%f %f',[1 2]); + Struct(count).nods=A(1); + Struct(count).density=A(2); + + %Get Info + A=fscanf(fid,'%s %s %s %s',5); + if ~strncmp(A,'#XposYpos',9), break; end + + %Get Coordinates + A=fscanf(fid,'%f %f',[2 Struct(count).nods]); + Struct(count).x=A(1,:)'; + Struct(count).y=A(2,:)'; + + if(Struct(count).nods~=length(Struct(count).x))error(['Profile ' num2str(count) ' reports incorrect length']); end; + + %Check if closed + if (Struct(count).nods > 1) && ... + (Struct(count).x(end) == Struct(count).x(1)) && ... + (Struct(count).y(end) == Struct(count).y(1)) + Struct(count).closed=true; + else + Struct(count).closed=false; + end +end + +%close file +fclose(fid); Index: sm/trunk-jpl/src/m/exp/expremovestraightsegments.m =================================================================== --- /issm/trunk-jpl/src/m/exp/expremovestraightsegments.m (revision 13009) +++ (revision ) @@ -1,31 +1,0 @@ -function expremovestraightsegments(newfilename,filename,cutoff) -%EXPREMOVESTRAIGHTSEGMENTS: remove straight segments connecting contours. -% -% Usage: expremovestraightsegments('argus.exp',100); -% -% - -a=expread(filename,1); -newcontours=a(1); - -for i=1:length(a), - contour=a(i); - - s=sqrt(contour.x.^2+contour.y.^2); - d=diff(s); - - pos=find(abs(d)>cutoff); - pos=[0;pos;length(contour.x)]; - - for j=1:length(pos)-1, - - newcontour=contour; - newcontour.x=contour.x(pos(j)+1:pos(j+1)); - newcontour.y=contour.y(pos(j)+1:pos(j+1)); - newcontour.nods=length(newcontour.x); - newcontours(end+1)=newcontour; - end -end -newcontours=newcontours(2:end); - -expwrite(newcontours,newfilename); Index: sm/trunk-jpl/src/m/exp/expsplit.m =================================================================== --- /issm/trunk-jpl/src/m/exp/expsplit.m (revision 13009) +++ (revision ) @@ -1,31 +1,0 @@ -function expsplit(domainoutline) -%EXPSPLIT - split exp file into sub-contours -% -% This routine reads in a domain outline file (Argus format) and plots all the contours -% This will create as many files there are contours in the domain, each file will be postfix with _i -% where i is the contour name. -% -% Usage: -% expsplit(domainoutline) -% -% Example: -% expsplit('Domain.exp'); -% -% See also EXPMASTER, EXPDOC - -%check nargin -if ~nargin | nargin>1 - help expsplit - error('expsplit error message: bad usage'); -end - -[path,root,ext]=fileparts(domainoutline); - -%Read file: -domains=expread(domainoutline,1); - -%split and write contours: -for i=1:length(domains), - subdomain=domains(i); - expwrite(subdomain,[root '_' num2str(i) ext]); -end Index: /issm/trunk-jpl/src/m/exp/expsquare.m =================================================================== --- /issm/trunk-jpl/src/m/exp/expsquare.m (revision 13009) +++ /issm/trunk-jpl/src/m/exp/expsquare.m (revision 13010) @@ -1,6 +1,6 @@ -function expbox(filename) +function expsquare(filename) %EXPBOX - Create a ARGUS file using to clicks % -% Two clicks on a plot are used to generate a rectangular box +% Two clicks on a plot are used to generate a square box % This box is written in EXP format on filename % @@ -18,5 +18,5 @@ %Get points -disp('Click twice to define a rectangular domain. First click for upper left corner, second for lower right corner'); +disp('Click twice to define a square domain. First click for upper left corner, second for lower right corner'); [x,y]=ginput(2); Index: /issm/trunk-jpl/src/m/exp/expwrite.m =================================================================== --- /issm/trunk-jpl/src/m/exp/expwrite.m (revision 13010) +++ /issm/trunk-jpl/src/m/exp/expwrite.m (revision 13010) @@ -0,0 +1,41 @@ +function expwrite(a,filename); +%EXPWRITE - write an Argus file from a structure given in input +% +% This routine write an Argus file form a structure containing the fields: +% x and y of the coordinates of the points. +% The first argument is the structure containing the points coordinates +% and the second one the file to be write. +% +% Usage: +% expwrite(a,filename) +% +% Example: +% expwrite(coordstruct,'domainoutline.exp') +% +% See also EXPDOC, EXPREAD, EXPWRITEASVERTICES + +fid=fopen(filename,'w'); +for n=1:length(a), + if(length(contours(n).x)~=length(contours(n).y)), + error('contours x and y coordinates must be of identical size'); + end + + if isfield(a,'name'), + if ~isempty(a(n).name), + fprintf(fid,'%s%s\n','## Name:',a(n).name); + else + fprintf(fid,'%s\n','## Name:'); + end + else + fprintf(fid,'%s\n','## Name:'); + end + + fprintf(fid,'%s\n','## Icon:0'); + fprintf(fid,'%s\n','# Points Count Value'); + fprintf(fid,'%i %f\n',[length(a(n).x) a(n).density]); + fprintf(fid,'%s\n','# X pos Y pos'); + fprintf(fid,'%10.10f %10.10f\n',[a(n).x a(n).y]'); + fprintf(fid,'\n',''); + +end +fclose(fid); Index: sm/trunk-jpl/src/m/exp/ginputquick.m =================================================================== --- /issm/trunk-jpl/src/m/exp/ginputquick.m (revision 13009) +++ (revision ) @@ -1,221 +1,0 @@ -function [out1,out2,out3] = ginput(arg1) -%GINPUT - Graphical input from mouse. -% -% [X,Y] = GINPUT(N) gets N points from the current axes and returns -% the X- and Y-coordinates in length N vectors X and Y. The cursor -% can be positioned using a mouse (or by using the Arrow Keys on some -% systems). Data points are entered by pressing a mouse button -% or any key on the keyboard except carriage return, which terminates -% the input before N points are entered. -% -% [X,Y] = GINPUT gathers an unlimited number of points until the -% return key is pressed. -% -% [X,Y,BUTTON] = GINPUT(N) returns a third result, BUTTON, that -% contains a vector of integers specifying which mouse button was -% used (1,2,3 from left) or ASCII numbers if a key on the keyboard -% was used. -% -% Usage: -% [out1,out2,out3] = ginput(arg1) - -% Copyright 1984-2005 The MathWorks, Inc. -% $Revision: 1.1 $ $Date: 2009/04/03 22:56:26 $ - -out1 = []; out2 = []; out3 = []; y = []; -c = computer; -if ~strcmp(c(1:2),'PC') - tp = get(0,'TerminalProtocol'); -else - tp = 'micro'; -end - -if ~strcmp(tp,'none') && ~strcmp(tp,'x') && ~strcmp(tp,'micro'), - if nargout == 1, - if nargin == 1, - out1 = trmginput(arg1); - else - out1 = trmginput; - end - elseif nargout == 2 || nargout == 0, - if nargin == 1, - [out1,out2] = trmginput(arg1); - else - [out1,out2] = trmginput; - end - if nargout == 0 - out1 = [ out1 out2 ]; - end - elseif nargout == 3, - if nargin == 1, - [out1,out2,out3] = trmginput(arg1); - else - [out1,out2,out3] = trmginput; - end - end -else - - fig = gcf; - figure(gcf); - - if nargin == 0 - how_many = -1; - b = []; - else - how_many = arg1; - b = []; - if ischar(how_many) ... - || size(how_many,1) ~= 1 || size(how_many,2) ~= 1 ... - || ~(fix(how_many) == how_many) ... - || how_many < 0 - error('MATLAB:ginput:NeedPositiveInt', 'Requires a positive integer.') - end - if how_many == 0 - ptr_fig = 0; - while(ptr_fig ~= fig) - ptr_fig = get(0,'PointerWindow'); - end - scrn_pt = get(0,'PointerLocation'); - loc = get(fig,'Position'); - pt = [scrn_pt(1) - loc(1), scrn_pt(2) - loc(2)]; - out1 = pt(1); y = pt(2); - elseif how_many < 0 - error('MATLAB:ginput:InvalidArgument', 'Argument must be a positive integer.') - end - end - - % Suspend figure functions - state = uisuspend(fig); - - toolbar = findobj(allchild(fig),'flat','Type','uitoolbar'); - if ~isempty(toolbar) - ptButtons = [uigettool(toolbar,'Plottools.PlottoolsOff'), ... - uigettool(toolbar,'Plottools.PlottoolsOn')]; - ptState = get (ptButtons,'Enable'); - set (ptButtons,'Enable','off'); - end - - set(fig,'pointer','fullcrosshair'); - fig_units = get(fig,'units'); - char = 0; - - % We need to pump the event queue on unix - % before calling WAITFORBUTTONPRESS - drawnow - - while how_many ~= 0 - % Use no-side effect WAITFORBUTTONPRESS - waserr = 0; - try - keydown = wfbp; - catch - waserr = 1; - end - if(waserr == 1) - if(ishandle(fig)) - set(fig,'units',fig_units); - uirestore(state); - error('MATLAB:ginput:Interrupted', 'Interrupted'); - else - error('MATLAB:ginput:FigureDeletionPause', 'Interrupted by figure deletion'); - end - end - - ptr_fig = get(0,'CurrentFigure'); - if(ptr_fig == fig) - if keydown - char = get(fig, 'CurrentCharacter'); - button = abs(get(fig, 'CurrentCharacter')); - scrn_pt = get(0, 'PointerLocation'); - set(fig,'units','pixels') - loc = get(fig, 'Position'); - pt = [scrn_pt(1) - loc(1), scrn_pt(2) - loc(2)]; - set(fig,'CurrentPoint',pt); - else - button = get(fig, 'SelectionType'); - if strcmp(button,'open') - button = 1; - elseif strcmp(button,'normal') - button = 1; - elseif strcmp(button,'extend') - button = 2; - elseif strcmp(button,'alt') - button = 3; - else - error('MATLAB:ginput:InvalidSelection', 'Invalid mouse selection.') - end - end - pt = get(gca, 'CurrentPoint'); - - how_many = how_many - 1; - - if(char == 13) % & how_many ~= 0) - % if the return key was pressed, char will == 13, - % and that's our signal to break out of here whether - % or not we have collected all the requested data - % points. - % If this was an early breakout, don't include - % the key info in the return arrays. - % We will no longer count it if it's the last input. - break; - end - - out1 = [out1;pt(1,1)]; - y = [y;pt(1,2)]; - b = [b;button]; - end - end - - uirestore(state); - if ~isempty(toolbar) && ~isempty(ptButtons) - set (ptButtons(1),'Enable',ptState{1}); - set (ptButtons(2),'Enable',ptState{2}); - end - set(fig,'units',fig_units); - - if nargout > 1 - out2 = y; - if nargout > 2 - out3 = b; - end - else - out1 = [out1 y]; - end - - line(out1,y); - line([out1(length(out1)) out1(1)],[y(length(y)) y(1)]); - -end - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -function key = wfbp -%WFBP Replacement for WAITFORBUTTONPRESS that has no side effects. - -fig = gcf; -current_char = []; - -% Now wait for that buttonpress, and check for error conditions -waserr = 0; -try - h=findall(fig,'type','uimenu','accel','C'); % Disabling ^C for edit menu so the only ^C is for - set(h,'accel',''); % interrupting the function. - keydown = waitforbuttonpress; - current_char = double(get(fig,'CurrentCharacter')); % Capturing the character. - if~isempty(current_char) && (keydown == 1) % If the character was generated by the - if(current_char == 3) % current keypress AND is ^C, set 'waserr'to 1 - waserr = 1; % so that it errors out. - end - end - - set(h,'accel','C'); % Set back the accelerator for edit menu. -catch - waserr = 1; -end -drawnow; -if(waserr == 1) - set(h,'accel','C'); % Set back the accelerator if it errored out. - error('MATLAB:ginput:Interrupted', 'Interrupted'); -end - -if nargout>0, key = keydown; end -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Index: /issm/trunk-jpl/src/m/geometry/SegIntersect.m =================================================================== --- /issm/trunk-jpl/src/m/geometry/SegIntersect.m (revision 13009) +++ /issm/trunk-jpl/src/m/geometry/SegIntersect.m (revision 13010) @@ -77,5 +77,3 @@ bool=0; return; - end - Index: sm/trunk-jpl/src/m/interp/FieldFindVarNames.m =================================================================== --- /issm/trunk-jpl/src/m/interp/FieldFindVarNames.m (revision 13009) +++ (revision ) @@ -1,153 +1,0 @@ -function Names=FieldFindVarNames(filename) -%FIELDFINDVARNAMES - find names of variables in a data set file -% -% This routines looks at the variables contained in a file and finds out -% the names of the variables that are needed for an interpolation (x,y,data) -% or (index,x,y,data) -% -% Usage: -% Names=FieldFindVarNames(filename) -% -% Example: -% Names=FieldFindVarNames('thickness.mat') -% -% See also: INTERPFROMFILE, GRIDDATA - -%some checks -if nargin~=1 | nargout~=1 - help FieldFindVarNames - error('FieldFindVarNames error message: bad usage'); -end -if ~exist(filename) - error(['FieldFindVarNames error message: file ' filename ' does not exist']); -end - -%Get variables -A=whos('-file',filename); - -%find x,y,vx and vy -xenum=NaN; yenum=NaN; dataenum=NaN; indexenum=NaN; -if length(A)==3, - isnode=1; - for i=1:3 - if strcmpi(A(i).name(1),'x'); - xenum=i; - elseif strcmpi(A(i).name(1),'y'); - yenum=i; - elseif (strncmpi(A(i).name,filename,3) | strncmpi(A(i).name,'data',4)), - dataenum=i; - else - %nothing - end - end -elseif length(A)==4, - isnode=0; - for i=1:4 - if strcmpi(A(i).name(1),'x'); - xenum=i; - elseif strcmpi(A(i).name(1),'y'); - yenum=i; - elseif (strncmpi(A(i).name,'index',5) | strncmpi(A(i).name,'elements',7)); - indexenum=i; - elseif (strncmpi(A(i).name,filename,3) | strncmpi(A(i).name,'data',4)), - dataenum=i; - else - %nothing - end - end -else - error(['FieldFindVarNames error message: file ' filename ' not supported yet (it should hold 3 variables x,y and data (for nodes) OR 4 variables x,y,index and data (for mesh))']); -end - -%2: if only one item is missing, find it by elimination -if ~isnode, - pos=find(isnan([xenum yenum indexenum dataenum])); - if length(pos)==1, - list=[xenum yenum indexenum dataenum]; list(pos)=[]; - if pos==1, - xenum=setdiff(1:4,list); - elseif pos==2, - yenum=setdiff(1:4,list); - elseif pos==3, - indexenum=setdiff(1:4,list); - elseif pos==4, - dataenum=setdiff(1:4,list); - end - end -else - pos=find(isnan([xenum yenum dataenum])); - if length(pos)==1, - list=[xenum yenum indexenum dataenum]; list(pos)=[]; - if pos==1, - xenum=setdiff(1:3,list); - elseif pos==2, - yenum=setdiff(1:3,list); - elseif pos==3, - dataenum=setdiff(1:3,list); - end - end -end - -%assum that we have found at least xenum and yenum -if ( isnan(xenum) | isnan(yenum)) - error(['FieldFindVarNames error message: file ' filename ' not supported yet (the coordinates vectors should be named x and y)']); -end - -%find index -if (~isnode & isnan(indexenum)), - for i=1:4 - lengthi=min(A(i).size); - if (lengthi==3), - indexenum=i; - end - end - if isnan(indexenum), - error(['FieldFindVarNames error message: file ' filename ' not supported yet (index not found)']); - end -end - -%4: last chance -if ~isnode, - pos=find(isnan([xenum yenum indexenum dataenum])); - if length(pos)==1, - list=[xenum yenum indexenum dataenum]; list(pos)=[]; - if pos==1, - xenum=setdiff(1:4,list); - elseif pos==2, - yenum=setdiff(1:4,list); - elseif pos==3, - indexenum=setdiff(1:4,list); - elseif pos==4, - dataenum=setdiff(1:4,list); - end - end -else - pos=find(isnan([xenum yenum dataenum])); - if length(pos)==1, - list=[xenum yenum indexenum dataenum]; list(pos)=[]; - if pos==1, - xenum=setdiff(1:3,list); - elseif pos==2, - yenum=setdiff(1:3,list); - elseif pos==3, - dataenum=setdiff(1:3,list); - end - end -end - -%last check -if isnan(dataenum) - error(['FieldFindVarNames error message: file ' filename ' not supported yet (data not found)']); -end - -%create output -Names=struct(); -Names.xname=A(xenum).name; -Names.yname=A(yenum).name; -Names.dataname=A(dataenum).name; -if ~isnode, - Names.indexname=A(indexenum).name; - Names.interp='mesh'; -else - Names.interp='node'; -end Index: sm/trunk-jpl/src/m/interp/FillHole.m =================================================================== --- /issm/trunk-jpl/src/m/interp/FillHole.m (revision 13009) +++ (revision ) @@ -1,39 +1,0 @@ -function field=FillHole(index,x,y,field) -%FILLHOLE - fill mesh data that has holes in it (hole is defined by field==NaN). Get the nearest neighboors to fill the holes. -% -% Usage: -% surface=FillHole(index,x,y,surface) -% -% Example: -% md.geometry.surface=FillHole(md.mesh.elements,x,md.mesh.y,md.geometry.surface) -% - -%some checks -if nargin~=4 | nargout~=1 - help FillHole - error('FillHole error message: bad usage'); -end - -if length(x)~=length(y), - error('plugdata error message: x and y should have the same length'); -end - - -pos_hole=find(isnan(field)); -pos_full=find(~isnan(field)); - -%reduce our field to not include any holes: -field_noholes=field(pos_full); -x_noholes=x(pos_full); -y_noholes=y(pos_full); - -for i=1:length(pos_hole) - - if (mod(i,100)==0), - fprintf('\b\b\b\b\b\b\b%5.2f%s',i/length(pos_hole)*100,' %'); - end - - %search the node on the closest to i, that is not in a hole either - [d posd]=min(sqrt((x(pos_hole(i))-x_noholes).^2+(y(pos_hole(i))-y_noholes).^2)); - field(pos_hole(i))=field_noholes(posd); -end Index: /issm/trunk-jpl/src/m/interp/InterpFromFile.m =================================================================== --- /issm/trunk-jpl/src/m/interp/InterpFromFile.m (revision 13009) +++ /issm/trunk-jpl/src/m/interp/InterpFromFile.m (revision 13010) @@ -43,2 +43,158 @@ data_out=InterpFromMeshToMesh2d(Data.(Names.indexname),Data.(Names.xname),Data.(Names.yname),Data.(Names.dataname),x,y); end + +end +function Names=FieldFindVarNames(filename) +%FIELDFINDVARNAMES - find names of variables in a data set file +% +% This routines looks at the variables contained in a file and finds out +% the names of the variables that are needed for an interpolation (x,y,data) +% or (index,x,y,data) +% +% Usage: +% Names=FieldFindVarNames(filename) +% +% Example: +% Names=FieldFindVarNames('thickness.mat') +% +% See also: INTERPFROMFILE, GRIDDATA + +%some checks +if nargin~=1 | nargout~=1 + help FieldFindVarNames + error('FieldFindVarNames error message: bad usage'); +end +if ~exist(filename) + error(['FieldFindVarNames error message: file ' filename ' does not exist']); +end + +%Get variables +A=whos('-file',filename); + +%find x,y,vx and vy +xenum=NaN; yenum=NaN; dataenum=NaN; indexenum=NaN; +if length(A)==3, + isnode=1; + for i=1:3 + if strcmpi(A(i).name(1),'x'); + xenum=i; + elseif strcmpi(A(i).name(1),'y'); + yenum=i; + elseif (strncmpi(A(i).name,filename,3) | strncmpi(A(i).name,'data',4)), + dataenum=i; + else + %nothing + end + end +elseif length(A)==4, + isnode=0; + for i=1:4 + if strcmpi(A(i).name(1),'x'); + xenum=i; + elseif strcmpi(A(i).name(1),'y'); + yenum=i; + elseif (strncmpi(A(i).name,'index',5) | strncmpi(A(i).name,'elements',7)); + indexenum=i; + elseif (strncmpi(A(i).name,filename,3) | strncmpi(A(i).name,'data',4)), + dataenum=i; + else + %nothing + end + end +else + error(['FieldFindVarNames error message: file ' filename ' not supported yet (it should hold 3 variables x,y and data (for nodes) OR 4 variables x,y,index and data (for mesh))']); +end + +%2: if only one item is missing, find it by elimination +if ~isnode, + pos=find(isnan([xenum yenum indexenum dataenum])); + if length(pos)==1, + list=[xenum yenum indexenum dataenum]; list(pos)=[]; + if pos==1, + xenum=setdiff(1:4,list); + elseif pos==2, + yenum=setdiff(1:4,list); + elseif pos==3, + indexenum=setdiff(1:4,list); + elseif pos==4, + dataenum=setdiff(1:4,list); + end + end +else + pos=find(isnan([xenum yenum dataenum])); + if length(pos)==1, + list=[xenum yenum indexenum dataenum]; list(pos)=[]; + if pos==1, + xenum=setdiff(1:3,list); + elseif pos==2, + yenum=setdiff(1:3,list); + elseif pos==3, + dataenum=setdiff(1:3,list); + end + end +end + +%assum that we have found at least xenum and yenum +if ( isnan(xenum) | isnan(yenum)) + error(['FieldFindVarNames error message: file ' filename ' not supported yet (the coordinates vectors should be named x and y)']); +end + +%find index +if (~isnode & isnan(indexenum)), + for i=1:4 + lengthi=min(A(i).size); + if (lengthi==3), + indexenum=i; + end + end + if isnan(indexenum), + error(['FieldFindVarNames error message: file ' filename ' not supported yet (index not found)']); + end +end + +%4: last chance +if ~isnode, + pos=find(isnan([xenum yenum indexenum dataenum])); + if length(pos)==1, + list=[xenum yenum indexenum dataenum]; list(pos)=[]; + if pos==1, + xenum=setdiff(1:4,list); + elseif pos==2, + yenum=setdiff(1:4,list); + elseif pos==3, + indexenum=setdiff(1:4,list); + elseif pos==4, + dataenum=setdiff(1:4,list); + end + end +else + pos=find(isnan([xenum yenum dataenum])); + if length(pos)==1, + list=[xenum yenum indexenum dataenum]; list(pos)=[]; + if pos==1, + xenum=setdiff(1:3,list); + elseif pos==2, + yenum=setdiff(1:3,list); + elseif pos==3, + dataenum=setdiff(1:3,list); + end + end +end + +%last check +if isnan(dataenum) + error(['FieldFindVarNames error message: file ' filename ' not supported yet (data not found)']); +end + +%create output +Names=struct(); +Names.xname=A(xenum).name; +Names.yname=A(yenum).name; +Names.dataname=A(dataenum).name; +if ~isnode, + Names.indexname=A(indexenum).name; + Names.interp='mesh'; +else + Names.interp='node'; +end +end Index: sm/trunk-jpl/src/m/interp/VelFindVarNames.m =================================================================== --- /issm/trunk-jpl/src/m/interp/VelFindVarNames.m (revision 13009) +++ (revision ) @@ -1,122 +1,0 @@ -function Names=VelFindVarNames(filename) -%VELFINDVARNAMES - find names of variables in a velocity data set file -% -% This routines looks at the variables contained in a file and finds out -% the names of the variables that are needed for an interpolation (x,y,vx,vy) -% or (index,x,y,vx,vy) -% -% Usage: -% Names=VelFindVarNames(filename) -% -% Example: -% Names=VelFindVarNames('velocities.mat') -% -% See also: INTERPFROMFILE, GRIDDATA - -%some checks -if nargin~=1 | nargout~=1 - help VelFindVarNames - error('VelFindVarNames error message: bad usage'); -end -if ~exist(filename) - error(['VelFindVarNames error message: file ' filename ' does not exist']); -end - -%Get variables -A=whos('-file',filename); - -%find x,y,vx and vy -xenum=NaN; yenum=NaN; vxenum=NaN; vyenum=NaN; indexenum=NaN; -if length(A)==4, - isnode=1; - for i=1:4 - if strcmpi(A(i).name(1),'x'); - xenum=i; - elseif strcmpi(A(i).name(1),'y'); - yenum=i; - else - if (strcmpi(A(i).name(end),'x') | strncmpi(A(i).name,'vx',2)); - vxenum=i; - elseif (strcmpi(A(i).name(end),'y') | strncmpi(A(i).name,'vy',2)); - vyenum=i; - end - end - end -elseif length(A)==5, - isnode=0; - for i=1:5 - if strcmpi(A(i).name(1),'x'); - xenum=i; - elseif strcmpi(A(i).name(1),'y'); - yenum=i; - elseif (strcmpi(A(i).name(1),'index') | strcmpi(A(i).name(1),'elements')); - indexenum=i; - else - if (strcmpi(A(i).name(end),'x') | strncmpi(A(i).name,'vx',2)); - vxenum=i; - elseif (strcmpi(A(i).name(end),'y') | strncmpi(A(i).name,'vy',2)); - vyenum=i; - end - end - end -else - error(['VelFindVarNames error message: file ' filename ' not supported yet (it should hold 4 variables x,y,vx and vy (for nodes) OR 5 variables x,y,index,vx and vy (for mesh))']); -end - -%assum that we have found at least vxenum and vyenum -if ( isnan(vxenum) | isnan(vyenum)) - error(['VelFindVarNames error message: file ' filename ' not supported yet (the velocities should be named vx and vy)']); -end - -%find index -if (~isnode & isnan(indexenum)), - for i=1:5 - lengthi=min(A(i).size); - if (lengthi==3), - indexenum=i; - end - end - if isnan(indexenum), - error(['VelFindVarNames error message: file ' filename ' not supported yet (index not found)']); - end -end - -%find x y -if (isnan(xenum) | isnan(yenum)) - - %check the size - if A(vxenum).size(1)==A(vxenum).size(2), - error(['VelFindVarNames error message: file ' filename ' not supported (velocities is a square matrix, save x and y with another name)']); - end - if ~(A(vxenum).size(1)==A(vyenum).size(1) & A(vxenum).size(2)==A(vyenum).size(2)), - error(['VelFindVarNames error message: file ' filename ' not supported (vx and vy matrices do not have the same size)']); - end - - %find xenum and yenum - for i=1:4 - lengthi=max(A(i).size); - if ((i~=vxenum) & (lengthi==A(vxenum).size(1) | lengthi==A(vxenum).size(1)+1)), - yenum=i; - elseif ((i~=vxenum) & (lengthi==A(vxenum).size(2) | lengthi==A(vxenum).size(2)+1)), - xenum=i; - end - end - - %last check - if (isnan(xenum) | isnan(yenum)) - error(['plugdata error message: file ' filename ' not supported yet']); - end -end - -%create output -Names=struct(); -Names.xname=A(xenum).name; -Names.yname=A(yenum).name; -Names.vxname=A(vxenum).name; -Names.vyname=A(vyenum).name; -if ~isnode, - Names.indexname=A(indexenum).name; - Names.interp='mesh'; -else - Names.interp='node'; -end Index: /issm/trunk-jpl/src/m/interp/plugvelocities.m =================================================================== --- /issm/trunk-jpl/src/m/interp/plugvelocities.m (revision 13009) +++ /issm/trunk-jpl/src/m/interp/plugvelocities.m (revision 13010) @@ -14,4 +14,5 @@ % See also: INTERPFROMFILE, GRIDDATA +disp('WARNING: deprecated functions (plugvelocities)'); %some checks if nargin~=3 | nargout~=1 @@ -40,2 +41,126 @@ md.initialization.vy=md.inversion.vy_obs; md.initialization.vel=md.inversion.vel_obs; +end + +function Names=VelFindVarNames(filename) +%VELFINDVARNAMES - find names of variables in a velocity data set file +% +% This routines looks at the variables contained in a file and finds out +% the names of the variables that are needed for an interpolation (x,y,vx,vy) +% or (index,x,y,vx,vy) +% +% Usage: +% Names=VelFindVarNames(filename) +% +% Example: +% Names=VelFindVarNames('velocities.mat') +% +% See also: INTERPFROMFILE, GRIDDATA + +%some checks +if nargin~=1 | nargout~=1 + help VelFindVarNames + error('VelFindVarNames error message: bad usage'); +end +if ~exist(filename) + error(['VelFindVarNames error message: file ' filename ' does not exist']); +end + +%Get variables +A=whos('-file',filename); + +%find x,y,vx and vy +xenum=NaN; yenum=NaN; vxenum=NaN; vyenum=NaN; indexenum=NaN; +if length(A)==4, + isnode=1; + for i=1:4 + if strcmpi(A(i).name(1),'x'); + xenum=i; + elseif strcmpi(A(i).name(1),'y'); + yenum=i; + else + if (strcmpi(A(i).name(end),'x') | strncmpi(A(i).name,'vx',2)); + vxenum=i; + elseif (strcmpi(A(i).name(end),'y') | strncmpi(A(i).name,'vy',2)); + vyenum=i; + end + end + end +elseif length(A)==5, + isnode=0; + for i=1:5 + if strcmpi(A(i).name(1),'x'); + xenum=i; + elseif strcmpi(A(i).name(1),'y'); + yenum=i; + elseif (strcmpi(A(i).name(1),'index') | strcmpi(A(i).name(1),'elements')); + indexenum=i; + else + if (strcmpi(A(i).name(end),'x') | strncmpi(A(i).name,'vx',2)); + vxenum=i; + elseif (strcmpi(A(i).name(end),'y') | strncmpi(A(i).name,'vy',2)); + vyenum=i; + end + end + end +else + error(['VelFindVarNames error message: file ' filename ' not supported yet (it should hold 4 variables x,y,vx and vy (for nodes) OR 5 variables x,y,index,vx and vy (for mesh))']); +end + +%assum that we have found at least vxenum and vyenum +if ( isnan(vxenum) | isnan(vyenum)) + error(['VelFindVarNames error message: file ' filename ' not supported yet (the velocities should be named vx and vy)']); +end + +%find index +if (~isnode & isnan(indexenum)), + for i=1:5 + lengthi=min(A(i).size); + if (lengthi==3), + indexenum=i; + end + end + if isnan(indexenum), + error(['VelFindVarNames error message: file ' filename ' not supported yet (index not found)']); + end +end + +%find x y +if (isnan(xenum) | isnan(yenum)) + + %check the size + if A(vxenum).size(1)==A(vxenum).size(2), + error(['VelFindVarNames error message: file ' filename ' not supported (velocities is a square matrix, save x and y with another name)']); + end + if ~(A(vxenum).size(1)==A(vyenum).size(1) & A(vxenum).size(2)==A(vyenum).size(2)), + error(['VelFindVarNames error message: file ' filename ' not supported (vx and vy matrices do not have the same size)']); + end + + %find xenum and yenum + for i=1:4 + lengthi=max(A(i).size); + if ((i~=vxenum) & (lengthi==A(vxenum).size(1) | lengthi==A(vxenum).size(1)+1)), + yenum=i; + elseif ((i~=vxenum) & (lengthi==A(vxenum).size(2) | lengthi==A(vxenum).size(2)+1)), + xenum=i; + end + end + + %last check + if (isnan(xenum) | isnan(yenum)) + error(['plugdata error message: file ' filename ' not supported yet']); + end +end + +%create output +Names=struct(); +Names.xname=A(xenum).name; +Names.yname=A(yenum).name; +Names.vxname=A(vxenum).name; +Names.vyname=A(vyenum).name; +if ~isnode, + Names.indexname=A(indexenum).name; + Names.interp='mesh'; +else + Names.interp='node'; +end Index: /issm/trunk-jpl/src/m/inversions/MisfitDeinterlace.m =================================================================== --- /issm/trunk-jpl/src/m/inversions/MisfitDeinterlace.m (revision 13010) +++ /issm/trunk-jpl/src/m/inversions/MisfitDeinterlace.m (revision 13010) @@ -0,0 +1,21 @@ +function Jstruct=MisfitDeinterlace(misfit,type) +%MISFITDEINTERLACE - deinterlace misfits that are mixed together, using type. +% +% Usage: +% Jstruct=MisfitDeinterlace(misfit,type) +% +% Example: +% Jstruct=MisfitDeinterlace(md.results.diagnostic.J,md.fit) +% +% +Jstruct=struct(); + +count=1; +for i=0:max(type), + pos=find(type==i); + if length(pos), + Jstruct(count).type=i; + Jstruct(count).J=misfit(pos); + count=count+1; + end +end Index: sm/trunk-jpl/src/m/latlong/pargenerate.m =================================================================== --- /issm/trunk-jpl/src/m/latlong/pargenerate.m (revision 13009) +++ (revision ) @@ -1,45 +1,0 @@ -function pargenerate(filename,xm,ym,data,varargin) -%PARGENERATE - generate parameter file for images ala Rignot -% -% Usage: -% pargenerate(filename,data,xm,ym,options); -% -% Supported options: -% - title: dataset title -% - latitude: standard latitude (degree) -% - meridian: meridian (degree) -% - format: binary format - -%process options -options=pairoptions(varargin{:}); -Title=getfieldvalue(options,'title','N/A'); -latitude=getfieldvalue(options,'latitude','N/A'); -meridian=getfieldvalue(options,'meridian','N/A'); -format=getfieldvalue(options,'format','single'); - -%Get data info -[nlines ncols]=size(data); -xmin=min(xm); -ymax=max(ym); -postx=abs(xm(2)-xm(1)); -posty=abs(ym(2)-ym(1)); - -%Open header file and get machine type -fid=fopen(filename,'wt'); -[filename, permission, machineformat, encoding] = fopen(fid); - -%write header file -fprintf(fid,'%s\n','ISSM gridded dataset parameter file'); -fprintf(fid,'%s%s\n','title: ',Title); -fprintf(fid,'%s \n','DEM_projection: PS'); -fprintf(fid,'%s%s\n','data_format: ',format); -fprintf(fid,'%s%s\n','endian: ',machineformat); -fprintf(fid,'%s%d\n','width: ',ncols); -fprintf(fid,'%s%d\n','nlines: ',nlines); -fprintf(fid,'%s%g%s\n','PS_secant_lat: ',latitude,' decimal degrees'); -fprintf(fid,'%s%g%s\n','PS_meridian_long: ',meridian,' decimal degrees'); -fprintf(fid,'%s%-15.10g\n','PS_corner_north: ',ymax); -fprintf(fid,'%s%-15.10g\n','PS_corner_east: ',xmin); -fprintf(fid,'%s%g%s\n','PS_post_north: ',postx,' m'); -fprintf(fid,'%s%g%s\n','PS_post_east: ',posty,' m'); -fclose(fid); Index: /issm/trunk-jpl/src/m/mesh/roundmesh.m =================================================================== --- /issm/trunk-jpl/src/m/mesh/roundmesh.m (revision 13009) +++ /issm/trunk-jpl/src/m/mesh/roundmesh.m (revision 13010) @@ -20,5 +20,5 @@ y_list=radius*y_list.*sin(theta); A=struct('x',x_list,'y',y_list,'density',1); -expgen('RoundDomainOutline.exp',A,1); +expwrite('RoundDomainOutline.exp',A); %Call Bamg Index: /issm/trunk-jpl/src/m/miscellaneous/fielddisplay.m =================================================================== --- /issm/trunk-jpl/src/m/miscellaneous/fielddisplay.m (revision 13010) +++ /issm/trunk-jpl/src/m/miscellaneous/fielddisplay.m (revision 13010) @@ -0,0 +1,154 @@ +function fielddisplay(md,name,comment) +%FIELDDISPLAY - display model field +% +% Usage: +% fielddisplay(md,offset,name,comment) + + %get field + field=md.(name); + + %disp corresponding line as a function of field type (offset set as 9 spaces) + parsedisplay(' ',name,field,comment); + +end %function + +function parsedisplay(offset,name,field,comment); %{{{ + + %string + if ischar(field), + + if length(field)>30; + displayunit(offset,name,'not displayed',comment), + else + displayunit(offset,name,['''' field ''''],comment), + end + + %numeric + elseif isnumeric(field) + + %get size + fieldsize=size(field); + + %double + if max(fieldsize)==1, + displayunit(offset,name,num2str(field),comment), + %matrix + else + displayunit(offset,name,['(' num2str(fieldsize(1)) 'x' num2str(fieldsize(2)) ')'],comment), + end + + %logical + elseif islogical(field) + + %get size + fieldsize=size(field); + + %single value + if max(fieldsize)==1, + if (field) + displayunit(offset,name,'true',comment), + else + displayunit(offset,name,'false',comment), + end + %matrix + else + displayunit(offset,name,['(' num2str(fieldsize(1)) 'x' num2str(fieldsize(2)) ')'],comment), + end + + %structure + elseif isstruct(field), + if ~isempty(fields(field)) + displayunit(offset,name,'(structure)',comment), + struct_display(field,[offset ' ']), + else + displayunit(offset,name,'N/A',comment), + end + + %cell + elseif iscell(field), + cell_display(offset,name,field,comment), + + else + displayunit(offset,name,'not displayed',comment), + + end +end%}}} + +function struct_display(structure,offset) % {{{ + + structure_fields=fields(structure); + + for i=1:length(structure_fields), + + %get current field + field=structure.(structure_fields{i}); + + %recursive call if necessary + if isstruct(field), + displayunit(offset,structure_fields{i},'(structure)',''), + struct_display(field,[offset ' ']); + + %display value + else + parsedisplay(offset,structure_fields{i},field,''); + end + end +end% }}} +function cell_display(offset,name,field,comment) % {{{ + + %initialization + string='{'; + + %go through the cell and fill string + if length(field)<5; + for i=1:length(field), + if ischar(field{i}), + string=[string '''' field{i} ''',']; + elseif (isnumeric(field{i}) & length(field{i})==1) + string=[string num2str(field{i}) ',' ]; + else + string='{'; + break + end + end + end + if strcmp(string,'{'), + string=['(' num2str(size(field,1)) 'x' num2str(size(field,2)) ')']; + else + string=[string(1:end-1) '}']; + end + + %call displayunit + displayunit(offset,name,string,comment); +end% }}} +function displayunit(offset,name,characterization,comment),% {{{ + + %take care of name + if length(name)>23, + name=[name(1:20) '...']; + end + + %take care of characterization + if (strcmp(characterization,['''' '''']) | strcmp(characterization,'NaN')), + characterization='N/A'; + end + if length(characterization)>15, + characterization=[characterization(1:12) '...']; + end + + %print + if isempty(comment) + disp(sprintf('%s%-23s: %-15s',offset,name,characterization)); + else + if ischar(comment), + disp(sprintf('%s%-23s: %-15s -- %s',offset,name,characterization,comment)); + elseif iscell(comment), + disp(sprintf('%s%-23s: %-15s -- %s',offset,name,characterization,comment{1})); + for i=2:length(comment), + disp(sprintf('%s%-23s %-15s %s',offset,'','',comment{i})); + end + else + error('fielddisplay error message: format for comment not supportet yet'); + end + end +end% }}} Index: /issm/trunk-jpl/src/m/miscellaneous/parallelrange.m =================================================================== --- /issm/trunk-jpl/src/m/miscellaneous/parallelrange.m (revision 13010) +++ /issm/trunk-jpl/src/m/miscellaneous/parallelrange.m (revision 13010) @@ -0,0 +1,23 @@ +function [i1,i2]=parallelrange(rank,numprocs,globalsize) +%PARALLELRANGE - from a rank, and a number of processors, figure out a range, for parallel tasks. +% +% Usage: +% [i1,i1]=parallelrange(rank,numprocs,globalsize) + +num_local_rows=zeros(numprocs,1); + +for i=1:numprocs, + %we use floor. we under distribute rows. The rows left are then redistributed, therefore resulting in a more even distribution. + num_local_rows(i)=floor(globalsize/numprocs); +end + + +%There may be some rows left. Distribute evenly. +row_rest=globalsize - numprocs*floor(globalsize/numprocs); + +for i=1:row_rest, + num_local_rows(i)=num_local_rows(i)+1; +end + +i1=sum(num_local_rows(1:rank-1))+1; +i2=i1+num_local_rows(rank)-1; Index: /issm/trunk-jpl/src/m/miscellaneous/parallelrange.py =================================================================== --- /issm/trunk-jpl/src/m/miscellaneous/parallelrange.py (revision 13010) +++ /issm/trunk-jpl/src/m/miscellaneous/parallelrange.py (revision 13010) @@ -0,0 +1,25 @@ +#! /usr/bin/env python +def parallelrange(rank,numprocs,globalsize): + """ + PARALLELRANGE - from a rank, and a number of processors, figure out a range, for parallel tasks. + + Usage: + i1,i2=parallelrange(rank,numprocs,globalsize) + """ + + #We use floor. we under distribute rows. The rows left are then redistributed, therefore resulting in a more even distribution. + num_local_rows=[int(globalsize/numprocs) for i in xrange(numprocs)] + + #There may be some rows left. Distribute evenly. + row_rest=globalsize - numprocs*int(globalsize/numprocs) + + for i in xrange(row_rest): + num_local_rows[i]=num_local_rows[i]+1 + + i1=0 + for i in xrange(rank-1): + i1+=num_local_rows[i] + i2=i1+num_local_rows[rank-1]-1 + + return i1,i2 + Index: /issm/trunk-jpl/src/m/regional/basinzoom.m =================================================================== --- /issm/trunk-jpl/src/m/regional/basinzoom.m (revision 13010) +++ /issm/trunk-jpl/src/m/regional/basinzoom.m (revision 13010) @@ -0,0 +1,94 @@ +function varargout=basinzoom(varargin) +%ANTZOOM - zoom on a basin in Antarctica or Greenland. +% +% This function zooms on an existing figure describing Antarctica or Greenland +% The zooming depends on the region name provided as input. +% +% Usage: +% varargout=basinzoom(options) + +%recover some options, and set defaults + +%is varargin an options database already? +if nargin==0, + options=pairoptions(varargin{:}); +elseif (isa(varargin{1},'plotoptions') | isa(varargin{1},'pairoptions')), + %do nothing to the varargin: + options=varargin{1}; +else + %process varargin for options: + options=pairoptions(varargin{:}); +end + +unitmultiplier=getfieldvalue(options,'unit',NaN); +basin=getfieldvalue(options,'basin'); + +if exist(options,'basindelta'), + + basindeltax=getfieldvalue(options,'basindelta',300); + basindeltay=getfieldvalue(options,'basindelta',300); +else + basindeltax=getfieldvalue(options,'basindeltax',300); + basindeltay=getfieldvalue(options,'basindeltay',300); +end + +%multiply by 1000 to get kms +basindeltax=basindeltax*1000; +basindeltay=basindeltay*1000; + +%Ok, find basin we are talking about: +load([jplsvn() '/projects/ModelData/Names/Names.mat']); + +%Go through names: +found=0; +for i=1:size(names,1), + if strcmpi(names{i,1},basin), + %ok, we've got the region. Get lat and long: + long=names{i,2}; + lat=names{i,3}; + hemisphere=names{i,4}; + found=1; + break; + end +end + +if ~found, + error(['basinzoom error message: cannot find basin ' basin '. Use isbasin to determine a basin name.']); +end + +if hemisphere==+1, + central_meridian=getfieldvalue(options,'central_meridian',45); + standard_parallel=getfieldvalue(options,'standard_parallel',70); +else + central_meridian=getfieldvalue(options,'central_meridian',0); + standard_parallel=getfieldvalue(options,'standard_parallel',71); +end + +%Transform lat long into x,y: +[xc,yc]=ll2xy(lat,long,hemisphere,central_meridian,standard_parallel); + +%compute x0,x1 and y0,y1 using basindeltax and basindeltay +x0=xc-basindeltax/2; +x1=xc+basindeltax/2; +y0=yc-basindeltay/2; +y1=yc+basindeltay/2; + +if ~isnan(unitmultiplier) + x0=x0*unitmultiplier; + x1=x1*unitmultiplier; + y0=y0*unitmultiplier; + y1=y1*unitmultiplier; +end + +%if output arguments are present, return the limits, +%otherwise, set them on the current graphic. +if nargout==2, + found=1; + varargout{1}=[x0 x1]; + varargout{2}=[y0 y1]; +else + xlim([x0 x1]); + ylim([y0 y1]); + found=1; + daspect([1;1;1]); +end Index: /issm/trunk-jpl/src/m/regional/isbasin.m =================================================================== --- /issm/trunk-jpl/src/m/regional/isbasin.m (revision 13010) +++ /issm/trunk-jpl/src/m/regional/isbasin.m (revision 13010) @@ -0,0 +1,17 @@ +function isbasin(name) +%ISBASIN: figure out if a basin name exists. +% +% +% Usage: index=isbasin('jks'); +% +% + +%First, load basin names: +load([jplsvn '/projects/ModelData/Names/Names.mat']); + +%go through names: +for i=1:length(names), + if ~isempty(strfind(names{i,1},name)), + disp(['''' names{i,1} ''' Long:' num2str(names{i,2}) ' Lat:' num2str(names{i,3}) ]); + end +end Index: /issm/trunk-jpl/src/m/regional/plotbasins.m =================================================================== --- /issm/trunk-jpl/src/m/regional/plotbasins.m (revision 13010) +++ /issm/trunk-jpl/src/m/regional/plotbasins.m (revision 13010) @@ -0,0 +1,18 @@ +%display all the domain outlines in a directory + +basins=listfiles; + +hold on +for i=1:length(basins), + %check whether this is a .exp file + basin=basins{i}; + if strcmpi(basin(end-3:end),'.exp'), + + contour=expread(basin,0); + x=contour(1).x; + y=contour(1).y; + x0=mean(x); y0=mean(y); + text(x0,y0,basin(1:end-4),'Fontsize',14); + expdisp(basin); + end +end Index: /issm/trunk-jpl/src/m/regional/showbasins.m =================================================================== --- /issm/trunk-jpl/src/m/regional/showbasins.m (revision 13010) +++ /issm/trunk-jpl/src/m/regional/showbasins.m (revision 13010) @@ -0,0 +1,69 @@ +function showbasins(varargin) +%SHOWBASINS - return basins that are within the xlim and ylim +% +% Usage: +% names=showbasins(options); +% Options: +% 'unit' default 1 +% 'hemisphere': default +1; +% 'central_meridian: 45 for Greenland and 0 for Antarctica +% 'standard_parallel: 70 for Greenland and 71 for Antarctica +% + +%is varargin an options database already? +if nargin==0, + options=pairoptions(varargin{:}); +elseif (isa(varargin{1},'plotoptions') | isa(varargin{1},'pairoptions')), + %do nothing to the varargin: + options=varargin{1}; +else + %process varargin for options: + options=pairoptions(varargin{:}); +end + + +%recover some options, and set defaults +unitmultiplier=getfieldvalue(options,'unit',1); +fontsize=getfieldvalue(options,'fontsize',12); +hemisphere=getfieldvalue(options,'hemisphere'); + +if strcmpi(hemisphere,'s'), + hemisphere=-1; +elseif strcmpi(hemisphere,'n'), + hemisphere=+1; +else + error('showbasins error message: hemispehre should be either ''n'' or ''s'''); + end + +if hemisphere==+1, + central_meridian=getfieldvalue(options,'central_meridian',45); + standard_parallel=getfieldvalue(options,'standard_parallel',70); +else + central_meridian=getfieldvalue(options,'central_meridian',0); + standard_parallel=getfieldvalue(options,'standard_parallel',71); +end + +%Ok, find basin we are talking about: +load([jplsvn '/projects/ModelData/Names/Names.mat']); + +%Get xlim and ylim, and convert into lat,long: +xlimits=xlim; x0=xlimits(1); x1=xlimits(2); +ylimits=ylim; y0=ylimits(1); y1=ylimits(2); + +%Convert names lat and long into x,y: +lat=cell2mat(names(:,3)); +long=cell2mat(names(:,2)); + +%Now, convert lat,long into x,y: +[x,y]=ll2xy(lat,long,hemisphere,central_meridian,standard_parallel); + +%Find x,y within xlimits and ylimits: +locations=find(x>x0 & xy0 & y