1 function [dev,d,h,i,f,x,y,z]=magdev(flat,flon,elevkm,nmax,igrf); |
1 function [dev,d,h,i,f,x,y,z]=magdev(flat,flon,elevkm,year); |
2 % 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,year); |
3 % |
3 % |
4 % compute magnetic deviation |
4 % compute magnetic deviation |
5 % input: |
5 % |
6 % flat latitude (degree) |
6 % based on IGRF12 model (observed until end of 2014, predicted until end of 2019) |
7 % flon longitude (degree) |
7 % |
8 % elevkm elevation above mean geoid (km) |
8 % input: flat - latitude (degree) |
9 % nmax number of harmonics |
9 % flon - longitude (degree) |
10 % igrf (need file, default 'IGRF95') |
10 % elevkm [0] - elevation above mean geoid (km) |
11 % |
11 % year - decimal year |
12 % output: |
12 % |
13 % dev mag. deviation (degree) |
13 % output: dev - mag. deviation (degree) |
14 % |
14 % |
15 % |
15 % |
16 % based of FORTRAN ROUTINE GEOMAG.FOR |
16 % based of FORTRAN ROUTINE GEOMAG.FOR |
17 % more info under http://fdd.gsfc.nasa.gov/IGRF.html |
17 % more info under http://www.ngdc.noaa.gov/IAGA/vmod/igrf.html |
18 % |
18 % |
|
19 % |
|
20 % version 0.5 last change 13.07.2015 |
|
21 |
|
22 |
19 % M. Visbeck, LDEO FEB 2000 |
23 % M. Visbeck, LDEO FEB 2000 |
20 % |
24 % rewritten for other epochs, G. Krahmann, IFM-GEOMAR Mar 2007 |
21 |
25 % |
22 if nargin<5 |
26 % modified for >=R14 GK, Jun 2007 0.1-->0.2 |
23 igrf='IGRF95'; |
27 % switch to IGRF11 GK, 08.11.2010 0.2-->0.3 |
24 igrf='IGRF00'; |
28 % replaced sind/cosd with sind/cosd GK, 31.05.2011 0.3-->0.4 |
|
29 % switch to IGRF12 GK, 12.07.2015 0.4-->0.5 |
|
30 |
|
31 %====================================================================== |
|
32 % M A G D E V . M |
|
33 % doc: Wed Sep 4 18:17:31 2019 |
|
34 % dlm: Wed Sep 4 18:18:16 2019 |
|
35 % (c) 2019 G. Krahmann |
|
36 % uE-Info: 40 49 NIL 0 0 72 0 2 8 NIL ofnI |
|
37 %====================================================================== |
|
38 |
|
39 % CHANGES BY ANT: |
|
40 % SEp 4, 2019: - changed sin/cos_d to sin/cosd |
|
41 |
|
42 |
|
43 % |
|
44 % read the coefficients |
|
45 % |
|
46 fname = 'igrf12coeffs.xls'; |
|
47 warning off % avoid non-sensical warning in >=R14 |
|
48 gh = xlsread(fname); |
|
49 warning on |
|
50 |
|
51 |
|
52 % |
|
53 % determine the maximum order of polynomials |
|
54 % |
|
55 % this might need modification in future versions of the file |
|
56 % |
|
57 if year<2000 |
|
58 nmax = 10; |
|
59 else |
|
60 nmax = 13; |
25 end |
61 end |
26 eval(['gh=',igrf,';']) |
62 |
27 |
63 |
28 if nargin<4 |
64 % |
29 nmax=10; |
65 % interpolate the coefficients in time |
|
66 % |
|
67 warning off % another workaround for annoying >=R14 warnings |
|
68 if year <= gh(1,end-1) |
|
69 gh2 = interp1( gh(1,3:end-1),gh(2:end,3:end-1)',year)'; |
|
70 else |
|
71 gh2 = interp1( gh(1,3:end-1),gh(2:end,3:end-1)',gh(1,end-1))'; |
|
72 ghc = gh2*0; |
|
73 dummy = gh(2:end,end); |
|
74 good = find(~isnan(dummy)); |
|
75 ghc(good) = dummy(good); |
|
76 gh2 = gh2+ghc*(year-gh(1,end-1)); |
30 end |
77 end |
31 |
78 warning on |
|
79 |
|
80 |
|
81 % |
|
82 % set default elevation |
|
83 % |
32 if nargin<3 |
84 if nargin<3 |
33 elevkm=0; |
85 elevkm=0; |
34 end |
86 end |
35 |
87 |
36 |
88 |
37 |
89 % |
38 [x,y,z] = shval3(flat,flon,elevkm,nmax,gh); |
90 % calculate field |
|
91 % |
|
92 [x,y,z] = shval3(flat,flon,elevkm,nmax,gh2); |
39 [d,i,h,f]=dihf (x,y,z); |
93 [d,i,h,f]=dihf (x,y,z); |
40 |
94 |
|
95 |
|
96 % |
|
97 % convert units |
|
98 % |
41 dev=d*180/pi; |
99 dev=d*180/pi; |
42 |
100 |
|
101 |
|
102 % |
|
103 % more detailed screen output, if output is not stored |
|
104 % |
43 if nargout<1 |
105 if nargout<1 |
44 disp([' model ',igrf,' harmonics: ',int2str(nmax)]) |
106 disp([' model ',fname,' harmonics: ',int2str(nmax)]) |
45 disp([' lat: ',num2str(flat)]) |
107 disp([' lat: ',num2str(flat)]) |
46 disp([' lon: ',num2str(flon)]) |
108 disp([' lon: ',num2str(flon)]) |
|
109 disp([' tim: ',num2str(year)]) |
47 disp([' mag dev [deg]: ',num2str(d*180/pi)]) |
110 disp([' mag dev [deg]: ',num2str(d*180/pi)]) |
48 end |
111 end |
49 |
112 |
50 %===================================================================== |
113 %===================================================================== |
51 |
114 |