Index: /issm/trunk-jpl/src/m/mech/strainrateuncert.m
===================================================================
--- /issm/trunk-jpl/src/m/mech/strainrateuncert.m	(revision 15829)
+++ /issm/trunk-jpl/src/m/mech/strainrateuncert.m	(revision 15830)
@@ -6,6 +6,6 @@
 %   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
+%	 'dvx' and 'dvy' are velocity errors in x and y components in m/yr.  
+%   These can either be scalars or arrays of length md.mesh.numberofvertices
 %
 %   Usage:
@@ -42,4 +42,6 @@
 index=md.mesh.elements;
 summation=[1;1;1];
+dvxlist=dvx(index);
+dvylist=dvy(index);
 
 %compute nodal functions coefficients N(x,y)=alpha x + beta y +gamma
@@ -48,17 +50,17 @@
 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));
+strainrateuncert.xx=sqrt((dvxlist.*alpha).^2*summation);
+strainrateuncert.yy=sqrt((dvylist.*beta).^2*summation);
+strainrateuncert.xy=0.5*sqrt((dvxlist.*beta).^2*summation+(dvylist.*alpha).^2*summation);
 
-exx=md.results.strainrate.xx;
+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);
+p1a=strainrateuncert.xx.*(0.5+0.25*(((exx-eyy)/2).^2+exy.^2).^(-1./2).*(exx-eyy));
+p2a=strainrateuncert.yy.*(0.5-0.25*(((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*(((exx-eyy)/2).^2+exy.^2).^(-1./2).*(exx-eyy));
+p2b=strainrateuncert.yy.*(0.5+0.25*(((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);