Index: /issm/branches/trunk-jpl-damage/src/m/model/initialdamage.m
===================================================================
--- /issm/branches/trunk-jpl-damage/src/m/model/initialdamage.m	(revision 12919)
+++ /issm/branches/trunk-jpl-damage/src/m/model/initialdamage.m	(revision 12919)
@@ -0,0 +1,41 @@
+function damage=initialdamage(md)
+%INITIALDAMAGE - compute initial damage guess for ICE SHELF inverse control method 
+%
+%	 This routine computes the initial damage guess as a function of
+%	 material properties, ice thickness, strain rate, and ice 
+%	 rigidity.  The model must contain computed strain rates.
+%
+%   Usage:
+%      md=initialdamage(md)
+%
+%   Example:
+%      md=initialdamage(md);
+
+%some checks
+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 macayeal elements. These will be treated like MacAyeal''s elements');
+end
+
+%average results onto vertices
+eps1=averaging(md,md.results.strainrate.principalvalue1,0)/md.constants.yts; % s^-1 for strain rates
+epsxx=averaging(md,md.results.strainrate.xx,0)/md.constants.yts;
+epsyy=averaging(md,md.results.strainrate.yy,0)/md.constants.yts;
+epseff=averaging(md,md.results.strainrate.effectivevalue,0)/md.constants.yts;
+n=averaging(md,md.materials.rheology_n,0);
+
+%lump material constants together
+const=md.materials.rho_ice*md.constants.g*(1-md.materials.rho_ice/md.materials.rho_water)/2;
+
+damage=1-0.5*const.*md.geometry.thickness./md.materials.rheology_B./eps1.^(1./n);
+%damage=1-const.*md.geometry.thickness.*epseff.^((n-1)./n)./md.materials.rheology_B./(eps1+epsxx+epsyy);
+
+pos=find(damage>1);
+damage(pos)=1;
+pos=find(damage<0);
+damage(pos)=0;