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

Last change on this file since 26740 was 26740, checked in by Mathieu Morlighem, 3 years ago

CHG: added 25834-26739
File size: 5.1 KB
RevLine 
[26740]1Index: ../trunk-jpl/src/m/classes/materials.m
2===================================================================
3--- ../trunk-jpl/src/m/classes/materials.m (revision 25991)
4+++ ../trunk-jpl/src/m/classes/materials.m (revision 25992)
5@@ -350,6 +350,164 @@
6
7
8 end % }}}
9+ function setlitho2prem(self) % {{{
10+ %based on materials.radius and materials.issolid, replaces materials.density, materials.lame_mu, materials.lame_lambda
11+ % with volumetric averages of a polynomial version of the PREM model
12+ if ~strcmpi(self.nature, 'litho')
13+ error('materials setlitho2prem error message: materials.nature must be ''litho''')
14+ end
15+
16+
17+ ra=self.radius(end);
18+ r(1) = 0.0; r(2) = 1221.5; r(3) = 3480.0; r(4) = 3630.;
19+ r(5) = 5600.0; r(6) = 5701.0; r(7) = 5771.0; r(8) = 5971.;
20+ r(9) = 6151.0; r(10) = 6291.0; r(11) = 6346.6;
21+ r(12) = 6356.0; r(13) = 6368.0; r(14) = 6371.;
22+
23+ if (r(14)*1e3 ~= ra)
24+ 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))]);
25+ end
26+
27+ d = zeros(12,4); %polynomial coef for density
28+ d(1,1) = 13.0885; d(1,3) = -8.8381 ;
29+ d(2,1) = 12.5815; d(2,2) = -1.2638; d(2,3) = -3.6426; d(2,4) = -5.5281;
30+ d(3,1) = 7.9565 ; d(3,2) = -6.4761; d(3,3) = 5.5283; d(3,4) = -3.0807;
31+ d(4,1) = 7.9565 ; d(4,2) = -6.4761; d(4,3) = 5.5283; d(4,4) = -3.0807;
32+ d(5,1) = 7.9565 ; d(5,2) = -6.4761; d(5,3) = 5.5283; d(5,4) = -3.0807;
33+ d(6,1) = 5.3197 ; d(6,2) = -1.4836;
34+ d(7,1) = 11.2494; d(7,2) = -8.0298;
35+ d(8,1) = 7.1089 ; d(8,2) = -3.8045;
36+ d(9,1) = 2.6910 ; d(9,2) = 0.6924;
37+ d(10,1) = 2.6910; d(10,2) = 0.6924;
38+ d(11,1) = 2.9 ;
39+ d(12,1) = 2.6 ;
40+
41+ % ocean
42+ if (~self.issolid(end))
43+ d(13,1) = 1.02 ;
44+
45+ % continental
46+ else
47+ d(13,1) = d(12,1);
48+ end
49+
50+ p = zeros(13,4); %polynomial coef for Vp, pressure wave velocity
51+ p(1,1) = 11.2622 ; p(1,3) = -6.3640;
52+ p(2,1) = 11.0487 ; p(2,2) = -4.0362; p(2,3) = 4.8023; p(2,4) = -13.5732;
53+ p(3,1) = 15.3891 ; p(3,2) = -5.3181; p(3,3) = 5.5242; p(3,4) = -2.5514;
54+ p(4,1) = 24.952 ; p(4,2) = -40.4673; p(4,3) = 51.4832; p(4,4) = -26.6419;
55+ p(5,1) = 29.2766 ; p(5,2) = -23.6027; p(5,3) = 5.5242; p(5,4) = -2.5514;
56+ p(6,1) = 19.0957 ; p(6,2) = -9.8672;
57+ p(7,1) = 39.7027 ; p(7,2) = -32.6166;
58+ p(8,1) = 20.3926 ; p(8,2) = -12.2569;
59+ p(9,1) = 4.1875 ; p(9,2) = 3.9382;
60+ p(10,1) = 4.1875 ; p(10,2) = 3.9382;
61+ p(11,1) = 6.8 ;
62+ p(12,1) = 5.8;
63+ %
64+ % ocean
65+ if (~self.issolid(end))
66+ p(13,1) = 1.45 ;
67+ %
68+ % continental
69+ else
70+ p(13,1) = p(12,1);
71+ end
72+ %----
73+ %
74+ s = zeros(13,4); %polynomial coef for Vs, shear wave velocity
75+ %
76+ s(1,1) = 3.6678; s(1,3) = -4.4475;
77+
78+ s(3,1) = 6.9254; s(3,2) = 1.4672; s(3,3) = -2.0834; s(3,4) = 0.9783;
79+ s(4,1) = 11.1671; s(4,2) = -13.7818; s(4,3) = 17.4575; s(4,4) = -9.2777;
80+ s(5,1) = 22.3459; s(5,2) = -17.2473; s(5,3) = -2.0834; s(5,4) = 0.9783;
81+ s(6,1) = 9.9839; s(6,2) = -4.9324;
82+ s(7,1) = 22.3512; s(7,2) = -18.5856 ;
83+ s(8,1) = 8.9496; s(8,2) = -4.4597;
84+ s(9,1) = 2.1519; s(9,2) = 2.3481;
85+ s(10,1) = 2.1519; s(10,2) = 2.3481;
86+ s(11,1) = 3.9 ;
87+ s(12,1) = 3.2 ;
88+ %
89+ % ocean (please don't modify)
90+ if (~self.issolid(end))
91+ %
92+ % continental
93+ else
94+ s(13,1) = s(12,1);
95+ end
96+ %
97+ %
98+ r = r*1e3;
99+
100+ %- handling the first layer : central sphere
101+ rad = self.radius;
102+ rad(1) = 0.;
103+
104+ for j = 1:self.numlayers
105+
106+ ro = 0.;
107+ vp = 0.;
108+ vs = 0.;
109+
110+ for i = 1:13
111+
112+ r1 = 0.;
113+ r2 = 0.;
114+ if ((rad(j) > r(i)) & (rad(j) <= r(i+1)))
115+ if (rad(j+1) <= r(i+1))
116+ r2 = rad(j+1);
117+ r1 = rad(j);
118+ else
119+ r2 = r(i+1);
120+ r1 = rad(j);
121+ end
122+ elseif (rad(j) <= r(i))
123+ if ((rad(j+1) > r(i)) & (rad(j+1) <= r(i+1)))
124+ r2 = rad(j+1);
125+ r1 = r(i);
126+ elseif (rad(j+1) > r(i+1))
127+ r2 = r(i+1);
128+ r1 = r(i);
129+ end
130+ end
131+
132+ t1 = d(i,1)/3.;
133+ t2 = d(i,2)/(ra*4.);
134+ t3 = d(i,3)/((ra^2)*5.);
135+ t4 = d(i,4)/((ra^3)*6.);
136+ 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) );
137+
138+ t1 = p(i,1)/3.;
139+ t2 = p(i,2)/(ra*4.);
140+ t3 = p(i,3)/((ra^2)*5.);
141+ t4 = p(i,4)/((ra^3)*6.);
142+ 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) );
143+
144+ t1 = s(i,1)/3.;
145+ t2 = s(i,2)/(ra*4.);
146+ t3 = s(i,3)/((ra^2)*5.);
147+ t4 = s(i,4)/((ra^3)*6.);
148+ 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) );
149+
150+ end
151+ ro = ro*3 / (rad(j+1)^3-rad(j)^3);
152+ vp = vp*3 /(rad(j+1)^3-rad(j)^3);
153+ vs = vs*3 / (rad(j+1)^3-rad(j)^3);
154+ mu = ro*vs.^2;
155+ la = ro*vp.^2 - 2.*mu;
156+ ro = ro*1e3;
157+ la = la*1e9;
158+ mu = mu*1e9;
159+
160+ self.density(j) = ro;
161+ self.lame_lambda(j) = la;
162+ self.lame_mu(j) = mu;
163+ end
164+
165+ end % }}}
166+
167 end
168 end
169
170@@ -376,3 +534,5 @@
171 end
172 end
173 end % }}}
174+
175+
Note: See TracBrowser for help on using the repository browser.