Index: /issm/trunk-jpl/src/m/mech/analyticaldamage.m
===================================================================
--- /issm/trunk-jpl/src/m/mech/analyticaldamage.m	(revision 16415)
+++ /issm/trunk-jpl/src/m/mech/analyticaldamage.m	(revision 16416)
@@ -13,17 +13,24 @@
 %								'Thomas' for a 2D ice shelf, taking into account full strain rate tensor (default)
 %		- 'smoothing'	: the amount of smoothing to be applied to the strain rate data.
-%								Type 'help averaging' for more information on its usage.
-%		- 'sigmab'		: a compressive backstress term to be subtracted from the driving stress 
-%								in the damage calculation
+%								Type 'help averaging' for more information on its
+%								usage. Defaults to 0.
+%		- 'sigmab'		: a-priori backstress in opposition to the driving stress.
+%								Defaults to 0 everywhere. 
+
+%		- 'coordsys'	: coordinate system for calculating the strain rate
+%							components. Must be one of: 
+%				'longitudinal': x axis aligned along a flowline at every point (default)
+%				'principal': x axis aligned along maximum principal strain rate
+%					at every point
+%				'xy': x and y axes same as in polar stereographic projection 
 %
 %   Return values:
 %		'damage' which is truncated in the range [0,1-1e-9]
 %
-%	   'B' is the rigidity, which is equal to md.materials.rheology_B in areas outside
-%		those defined by 'mask.'  Within areas defined by 'mask,' where negative damage 
-%		is inferred, 'B' is updated to make damage equal to zero.  
+%	   'B' is the ice rigidity calculated assuming D=0 everywhere. 
 %
-%		'backstress' is the inferred backstress necessary to balance the analytical solution
-%		(keeping damage within its appropriate limits, e.g. D in [0,1]).
+%		'backstress' is the inferred backstress necessary to balance the
+%		analytical solution (keeping damage within its appropriate limits, e.g. D
+%		in [0,1]).
 %
 %   Usage:
@@ -31,5 +38,5 @@
 %
 %   Example:
-%      [damage,B,backstress]=analyticaldamage(md,'eq','Weertman2D','smoothing',2,'backstress',10e3);
+%      [damage,B,backstress]=analyticaldamage(md,'eq','Weertman2D','smoothing',2,'sigmab',10e3,'coordsys','longitudinal');
 
 % check inputs
@@ -53,60 +60,10 @@
 smoothing = getfieldvalue(options,'smoothing',0);
 sigmab = getfieldvalue(options,'sigmab',0);
+coordsys = getfieldvalue(options,'coordsys','longitudinal');
 if length(sigmab)==1,
 	sigmab=sigmab*ones(md.mesh.numberofelements,1);
 end
 
-% average element strain rates onto vertices
-e1=averaging(md,md.results.strainrate.principalvalue1,smoothing)/md.constants.yts; % convert to s^-1
-e2=averaging(md,md.results.strainrate.principalvalue2,smoothing)/md.constants.yts;
-exx=averaging(md,md.results.strainrate.xx,smoothing)/md.constants.yts;
-eyy=averaging(md,md.results.strainrate.yy,smoothing)/md.constants.yts;
-exy=averaging(md,md.results.strainrate.xy,smoothing)/md.constants.yts;
-
-% checks: any of e1 or e2 equal to zero?
-pos=find(e1==0);
-if any(pos==1)
-	disp('WARNING: first principal strain rate equal to zero.  Value set to 1e-13 s^-1');
-	e1(pos)=1e-13;
-end
-pos=find(e2==0);
-if any(pos==1)
-	disp('WARNING: second principal strain rate equal to zero.  Value set to 1e-13 s^-1');
-	e2(pos)=1e-13;
-end
-
-%% old method using principal strain rates {{{
-%% ex=maximum principal tensile strain rate
-%ex=e1;
-%a=e2./e1;
-%pos=find(e1<0 & e2>0); % longitudinal compression and lateral tension
-%a(pos)=e1(pos)./e2(pos);
-%ex(pos)=e2(pos);
-%pos2=find(e1<0 & e2<0 & abs(e1)<abs(e2)); % lateral and longitudinal compression
-%a(pos2)=e1(pos2)./e2(pos2);
-%ex(pos2)=e2(pos2);
-%pos3=find(e1>0 & e2>0 & abs(e1)<abs(e2)); % lateral and longitudinal tension 
-%a(pos3)=e1(pos3)./e2(pos3);
-%ex(pos3)=e2(pos3);
-%id=find(e1<0 & e2<0);
-%a(id)=-a(id); % where both strain rates are compressive, enforce negative alpha
-%
-%% }}}
-
-% new method using longitudinal strain rates defined by observed velocity vector
-velangle=atan(md.initialization.vy./md.initialization.vx);
-ex=0.5*(exx+eyy)+0.5*(exx-eyy).*cos(2*velangle)+exy.*sin(2*velangle);
-ey=exx+eyy-ex; % trace of strain rate tensor is invariant
-exy=-0.5*(exx-eyy).*sin(2*velangle)+exy.*cos(2*velangle);
-a=ey./ex;
-b=exy./ex;
-pos=find(ex<0 & ey<0);
-%length(pos)
-a(pos)=-a(pos);
-
-% a < -1 in areas of strong lateral compression or longitudinal compression
-% and theta is undefined at a = -2
-pos=find(abs((abs(a)-2))<1e-3);
-a(pos)=-2+1e-3;
+[a,b,theta,ex]=thomasparams(md,'eq',eq,'smoothing',smoothing,'coordsys',coordsys);
 
 % spreading stress
@@ -120,20 +77,5 @@
 n=averaging(md,md.materials.rheology_n,0);
 
-switch eq
-	case 'Weertman1D'
-		theta=1./8*ones(md.mesh.numberofvertices,1);
-		a=zeros(md.mesh.numberofvertices,1);
-	case 'Weertman2D'
-		theta=1./9*ones(md.mesh.numberofvertices,1);
-		a=ones(md.mesh.numberofvertices,1);
-	case 'Thomas'
-		theta=((1+a+a.^2+b.^2).^((n-1)/2))./(abs(2+a).^n);
-	otherwise
-		error('argument passed to "eq" not valid.  Type "help analyticaldamage" for usage');
-end
-
 D=1-(1+a+a.^2+b.^2).^((n-1)./(2*n))./abs(ex).^(1./n).*(T-sigmab)./B./(2+a)./sign(ex);
-
-backstress=zeros(md.mesh.numberofvertices,1);
 
 % D>1 where (2+a).*sign(ex)<0, compressive regions where high backstress needed
@@ -141,5 +83,7 @@
 D(pos)=0;
 
-% backstress to bring damage to zero
+backstress=zeros(md.mesh.numberofvertices,1);
+
+% backstress to bring D down to one 
 backstress(pos)=T(pos)-(1-D(pos)).*B(pos).*sign(ex(pos)).*(2+a(pos)).*abs(ex(pos)).^(1./n(pos))./(1+a(pos)+a(pos).^2).^((n(pos)-1)/2./n(pos));
 
@@ -154,6 +98,8 @@
 backstress(pos)=0;
 
-% increased rigidity to bring negative damage to zero
-B(pos)=sign(ex(pos))./(2+a(pos)).*(1+a(pos)+a(pos).^2).^((n(pos)-1)/2./n(pos)).*T(pos)./abs(ex(pos)).^(1./n(pos));
+% rigidity from Thomas relation for D=0 and backstress=0
+B=sign(ex)./(2+a).*(1+a+a.^2).^((n-1)/2./n).*T./(abs(ex).^(1./n));
+pos=find(B<0);
+B(pos)=md.materials.rheology_B(pos);
 
 damage=D;
Index: /issm/trunk-jpl/src/m/mech/backstressfrominversion.m
===================================================================
--- /issm/trunk-jpl/src/m/mech/backstressfrominversion.m	(revision 16416)
+++ /issm/trunk-jpl/src/m/mech/backstressfrominversion.m	(revision 16416)
@@ -0,0 +1,71 @@
+function backstress=backstressfrominversion(md,varargin)
+%BACKSTRESSFROMINVERSION - compute ice shelf backstress from inversion results 
+%
+%	 This routine computes backstress based on the analytical formalism of
+%	 Thomas (1973) and Borstad et al. (2013).  The model must contain inversion
+%	 results for ice rigidity.  Strain rates must also be included, either from
+%	 observed or modeled velocities.  Ice rigidity B is assumed to be
+%	 parameterized by the ice temperature in md.materials.rheology_B.
+%
+%   Available options:
+%		- 'tempmask'	: mask the inverted rigidity to be no more than
+%							appropriate for the temperature of the ice?  
+%							Boolean, defaults to false.
+%		- 'smoothing'	: the amount of smoothing to be applied to the strain rate data.
+%								Type 'help averaging' for more information on its
+%								usage. Defaults to 0.
+%		- 'coordsys'	: coordinate system for calculating the strain rate
+%							components. Must be one of: 
+%				'longitudinal': x axis aligned along a flowline at every point (default)
+%				'principal': x axis aligned along maximum principal strain rate
+%					at every point
+%				'xy': x and y axes same as in polar stereographic projection 
+%
+%   Return values:
+%		'backstress' is the inferred backstress necessary to balance the
+%		analytical solution (keeping damage within its appropriate limits, e.g. D
+%		in [0,1]).
+%
+%   Usage:
+%      backstress=backstressfrominversion(md,options)
+%
+%   Example:
+%      backstress=backstressfrominversion(md,'smoothing',2,'coordsys','longitudinal','tempmask',true);
+
+% check inputs
+if (nargin<1),
+	help backstressfrominversion
+	error('bad usage');
+end
+if isempty(fieldnames(md.results)),
+	error(['md.results.strainrate is not present.  Calculate using md=mechanicalproperties(md,vx,vy)']);
+end
+if ~(md.mesh.dimension==2),
+	error('only 2d model supported currently');
+end
+if any(md.flowequation.element_equation~=2),
+	disp('Warning: the model has some non SSA elements. These will be treated like SSA elements');
+end
+
+% process options
+options = pairoptions(varargin{:});
+smoothing = getfieldvalue(options,'smoothing',0);
+coordsys = getfieldvalue(options,'coordsys','longitudinal');
+tempmask = getfieldvalue(options,'tempmask',false);
+
+T=0.5*md.materials.rho_ice*md.constants.g*(1-md.materials.rho_ice/md.materials.rho_water)*md.geometry.thickness;
+n=averaging(md,md.materials.rheology_n,0);
+B=md.materials.rheology_B;
+Bi=md.results.StressbalanceSolution.MaterialsRheologyBbar;
+
+[a0,b0,theta0,ex0]=thomasparams(md,'eq','Thomas','smoothing',smoothing,'coordsys',coordsys);
+
+if tempmask
+	Bi=md.results.StressbalanceSolution.MaterialsRheologyBbar;
+   pos=find(Bi>md.materials.rheology_B);
+   Bi(pos)=md.materials.rheology_B(pos);
+end
+
+% analytical backstress solution
+backstress=T-Bi.*sign(ex0).*(2+a0).*abs(ex0).^(1./n)./((1+a0+a0.^2+b0.^2).^((n-1)/2./n));
+backstress(find(backstress<0))=0;
Index: /issm/trunk-jpl/src/m/mech/damagefrominversion.m
===================================================================
--- /issm/trunk-jpl/src/m/mech/damagefrominversion.m	(revision 16416)
+++ /issm/trunk-jpl/src/m/mech/damagefrominversion.m	(revision 16416)
@@ -0,0 +1,33 @@
+function damage=damagefrominversion(md)
+%DAMAGEFROMINVERSION - compute ice shelf damage from inversion results
+%
+%	This routine computes damage based on the analytical formalism of Borstad et
+%	al. (2013).  The model must contain inversion results for ice rigidity.  Ice
+%	rigidity B is assumed to be parameterized by the ice temperature in
+%	md.materials.rheology_B. 
+%
+%	Usage:
+%		damage=damagefrominversion(md)
+%
+%	Example:
+%		damage=damagefrominversion(md)
+
+% check inputs
+if (nargin<1),
+	help backstressfrominversion
+	error('bad usage');
+end
+if isempty(fieldnames(md.results)),
+	error(['md.results.strainrate is not present.  Calculate using md=mechanicalproperties(md,vx,vy)']);
+end
+if ~(md.mesh.dimension==2),
+	error('only 2d model supported currently');
+end
+if any(md.flowequation.element_equation~=2),
+	disp('Warning: the model has some non SSA elements. These will be treated like SSA elements');
+end
+
+damage=zeros(md.mesh.numberofvertices,1);
+% Damage where Bi softer than B(T)
+pos=find(md.results.StressbalanceSolution.MaterialsRheologyBbar<md.materials.rheology_B);
+damage(pos)=1-md.results.StressbalanceSolution.MaterialsRheologyBbar(pos)./md.materials.rheology_B(pos);
Index: /issm/trunk-jpl/src/m/mech/thomasparams.m
===================================================================
--- /issm/trunk-jpl/src/m/mech/thomasparams.m	(revision 16416)
+++ /issm/trunk-jpl/src/m/mech/thomasparams.m	(revision 16416)
@@ -0,0 +1,147 @@
+function [alpha,beta,theta,ex]=thomasparams(md,varargin)
+%THOMASPARAMS - compute Thomas' geometric parameters for an ice shelf 
+%
+%	 This routine computes geometric parameters representing ratios between
+%	 components of the horizontal strain rate tensor for an ice shelf, as
+%	 originally developed in Thomas (1973).  The model must contain computed
+%	 strain rates, either from observed or modeled ice velocities.
+%
+%   Available options:
+%	 -'eq'			: analytical equation to use in the calculation.  Must be one of:
+%				'Thomas' for a 2D ice shelf, taking into account full strain rate
+%					tensor (default)
+%				'Weertman1D' for a confined ice shelf free to flow in one direction
+%				'Weertman2D' for an unconfined ice shelf free to spread in any direction
+%
+%	 -'smoothing'	: an integer smoothing parameter for the averaging function
+%						(default 0) Type 'help averaging' for more information on its usage.
+%
+%	 -'coordsys'	: coordinate system for calculating the strain rate
+%						components. Must be one of:
+%				'longitudinal': x axis aligned along a flowline at every point (default)
+%				'principal': x axis aligned along maximum principal strain rate
+%					at every point
+%				'xy': x and y axes same as in polar stereographic projection 
+%
+%   Return values: 
+%
+%		'alpha' which is the ratio e_yy/e_xx between components of the strain
+%		rate tensor
+%
+%		'beta' which is the ratio e_xy/e_xx between components of the strain rate
+%		tensor
+%
+%		'theta' which is a combination of alpha and beta arising from the form of
+%		the equivalent stress
+%
+%		'exx' is the longitudinal strain rate along a coordinate system defined
+%		by 'coordsys' 
+%
+%   Usage: [alpha,beta,theta,exx]=ThomasParams(md,options)
+%
+%   Example: [alpha,beta,theta,exx]=ThomasParams(md,'eq','Thomas','smoothing',2,'coordsys','longitudinal')
+
+%some checks
+if (nargin<4)
+	help ThomasParams
+	error('bad usage');
+end
+if isempty(fieldnames(md.results)),
+	error(['md.results.strainrate is not present.  Calculate using md=mechanicalproperties(md,vx,vy)'])
+end
+if ~(md.mesh.dimension==2)
+	error('only 2d model supported currently');
+end
+if any(md.flowequation.element_equation~=2),
+	disp('Warning: the model has some non SSA elements. These will be treated like SSA elements');
+end
+
+% process options
+options = pairoptions(varargin{:});
+eq = getfieldvalue(options,'eq','Thomas');
+smoothing = getfieldvalue(options,'smoothing',0);
+coordsys = getfieldvalue(options,'coordsys','longitudinal');
+
+% average element strain rates onto vertices
+e1=averaging(md,md.results.strainrate.principalvalue1,smoothing)/md.constants.yts; % convert to s^-1
+e2=averaging(md,md.results.strainrate.principalvalue2,smoothing)/md.constants.yts;
+exx=averaging(md,md.results.strainrate.xx,smoothing)/md.constants.yts;
+eyy=averaging(md,md.results.strainrate.yy,smoothing)/md.constants.yts;
+exy=averaging(md,md.results.strainrate.xy,smoothing)/md.constants.yts;
+
+% checks: any of e1 or e2 equal to zero?
+pos=find(e1==0);
+if any(pos==1)
+	disp('WARNING: first principal strain rate equal to zero.  Value set to 1e-13 s^-1');
+	e1(pos)=1e-13;
+end
+pos=find(e2==0);
+if any(pos==1)
+	disp('WARNING: second principal strain rate equal to zero.  Value set to 1e-13 s^-1');
+	e2(pos)=1e-13;
+end
+
+switch coordsys
+	case 'principal'
+		b=zeros(md.mesh.numberofvertices,1);
+		ex=e1;
+		a=e2./e1;
+		pos=find(e1<0 & e2>0); % longitudinal compression and lateral tension
+		a(pos)=e1(pos)./e2(pos);
+		ex(pos)=e2(pos);
+		pos2=find(e1<0 & e2<0 & abs(e1)<abs(e2)); % lateral and longitudinal compression
+		a(pos2)=e1(pos2)./e2(pos2);
+		ex(pos2)=e2(pos2);
+		pos3=find(e1>0 & e2>0 & abs(e1)<abs(e2)); % lateral and longitudinal tension
+		a(pos3)=e1(pos3)./e2(pos3);
+		ex(pos3)=e2(pos3);
+		id=find(e1<0 & e2<0);
+		a(id)=-a(id); % where both strain rates are compressive, enforce negative alpha
+		%mask=ismember(1:md.mesh.numberofvertices,id);
+		%plotmodel(md,'data',ex,'mask',mask)
+	case 'xy'
+		ex=exx;
+		a=eyy./exx;
+		b=exy./exx;
+	case 'longitudinal'
+		% using longitudinal strain rates defined by observed velocity vector
+		velangle=atan(md.initialization.vy./md.initialization.vx);
+		pos=find(md.initialization.vx==0);
+		velangle(pos)=pi/2;
+		ex=0.5*(exx+eyy)+0.5*(exx-eyy).*cos(2*velangle)+exy.*sin(2*velangle);
+		ey=exx+eyy-ex; % trace of strain rate tensor is invariant
+		exy=-0.5*(exx-eyy).*sin(2*velangle)+exy.*cos(2*velangle);
+		a=ey./ex;
+		b=exy./ex;
+		pos=find(ex<0 & ey<0);
+		%length(pos)
+		a(pos)=-a(pos);
+	otherwise
+		error('argument passed to "coordsys" not valid');
+end
+
+% a < -1 in areas of strong lateral compression or longitudinal compression and
+% theta flips sign at a = -2
+pos=find(abs((abs(a)-2))<1e-3);
+if length(pos)>0,
+	disp(['Warning: ', num2str(length(pos)), ' vertices have alpha within 1e-3 of -2'])
+end
+a(pos)=-2+1e-3;
+
+% rheology
+n=averaging(md,md.materials.rheology_n,0);
+
+switch eq
+	case 'Weertman1D'
+		theta=1./8;
+		a=zeros(md.mesh.numberofvertices,1);
+	case 'Weertman2D'
+		theta=1./9;
+		a=ones(md.mesh.numberofvertices,1);
+	case 'Thomas'
+		theta=((1+a+a.^2+b.^2).^((n-1)/2))./(abs(2+a).^n);
+	otherwise
+		error('argument passed to "eq" not valid.  Type "help zinv" for usage');
+end
+alpha=a;
+beta=b;