Index: /issm/trunk-jpl/src/m/mech/analyticaldamage.m
===================================================================
--- /issm/trunk-jpl/src/m/mech/analyticaldamage.m	(revision 16056)
+++ /issm/trunk-jpl/src/m/mech/analyticaldamage.m	(revision 16056)
@@ -0,0 +1,248 @@
+function [damage,B,backstress]=analyticaldamage(md,varargin)
+%ANALYTICALDAMAGE - compute damage for an ice shelf 
+%
+%	 This routine computes damage as a function of water/ice
+%	 material properties, ice thickness, strain rate, and ice 
+%	 rigidity.  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:
+%								'Weertman1D' for a confined ice shelf free to flow in one direction
+%								'Weertman2D' for an unconfined ice shelf free to spread in any direction
+%								'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
+%
+%   Return values:
+%		'damage' which is truncated in the range [0,1-1e-9]
+%
+%		'mask' is a mask defining areas where damage was calculated as negative prior
+%		to truncation, indicating that the ice is too warm or a back stress has not
+%		been accounted for.
+%
+%	   '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.  
+%
+%		'backstress' is the inferred backstress necessary to balance the analytical solution
+%		(keeping damage within its appropriate limits, e.g. D in [0,1]).
+%
+%   Usage:
+%      [damage,mask,B,backstress]=analyticaldamage(md,options)
+%
+%   Example:
+%      [damage,mask,B,backstress]=analyticaldamage(md,'eq','Weertman2D','smoothing',2,'backstress',10e3);
+
+% check inputs
+if (nargin<1),
+	help analyticaldamage
+	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);
+sigmab = getfieldvalue(options,'sigmab',0);
+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;
+
+%plotmodel(md,'data',a,'caxis',[-2,1])
+%plotmodel(md,'data',ex,'caxis',[-0.5e-9,0.5e-9])
+
+% spreading stress
+rhoi=md.materials.rho_ice;
+rhow=md.materials.rho_water;
+C=0.5*rhoi*md.constants.g*(1-rhoi/rhow);
+T=C*md.geometry.thickness;
+
+% rheology
+B=md.materials.rheology_B;
+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
+
+%D2=1-((theta./ex).^(1./n)).*(T-sigmab)./B;
+%length(find(D2<0))
+%length(find(ex<0))
+%length(find(D2>1))
+%length(find(theta<0))
+%pos=find(ex<0);
+%mask2=ismember(1:md.mesh.numberofvertices,pos);
+D=1-(1+a+a.^2+b.^2).^((n-1)./(2*n))./abs(ex).^(1./n).*(T-sigmab)./B./(2+a)./sign(ex);
+
+%pos=find(D2>1);
+%mask1=ismember(1:md.mesh.numberofvertices,pos);
+%pos=find(a<-2);
+%mask2=ismember(1:md.mesh.numberofvertices,pos);
+%pos=find(ex<-0);
+%mask3=ismember(1:md.mesh.numberofvertices,pos);
+%pos=find(((2+a).*sign(ex))<0);
+%mask4=ismember(1:md.mesh.numberofvertices,pos);
+%plotmodel(md,'nlines',2,'ncols',2,'data',D2,'data',a,'data',ex,'data',(2+a).*sign(ex),...
+%	'mask#1',mask1,'mask#2',mask2,'mask#3',mask3,'mask#4',mask4,...
+%	'expdisp#all','./Exp/LarsenC_DomainOutline_Bedmap2extent.exp',...
+%	'caxis#1',[0 1],'caxis#2',[-2,1],'caxis#3',[0 4e-9],'caxis#4',[0 1]);
+
+backstress=zeros(md.mesh.numberofvertices,1);
+
+% new corrections
+% D>1 where (2+a).*sign(ex)<0, compressive regions where high backstress needed
+pos=find(D>1);
+D(pos)=0;
+% backstress to bring damage to zero
+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));
+%backstress(pos)=T(pos)-(1-D(pos)).*B(pos).*(ex(pos)./theta(pos)).^(1./n(pos));
+
+pos=find(D<0);
+mask=ismember(1:md.mesh.numberofvertices,pos);
+%plotmodel(md,'data',D,'caxis',[0,1],'mask',mask)
+D(pos)=0;
+% backstress to bring negative damage to zero
+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));
+
+pos=find(backstress<0);
+%length(pos)
+backstress(pos)=0;
+
+%plotmodel(md,'data',backstress,'caxis',[0 2e5],'edgecolor','k')
+
+% 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));
+%B(pos)=T(pos)./(1-D(pos)).*(theta(pos)./ex(pos)).^(1./n(pos)); % same calculated both ways
+
+% enforce maximum B corresponding to -50 deg C ice
+%Bmax=paterson(273.15-50);
+%posB=find(B>Bmax);
+%B(posB)=Bmax;
+
+%plotmodel(md,'nlines',2,'ncols',1,'data',D,'data',D2,'caxis#all',[0,1])
+%plotmodel(md,'data',paterson(md.initialization.temperature)-B,'caxis',[-1e8 1e8])
+%plotmodel(md,'data',B,'caxis',[0.8e8 5e8])
+%plotmodel(md,'data',D2,'caxis',[0,1])
+%plotmodel(md,'data',(1-D2).*B,'caxis',[0.1e8 2e8])
+
+%plotmodel(md,'data',backstress2./T,'caxis',[1 2],'data',D2)
+%max(D2)
+%min(D2)
+
+%pos=find(D>1);
+%D(pos)=1-1e-9;
+%
+%pos=find(isnan(D));
+%D(pos)=0;
+%
+%% additional calculations here where negative damage is inferred
+%pos=find(D<0);
+%D(pos)=0;
+%
+%% magnitude of back stress to bring negative damage to zero
+%backstress=zeros(md.mesh.numberofvertices,1);
+%backstress(pos)=T(pos)-(1-D(pos)).*B(pos).*(ex(pos)./theta(pos)).^(1./n(pos));
+%
+%mask=ismember(1:md.mesh.numberofvertices,pos);
+%
+%%% enforce max back stress equal to spreading stress
+%%posS=find(backstress>T); 
+%%backstress(posS)=T(posS);
+%
+%% decreased thickness in areas of negative damage
+%dH=zeros(md.mesh.numberofvertices,1);
+%dH(pos)=1/C*((1-D(pos)).*B(pos).*(ex(pos)./theta(pos)).^(1./n(pos))-T(pos));
+%
+%% higher rigidity to eliminate areas of negative damage
+%B(pos)=T(pos)./(1-D(pos)).*(theta(pos)./ex(pos)).^(1./n(pos));
+%
+%
+%temp=[0:-0.1:-50]+273.15;
+%rig=paterson(temp);
+%newtemp=md.initialization.temperature;
+%
+%for i=1:length(pos)
+%	[mindiff,index]=min(abs(rig-B(pos(i))));
+%	newtemp(pos(i))=temp(index);
+%end
+%
+%dT=newtemp-md.initialization.temperature;
+
+%any(D<0)
+%any(D>1)
+%any(backstress<0)
+damage=D;
Index: /issm/trunk-jpl/src/m/mech/analyticaldamage.py
===================================================================
--- /issm/trunk-jpl/src/m/mech/analyticaldamage.py	(revision 16056)
+++ /issm/trunk-jpl/src/m/mech/analyticaldamage.py	(revision 16056)
@@ -0,0 +1,245 @@
+import numpy as npy
+from pairoptions import pairoptions
+from averaging import averaging
+from plotmodel import plotmodel
+
+def analyticaldamage(md,*args):
+	'''
+	ANALYTICALDAMAGE - compute damage for an ice shelf 
+	
+		 This routine computes damage as a function of water/ice
+		 material properties, ice thickness, strain rate, and ice 
+		 rigidity.  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:
+									'Weertman1D' for a confined ice shelf free to flow in one direction
+									'Weertman2D' for an unconfined ice shelf free to spread in any direction
+									'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
+	
+	   Return values:
+			'damage' which is truncated in the range [0,1-1e-9]
+	
+			'mask' is a mask defining areas where damage was calculated as negative prior
+			to truncation, indicating that the ice is too warm or a back stress has not
+			been accounted for.
+	
+		   '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.  
+	
+			'backstress' is the inferred backstress necessary to balance the analytical solution
+			(keeping damage within its appropriate limits, e.g. D in [0,1]).
+	
+	   Usage:
+	      [damage,mask,B,backstress]=analyticaldamage(md,options)
+	
+	   Example:
+	      [damage,mask,B,backstress]=analyticaldamage(md,'eq','Weertman2D','smoothing',2,'backstress',10e3)
+	'''
+
+	# check inputs
+	if 'strainrate' not in md.results.__dict__:
+		raise StandardError('md.results.strainrate not present.  Calculate using md=mechanicalproperties(md,vx,vy)')
+	if md.mesh.dimension!=2:
+		raise StandardError('only 2D model supported currently')
+	if npy.any(md.flowequation.element_equation!=2):
+		print 'Warning: the model has some non SSA elements. These will be treated like SSA elements'
+
+	# process options
+	options = pairoptions(*args)
+	eq = options.getfieldvalue('eq','Thomas')
+	smoothing = options.getfieldvalue('smoothing',0)
+	sigmab = options.getfieldvalue('sigmab',0)
+	if len(sigmab==1):
+		sigmab=sigmab*npy.ones(md.mesh.numberofvertices,)
+	
+	# 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=npy.nonzero(e1==0)
+	if npy.any(pos==1):
+		print 'WARNING: first principal strain rate equal to zero.  Value set to 1e-13 s^-1'
+		e1[pos]=1e-13
+	pos=npy.nonzero(e2==0)
+	if npy.any(pos==1):
+		disp('WARNING: second principal strain rate equal to zero.  Value set to 1e-13 s^-1');
+		e2[pos]=1e-13
+	
+	## 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=npy.arctan(md.initialization.vy/md.initialization.vx)
+	ex=0.5*(exx+eyy)+0.5*(exx-eyy)*npy.cos(2.*velangle)+exy*npy.sin(2.*velangle)
+	ey=exx+eyy-ex # trace of strain rate tensor is invariant
+	exy=-0.5*(exx-eyy)*npy.sin(2.*velangle)+exy*npy.cos(2.*velangle)
+	a=ey/ex
+	b=exy/ex
+	pos=npy.nonzero(npy.logical_or(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=npy.nonzero(abs((abs(a)-2))<1e-3)
+	a[pos]=-2+1e-3
+	
+	#plotmodel(md,'data',a,'caxis',[-2,1])
+	#plotmodel(md,'data',ex,'caxis',[-0.5e-9,0.5e-9])
+	
+	# spreading stress
+	rhoi=md.materials.rho_ice
+	rhow=md.materials.rho_water
+	C=0.5*rhoi*md.constants.g*(1.-rhoi/rhow)
+	T=C*md.geometry.thickness
+	
+	# rheology
+	B=md.materials.rheology_B
+	n=averaging(md,md.materials.rheology_n,0)
+	
+	if eq=='Weertman1D':
+		theta=1./8*npy.ones(md.mesh.numberofvertices,)
+		a=npy.zeros(md.mesh.numberofvertices,)
+	elif eq=='Weertman2D':
+		theta=1./9*npy.ones(md.mesh.numberofvertices,1)
+		a=npy.ones(md.mesh.numberofvertices,)
+	elif eq=='Thomas':
+		theta=((1+a+a**2+b**2)**((n-1)/2))/(abs(2+a)**n)
+	else:
+		raise StandardError('argument passed to "eq" not valid.  Type "help analyticaldamage" for usage')
+	
+	#D2=1-((theta./ex).^(1./n)).*(T-sigmab)./B;
+	#length(find(D2<0))
+	#length(find(ex<0))
+	#length(find(D2>1))
+	#length(find(theta<0))
+	#pos=find(ex<0);
+	#mask2=ismember(1:md.mesh.numberofvertices,pos);
+	D=1-(1+a+a**2+b**2)**((n-1)/(2*n))/abs(ex)**(1./n)*(T-sigmab)/B/(2+a)/npy.sign(ex)
+	
+	#pos=find(D2>1);
+	#mask1=ismember(1:md.mesh.numberofvertices,pos);
+	#pos=find(a<-2);
+	#mask2=ismember(1:md.mesh.numberofvertices,pos);
+	#pos=find(ex<-0);
+	#mask3=ismember(1:md.mesh.numberofvertices,pos);
+	#pos=find(((2+a).*sign(ex))<0);
+	#mask4=ismember(1:md.mesh.numberofvertices,pos);
+	#plotmodel(md,'nlines',2,'ncols',2,'data',D2,'data',a,'data',ex,'data',(2+a).*sign(ex),...
+	#	'mask#1',mask1,'mask#2',mask2,'mask#3',mask3,'mask#4',mask4,...
+	#	'expdisp#all','./Exp/LarsenC_DomainOutline_Bedmap2extent.exp',...
+	#	'caxis#1',[0 1],'caxis#2',[-2,1],'caxis#3',[0 4e-9],'caxis#4',[0 1]);
+	
+	backstress=npy.zeros(md.mesh.numberofvertices,)
+	
+	# new corrections
+	# D>1 where (2+a).*sign(ex)<0, compressive regions where high backstress needed
+	pos=npy.nonzero(D>1)
+	D[pos]=0
+	# backstress to bring damage to zero
+	backstress[pos]=T[pos]-(1-D[pos])*B[pos]*npy.sign(ex[pos])*(2+a[pos])*abs(ex[pos])**(1./n[pos])/(1+a[pos]+a[pos]**2)**((n[pos]-1)/2./n[pos])
+	#backstress[pos]=T[pos]-(1-D[pos]).*B[pos].*(ex[pos]./theta[pos]).^(1./n[pos]);
+	
+	pos=npy.nonzero(D<0)
+	#mask=ismember(1:md.mesh.numberofvertices,pos)
+	#plotmodel(md,'data',D,'caxis',[0,1],'mask',mask)
+	D[pos]=0
+	# backstress to bring negative damage to zero
+	backstress[pos]=T[pos]-(1-D[pos])*B[pos]*npy.sign(ex[pos])*(2+a[pos])*abs(ex[pos])**(1./n[pos])/(1+a[pos]+a[pos]**2)**((n[pos]-1)/2./n[pos])
+	
+	pos=npy.nonzero(backstress<0)
+	#length[pos]
+	backstress[pos]=0
+	
+	#plotmodel(md,'data',backstress,'caxis',[0 2e5],'edgecolor','k')
+	
+	# increased rigidity to bring negative damage to zero
+	B[pos]=npy.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]);
+	#B[pos]=T[pos]./(1-D[pos]).*(theta[pos]./ex[pos]).^(1./n[pos]); # same calculated both ways
+	
+	# enforce maximum B corresponding to -50 deg C ice
+	#Bmax=paterson(273.15-50);
+	#posB=find(B>Bmax);
+	#B(posB)=Bmax;
+	
+	#plotmodel(md,'nlines',2,'ncols',1,'data',D,'data',D2,'caxis#all',[0,1])
+	#plotmodel(md,'data',paterson(md.initialization.temperature)-B,'caxis',[-1e8 1e8])
+	#plotmodel(md,'data',B,'caxis',[0.8e8 5e8])
+	#plotmodel(md,'data',D2,'caxis',[0,1])
+	#plotmodel(md,'data',(1-D2).*B,'caxis',[0.1e8 2e8])
+	
+	#plotmodel(md,'data',backstress2./T,'caxis',[1 2],'data',D2)
+	#max(D2)
+	#min(D2)
+	
+	#pos=find(D>1);
+	#D[pos]=1-1e-9;
+	#
+	#pos=find(isnan(D));
+	#D[pos]=0;
+	#
+	## additional calculations here where negative damage is inferred
+	#pos=find(D<0);
+	#D[pos]=0;
+	#
+	## magnitude of back stress to bring negative damage to zero
+	#backstress=zeros(md.mesh.numberofvertices,1);
+	#backstress[pos]=T[pos]-(1-D[pos]).*B[pos].*(ex[pos]./theta[pos]).^(1./n[pos]);
+	#
+	#mask=ismember(1:md.mesh.numberofvertices,pos);
+	#
+	## enforce max back stress equal to spreading stress
+	##posS=find(backstress>T); 
+	##backstress(posS)=T(posS);
+	#
+	## decreased thickness in areas of negative damage
+	#dH=zeros(md.mesh.numberofvertices,1);
+	#dH[pos]=1/C*((1-D[pos]).*B[pos].*(ex[pos]./theta[pos]).^(1./n[pos])-T[pos]);
+	#
+	## higher rigidity to eliminate areas of negative damage
+	#B[pos]=T[pos]./(1-D[pos]).*(theta[pos]./ex[pos]).^(1./n[pos]);
+	#
+	#
+	#temp=[0:-0.1:-50]+273.15;
+	#rig=paterson(temp);
+	#newtemp=md.initialization.temperature;
+	#
+	#for i=1:length(pos)
+	#	[mindiff,index]=min(abs(rig-B(pos(i))));
+	#	newtemp(pos(i))=temp(index);
+	#end
+	#
+	#dT=newtemp-md.initialization.temperature;
+	
+	#any(D<0)
+	#any(D>1)
+	#any(backstress<0)
+	damage=D
+	
+	return damage, B, backstress