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