Index: /issm/trunk-jpl/src/m/mech/mechanicalproperties.m
===================================================================
--- /issm/trunk-jpl/src/m/mech/mechanicalproperties.m	(revision 14295)
+++ /issm/trunk-jpl/src/m/mech/mechanicalproperties.m	(revision 14296)
@@ -119,10 +119,10 @@
 
 strainrate=struct('xx',[],'yy',[],'xy',[],'principalvalue1',[],'principalaxis1',[],'principalvalue2',[],'principalaxis2',[],'effectivevalue',[]);
-strainrate.xx=ux;
-strainrate.yy=vy;
-strainrate.xy=uyvx;
-strainrate.principalvalue1=valuesstrain(:,1)*(365.25*24*3600); %strain rate in 1/a instead of 1/s
+strainrate.xx=ux*md.constants.yts; %strain rate in 1/a instead of 1/s
+strainrate.yy=vy*md.constants.yts; 
+strainrate.xy=uyvx*md.constants.yts; 
+strainrate.principalvalue1=valuesstrain(:,1)*md.constants.yts; 
 strainrate.principalaxis1=directionsstrain(:,1:2);
-strainrate.principalvalue2=valuesstrain(:,2)*(365.25*24*3600); %strain rate in 1/a instead of 1/s
+strainrate.principalvalue2=valuesstrain(:,2)*md.constants.yts; 
 strainrate.principalaxis2=directionsstrain(:,3:4);
 strainrate.effectivevalue=1/sqrt(2)*sqrt(strainrate.xx.^2+strainrate.yy.^2+2*strainrate.xy.^2);
Index: /issm/trunk-jpl/src/m/mech/strainrateuncert.m
===================================================================
--- /issm/trunk-jpl/src/m/mech/strainrateuncert.m	(revision 14296)
+++ /issm/trunk-jpl/src/m/mech/strainrateuncert.m	(revision 14296)
@@ -0,0 +1,71 @@
+function md=strainratuncert(md,vx,vy,dvx,dvy)
+%STRAINRATEUNCERT - compute uncertainty in strain rate components
+%
+%   this routine computes the uncertainties in the strain rate tensor
+%	 components given the uncertainty in surface velocity data.
+%   The results are stored in md.results
+%
+%	 'dvx' and 'dvy' are velocity errors in x and y components.  These 
+%	 can either be scalars or arrays of length md.mesh.numberofvertices
+%
+%   Usage:
+%      md=strainrateuncert(md,vx,vy,dv)
+%
+%   Example:
+%      md=mechanicalproperties(md,md.initialization.vx,md.initialization.vy,5);
+%      md=mechanicalproperties(md,md.inversion.vx_obs,md.inversion.vy_obs,dv);
+
+%some checks
+if length(vx)~=md.mesh.numberofvertices | length(vy)~=md.mesh.numberofvertices,
+	error(['the input velocity should be of size ' num2str(md.mesh.numberofvertices) '!'])
+end
+if length(dvx)==1,
+	dvx=dvx*ones(md.mesh.numberofelements,1);
+end
+if length(dvx)~=md.mesh.numberofelements,
+	error(['the velocity error dvx should be of size ' num2str(md.mesh.numberofelements) ' or 1!'])
+end
+if length(dvy)==1,
+	dvy=dvy*ones(md.mesh.numberofelements,1);
+end
+if length(dvy)~=md.mesh.numberofelements,
+	error(['the velocity error dvy should be of size ' num2str(md.mesh.numberofelements) ' or 1!'])
+end
+if ~(md.mesh.dimension==2)
+	error('only 2d model supported yet');
+end
+if any(md.flowequation.element_equation~=2),
+	disp('Warning: the model has some non macayeal elements. These will be treated like MacAyeal''s elements');
+end
+
+%initialization
+index=md.mesh.elements;
+summation=[1;1;1];
+
+%compute nodal functions coefficients N(x,y)=alpha x + beta y +gamma
+[alpha beta]=GetNodalFunctionsCoeff(index,md.mesh.x,md.mesh.y);
+
+strainrateuncert=struct('xx',[],'yy',[],'xy',[],'principalvalue1',[],'principalvalue2',[],'effectivevalue',[]);
+
+strainrateuncert.xx=dvx.*sqrt(alpha.^2*summation);
+strainrateuncert.yy=dvy.*sqrt(beta.^2*summation);
+strainrateuncert.xy=0.5*sqrt(dvx.^2.*(beta.^2*summation)+dvy.^2.*(alpha.^2*summation));
+
+exx=md.results.strainrate.xx;
+eyy=md.results.strainrate.yy;
+exy=md.results.strainrate.xy;
+p1a=strainrateuncert.xx.*(0.5+0.25*(0.5*((exx-eyy)/2).^2+exy.^2).^(-1./2).*(exx-eyy));
+p2a=strainrateuncert.yy.*(0.5-0.25*(0.5*((exx-eyy)/2).^2+exy.^2).^(-1./2).*(exx-eyy));
+p3a=strainrateuncert.xy.*(((exx-eyy)/2).^(2)+exy.^2).^(-1./2).*exy;
+p1b=strainrateuncert.xx.*(0.5-0.25*(0.5*((exx-eyy)/2).^2+exy.^2).^(-1./2).*(exx-eyy));
+p2b=strainrateuncert.yy.*(0.5+0.25*(0.5*((exx-eyy)/2).^2+exy.^2).^(-1./2).*(exx-eyy));
+p3b=strainrateuncert.xy.*(-(((exx-eyy)/2).^(2)+exy.^2).^(-1./2).*exy);
+strainrateuncert.principalvalue1=sqrt(p1a.^2+p2a.^2+p3a.^2);
+strainrateuncert.principalvalue2=sqrt(p1b.^2+p2b.^2+p3b.^2);
+
+effa=strainrateuncert.xx/sqrt(2).*(exx.^2+eyy.^2+2*exy.^2).^(-1./2).*exx;
+effb=strainrateuncert.yy/sqrt(2).*(exx.^2+eyy.^2+2*exy.^2).^(-1./2).*eyy;
+effc=2*strainrateuncert.xy/sqrt(2).*(exx.^2+eyy.^2+2*exy.^2).^(-1./2).*exy;
+strainrateuncert.effectivevalue=sqrt(effa.^2+effb.^2+effc.^2);
+
+md.results.strainrateuncert=strainrateuncert;