source: issm/oecreview/Archive/25834-26739/ISSM-25991-25992.diff@ 28275

../trunk-jpl/src/m/classes/materials.m

@@ -350 +350 @@
 
 
                 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
 
@@ -376 +534 @@
                 end
         end
 end % }}}
+
+