|
1 function [dev,d,h,i,f,x,y,z]=magdev(flat,flon,elevkm,nmax,igrf); |
|
2 % function [dev,d,h,i,f,x,y,z]=magdev(flat,flon,elevkm,nmax,igrf); |
|
3 % |
|
4 % compute magnetic deviation |
|
5 % input: |
|
6 % flat latitude (degree) |
|
7 % flon longitude (degree) |
|
8 % elevkm elevation above mean geoid (km) |
|
9 % nmax number of harmonics |
|
10 % igrf (need file, default 'IGRF95') |
|
11 % |
|
12 % output: |
|
13 % dev mag. deviation (degree) |
|
14 % |
|
15 % |
|
16 % based of FORTRAN ROUTINE GEOMAG.FOR |
|
17 % more info under http://fdd.gsfc.nasa.gov/IGRF.html |
|
18 % |
|
19 % M. Visbeck, LDEO FEB 2000 |
|
20 % |
|
21 |
|
22 if nargin<5 |
|
23 igrf='IGRF95'; |
|
24 igrf='IGRF00'; |
|
25 end |
|
26 eval(['gh=',igrf,';']) |
|
27 |
|
28 if nargin<4 |
|
29 nmax=10; |
|
30 end |
|
31 |
|
32 if nargin<3 |
|
33 elevkm=0; |
|
34 end |
|
35 |
|
36 |
|
37 |
|
38 [x,y,z] = shval3(flat,flon,elevkm,nmax,gh); |
|
39 [d,i,h,f]=dihf (x,y,z); |
|
40 |
|
41 dev=d*180/pi; |
|
42 |
|
43 if nargout<1 |
|
44 disp([' model ',igrf,' harmonics: ',int2str(nmax)]) |
|
45 disp([' lat: ',num2str(flat)]) |
|
46 disp([' lon: ',num2str(flon)]) |
|
47 disp([' mag dev [deg]: ',num2str(d*180/pi)]) |
|
48 end |
|
49 |
|
50 %===================================================================== |
|
51 |
|
52 function [x,y,z] = shval3(flat,flon,elevkm,nmax,gh) |
|
53 % function [x,y,z] = shval3(flat,flon,elevkm,nmax,gh) |
|
54 % ================================================================ |
|
55 % |
|
56 % version 1.01 |
|
57 % |
|
58 % calculates field components from spherical harmonic (sh) |
|
59 % models. |
|
60 % |
|
61 % input: |
|
62 % flat - north latitude, in degrees |
|
63 % flon - east longitude, in degrees |
|
64 % elevkm - elevation above mean sea level |
|
65 % nmax - maximum degree and order of coefficients |
|
66 % gh - schmidt quasi-normal internal spherical |
|
67 % harmonic coefficients |
|
68 % iext - external coefficients flag (= 0 if none) |
|
69 % ext - the three 1st-degree external coefficients |
|
70 % (not used if iext = 0) |
|
71 % |
|
72 % output: |
|
73 % x - northward component |
|
74 % y - eastward component |
|
75 % z - vertically-downward component |
|
76 % |
|
77 % based on subroutine 'igrf' by d. r. barraclough and |
|
78 % s. r. c. malin, report no. 71/1, institute of geological |
|
79 % sciences, u.k. |
|
80 % |
|
81 % norman w. peddie, u.s. geological survey, mail stop 964, |
|
82 % federal center, box 25046, denver, colorado 80225 |
|
83 % |
|
84 % ================================================================ |
|
85 % the required sizes of the arrays used in this subroutine |
|
86 % depend on the value of nmax. the minimum dimensions |
|
87 % needed are indicated in the table below. (note that this |
|
88 % version is dimensioned for nmax of 14 or less). |
|
89 % |
|
90 % ================================================================ |
|
91 dtr = pi/180; |
|
92 r=elevkm; |
|
93 erad=6371.2; |
|
94 a2=40680925.; |
|
95 b2=40408588.; |
|
96 |
|
97 |
|
98 |
|
99 slat = sin(flat*dtr); |
|
100 aa = min(89.999,max(-89.999,flat)); |
|
101 clat = cos(aa*dtr); |
|
102 sl(1) = sin(flon*dtr); |
|
103 cl(1) = cos(flon*dtr); |
|
104 |
|
105 x=0.0; |
|
106 y=0.0; |
|
107 z=0.0; |
|
108 sd = 0.0; |
|
109 cd = 1.0; |
|
110 n=0; |
|
111 l=1; |
|
112 m=1; |
|
113 |
|
114 npq = (nmax*(nmax+3))/2; |
|
115 aa = a2*clat*clat; |
|
116 bb = b2*slat*slat; |
|
117 cc = aa+bb; |
|
118 dd = sqrt(cc); |
|
119 r=sqrt(elevkm*(elevkm+2.0*dd)+(a2*aa+b2*bb)/cc); |
|
120 cd = (elevkm+dd)/r; |
|
121 sd = (a2-b2)/dd*slat*clat/r; |
|
122 aa = slat; |
|
123 slat = slat*cd-clat*sd; |
|
124 clat = clat*cd+aa*sd; |
|
125 ratio = erad/r; |
|
126 aa = sqrt(3.0); |
|
127 p(1) = 2.0*slat; |
|
128 p(2) = 2.0*clat; |
|
129 p(3) = 4.5*slat*slat-1.5; |
|
130 p(4) = 3.0*aa*clat*slat; |
|
131 q(1) = -clat; |
|
132 q(2) = slat; |
|
133 q(3) = -3.0*clat*slat; |
|
134 q(4) = aa*(slat*slat-clat*clat); |
|
135 |
|
136 |
|
137 for k = 1: npq |
|
138 |
|
139 if (n<m) |
|
140 m=0; |
|
141 n=n+1; |
|
142 rr = ratio^(n+2); |
|
143 fn=n; |
|
144 end; |
|
145 |
|
146 fm=m; |
|
147 |
|
148 if (k>=5) |
|
149 if (m==n) |
|
150 aa = sqrt(1.0-.5/fm); |
|
151 j=k-n-1; |
|
152 p(k) = (1.0+1.0/fm)*aa*clat*p(j); |
|
153 q(k) = aa*(clat*q(j)+slat/fm*p(j)); |
|
154 sl(m) = sl(m-1)*cl(1)+cl(m-1)*sl(1); |
|
155 cl(m) = cl(m-1)*cl(1)-sl(m-1)*sl(1); |
|
156 else |
|
157 aa = sqrt(fn*fn-fm*fm); |
|
158 bb = sqrt((fn-1.0)^2-fm*fm)/aa; |
|
159 cc = (2.0*fn-1.0)/aa; |
|
160 i=k-n; |
|
161 j=k-2*n+1; |
|
162 p(k) = (fn+1.0)*(cc*slat/fn*p(i)-bb/(fn-1.0)*p(j)); |
|
163 q(k) = cc*(slat*q(i)-clat/fn*p(i))-bb*q(j); |
|
164 end; |
|
165 end; |
|
166 |
|
167 aa = rr*gh(l); |
|
168 |
|
169 if (m==0) |
|
170 x=x+aa*q(k); |
|
171 z=z-aa*p(k); |
|
172 l=l+1; |
|
173 else |
|
174 bb = rr*gh(l+1); |
|
175 cc = aa*cl(m)+bb*sl(m); |
|
176 x=x+cc*q(k); |
|
177 z=z-cc*p(k); |
|
178 if (clat>0.0) |
|
179 y=y+(aa*sl(m)-bb*cl(m))*fm*p(k)/((fn+1.0)*clat); |
|
180 else |
|
181 y=y+(aa*sl(m)-bb*cl(m))*q(k)*slat; |
|
182 end; |
|
183 l=l+2; |
|
184 end; |
|
185 |
|
186 m=m+1; |
|
187 |
|
188 end; |
|
189 |
|
190 aa=x; |
|
191 x=x*cd+z*sd; |
|
192 z=z*cd-aa*sd; |
|
193 |
|
194 return; |
|
195 |
|
196 %================================================================== |
|
197 function [d,i,h,f] = dihf(x,y,z) |
|
198 % function [d,i,h,f] = dihf(x,y,z) |
|
199 % =============================================================== |
|
200 % |
|
201 % version 1.01 |
|
202 % |
|
203 % computes the geomagnetic elements d, i, h, and f from |
|
204 % x, y, and z. |
|
205 % |
|
206 % input: |
|
207 % x - northward component |
|
208 % y - eastward component |
|
209 % z - vertically-downward component |
|
210 % |
|
211 % output: |
|
212 % d - declination |
|
213 % i - inclination |
|
214 % h - horizontal intensity |
|
215 % f - total intensity |
|
216 % |
|
217 % a. zunde |
|
218 % usgs, ms 964, box 25046 federal center, denver, co 80225 |
|
219 % |
|
220 % =============================================================== |
|
221 sn=0.0001; |
|
222 |
|
223 % --------------------------------------------------------------- |
|
224 % if d and i cannot be determined, set equal to NaN |
|
225 % --------------------------------------------------------------- |
|
226 h=sqrt(x*x+y*y); |
|
227 f=sqrt(x*x+y*y+z*z); |
|
228 if (f<sn) |
|
229 d=NaN; |
|
230 i=NaN; |
|
231 else |
|
232 i=atan2(z,h); |
|
233 if (h<sn) |
|
234 d=NaN; |
|
235 else |
|
236 hpx = h+x; |
|
237 if (hpx<200) |
|
238 d=pi; |
|
239 else |
|
240 d=2.0*atan2(y,hpx); |
|
241 end; |
|
242 end; |
|
243 end; |
|
244 return; |
|
245 |