Index: /issm/trunk-jpl/src/m/classes/materials.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/materials.m	(revision 25991)
+++ /issm/trunk-jpl/src/m/classes/materials.m	(revision 25992)
@@ -351,4 +351,162 @@
 
 		end % }}}
+		function setlitho2prem(self) % {{{
+			%based on materials.radius and materials.issolid, replaces materials.density, materials.lame_mu, materials.lame_lambda
+			% with volumetric averages of a polynomial version of the PREM model
+			if ~strcmpi(self.nature, 'litho')
+			error('materials setlitho2prem error message: materials.nature must be ''litho''')
+			end
+		   
+		   
+		   ra=self.radius(end);
+		   r(1) = 0.0;     r(2) = 1221.5; r(3) = 3480.0; r(4) = 3630.;
+		   r(5) = 5600.0;  r(6) = 5701.0;  r(7) = 5771.0; r(8) = 5971.;
+		   r(9) = 6151.0; r(10) = 6291.0; r(11) = 6346.6;
+		   r(12) = 6356.0; r(13) = 6368.0; r(14) = 6371.;
+		  
+		   if (r(14)*1e3 ~= ra) 
+			error(['materials setlitho2prem error message: Earth surface radius in PREM r=' num2str(r(end)*1e3) ' does not match materials.radius(end)=' num2str(self.radius(end))]); 
+		   end
+
+		  d = zeros(12,4); %polynomial coef for density
+		  d(1,1) = 13.0885; 		         	  d(1,3) = -8.8381 ;
+		  d(2,1) = 12.5815; d(2,2) = -1.2638; d(2,3) = -3.6426; d(2,4) = -5.5281;
+		  d(3,1) = 7.9565 ; d(3,2) = -6.4761;   d(3,3) = 5.5283;  d(3,4) = -3.0807;
+		  d(4,1) = 7.9565 ; d(4,2) = -6.4761;   d(4,3) = 5.5283;  d(4,4) = -3.0807;
+		  d(5,1) = 7.9565 ; d(5,2) = -6.4761;   d(5,3) = 5.5283;  d(5,4) = -3.0807;
+		  d(6,1) = 5.3197 ; d(6,2) = -1.4836;
+		  d(7,1) = 11.2494; d(7,2) = -8.0298;
+		  d(8,1) = 7.1089 ; d(8,2) = -3.8045;
+		  d(9,1) = 2.6910 ; d(9,2) = 0.6924;
+		  d(10,1) = 2.6910; d(10,2) = 0.6924;
+		  d(11,1) = 2.9  ;
+		  d(12,1) = 2.6  ;
+
+		% ocean
+		  if (~self.issolid(end)) 
+		  	d(13,1) = 1.02 ;
+
+		% continental
+		  else
+		  d(13,1) = d(12,1);
+		  end
+
+		  p = zeros(13,4); %polynomial coef for Vp, pressure wave velocity
+		  p(1,1) = 11.2622 ; p(1,3) = -6.3640;
+		  p(2,1) = 11.0487 ; p(2,2) = -4.0362; p(2,3)  = 4.8023; p(2,4) = -13.5732;
+		  p(3,1) = 15.3891 ; p(3,2) = -5.3181; p(3,3)  = 5.5242; p(3,4) = -2.5514;
+		  p(4,1) = 24.952 ; p(4,2)  = -40.4673; p(4,3) = 51.4832; p(4,4) = -26.6419;
+		  p(5,1) = 29.2766 ; p(5,2) = -23.6027; p(5,3) = 5.5242; p(5,4) = -2.5514;
+		  p(6,1) = 19.0957 ; p(6,2)  = -9.8672;
+		  p(7,1) = 39.7027 ; p(7,2)  = -32.6166;
+		  p(8,1) = 20.3926 ; p(8,2)  = -12.2569;
+		  p(9,1) = 4.1875 ; p(9,2)  = 3.9382;
+		  p(10,1) = 4.1875 ; p(10,2) = 3.9382;
+		  p(11,1) = 6.8 ;
+		  p(12,1) = 5.8;
+		%
+		% ocean
+		  if (~self.issolid(end)) 
+		  p(13,1) = 1.45 ;
+		%
+		% continental
+		  else
+		  p(13,1) = p(12,1);
+		  end
+		%----
+		%
+		  s = zeros(13,4); %polynomial coef for Vs, shear wave velocity
+		%
+		  s(1,1) = 3.6678; s(1,3) = -4.4475;
+
+		  s(3,1) = 6.9254; s(3,2) = 1.4672; s(3,3) = -2.0834; s(3,4) = 0.9783;
+		  s(4,1) = 11.1671; s(4,2) = -13.7818; s(4,3) = 17.4575; s(4,4) = -9.2777;
+		  s(5,1) = 22.3459; s(5,2) = -17.2473; s(5,3) = -2.0834; s(5,4) = 0.9783;
+		  s(6,1) = 9.9839; s(6,2) = -4.9324;
+		  s(7,1) = 22.3512; s(7,2) = -18.5856 ;
+		  s(8,1) = 8.9496; s(8,2) = -4.4597;
+		  s(9,1) = 2.1519; s(9,2) = 2.3481;
+		  s(10,1) = 2.1519; s(10,2) = 2.3481;
+		  s(11,1) = 3.9 ;
+		  s(12,1) = 3.2 ;
+		%
+		% ocean (please don't modify)
+		  if (~self.issolid(end))
+		%
+		% continental
+		  else
+		  s(13,1) = s(12,1);
+		  end
+		%
+		%
+		  r = r*1e3;
+		  
+		  %- handling the first layer : central sphere
+		  rad = self.radius;
+		  rad(1) = 0.;
+		  
+		  for j = 1:self.numlayers
+		  
+			ro = 0.;
+			vp = 0.;
+			vs = 0.;
+
+			for i = 1:13
+
+				r1 = 0.;
+				r2 = 0.;
+				if ((rad(j) > r(i)) & (rad(j) <= r(i+1))) 
+					if (rad(j+1) <= r(i+1)) 
+						r2 = rad(j+1);
+						r1 = rad(j);
+					else
+						r2 = r(i+1);
+						r1 = rad(j);
+					end
+				elseif (rad(j) <= r(i))
+					if ((rad(j+1) > r(i)) & (rad(j+1) <= r(i+1)))
+						r2 = rad(j+1);
+						r1 = r(i);
+					elseif (rad(j+1) > r(i+1))
+						r2 = r(i+1);
+						r1 = r(i);
+					end
+				end
+
+				t1 = d(i,1)/3.;
+				t2 = d(i,2)/(ra*4.);
+				t3 = d(i,3)/((ra^2)*5.);
+				t4 = d(i,4)/((ra^3)*6.);
+				ro =  ro + t1*(r2^3) + t2*(r2^4) + t3*(r2^5) + t4*(r2^6) - ( t1*(r1^3) + t2*(r1^4) + t3*(r1^5) + t4*(r1^6) );
+					  
+				t1 = p(i,1)/3.;
+				t2 = p(i,2)/(ra*4.);
+				t3 = p(i,3)/((ra^2)*5.);
+				t4 = p(i,4)/((ra^3)*6.);
+				vp =  vp + t1*(r2^3) + t2*(r2^4) + t3*(r2^5) + t4*(r2^6) - ( t1*(r1^3) + t2*(r1^4) + t3*(r1^5) + t4*(r1^6) );
+					  
+				t1 = s(i,1)/3.;
+				t2 = s(i,2)/(ra*4.);
+				t3 = s(i,3)/((ra^2)*5.);
+				t4 = s(i,4)/((ra^3)*6.);
+				vs =  vs + t1*(r2^3) + t2*(r2^4) + t3*(r2^5) + t4*(r2^6) - ( t1*(r1^3) + t2*(r1^4) + t3*(r1^5) + t4*(r1^6) );
+
+			end 
+			ro = ro*3 / (rad(j+1)^3-rad(j)^3);
+			vp = vp*3 /(rad(j+1)^3-rad(j)^3);
+			vs = vs*3 / (rad(j+1)^3-rad(j)^3);
+			mu = ro*vs.^2;
+			la = ro*vp.^2 - 2.*mu;
+			ro = ro*1e3;
+			la = la*1e9;
+			mu = mu*1e9;
+
+			self.density(j) = ro;
+			self.lame_lambda(j) = la;
+			self.lame_mu(j) = mu;
+		  end
+
+		end % }}}
+   
 	end
 end
@@ -377,2 +535,4 @@
 	end
 end % }}}
+
+