/[PAMELA software]/yodaUtility/igrf/shellig.for
ViewVC logotype

Contents of /yodaUtility/igrf/shellig.for

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.1 - (show annotations) (download)
Sat Jun 10 06:33:31 2006 UTC (18 years, 6 months ago) by kusanagi
Branch: MAIN
CVS Tags: yodaUtility2_0/00, yodaUtility2_2/00, yodaUtility2_1/00, HEAD
Software for geomagnetic field calculation. The validity time interval is 2000-2010.

1 C SHELLIG.FOR, Version 2.0, January 1992
2 C
3 C 11/01/91-DKB- SHELLG: lowest starting point for B0 search is 2
4 C 1/27/92-DKB- Adopted to IGRF-91 coeffcients model
5 C 2/05/92-DKB- Reduce variable-names: INTER(P)SHC,EXTRA(P)SHC,INITI(ALI)ZE
6 C 8/08/95-DKB- Updated to IGRF-45-95; new coeff. DGRF90, IGRF95, IGRF95S
7 C 5/31/00-DKB- Updated to IGRF-45-00; new coeff.: IGRF00, IGRF00s
8 C 3/24/05-DKB- Updated to IGRF-45-10; new coeff.: IGRF05, IGRF05s
9 C
10 C*********************************************************************
11 C SUBROUTINES FINDB0, SHELLG, STOER, FELDG, FELDCOF, GETSHC, *
12 C INTERSHC, EXTRASHC, INITIZE *
13 C*********************************************************************
14 C*********************************************************************
15 C
16 C
17 SUBROUTINE FINDB0(STPS,BDEL,VALUE,BEQU,RR0)
18 C--------------------------------------------------------------------
19 C FINDS SMALLEST MAGNETIC FIELD STRENGTH ON FIELD LINE
20 C
21 C INPUT: STPS STEP SIZE FOR FIELD LINE TRACING
22 C COMMON/FIDB0/
23 C SP DIPOLE ORIENTED COORDINATES FORM SHELLG; P(1,*),
24 C P(2,*), P(3,*) CLOSEST TO MAGNETIC EQUATOR
25 C BDEL REQUIRED ACCURACY = [ B(LAST) - BEQU ] / BEQU
26 C B(LAST) IS FIELD STRENGTH BEFORE BEQU
27 C
28 C OUTPUT: VALUE =.FALSE., IF BEQU IS NOT MINIMAL VALUE ON FIELD LINE
29 C BEQU MAGNETIC FIELD STRENGTH AT MAGNETIC EQUATOR
30 C RR0 EQUATORIAL RADIUS NORMALIZED TO EARTH RADIUS
31 C BDEL FINAL ACHIEVED ACCURACY
32 C--------------------------------------------------------------------
33 DIMENSION P(8,4),SP(3)
34 LOGICAL VALUE
35 COMMON/FIDB0/ SP
36 C
37 STEP=STPS
38 IRUN=0
39 7777 IRUN=IRUN+1
40 IF(IRUN.GT.5) THEN
41 VALUE=.FALSE.
42 GOTO 8888
43 ENDIF
44 C*********************FIRST THREE POINTS
45 P(1,2)=SP(1)
46 P(2,2)=SP(2)
47 P(3,2)=SP(3)
48 STEP=-SIGN(STEP,P(3,2))
49 CALL STOER(P(1,2),BQ2,R2)
50 P(1,3)=P(1,2)+0.5*STEP*P(4,2)
51 P(2,3)=P(2,2)+0.5*STEP*P(5,2)
52 P(3,3)=P(3,2)+0.5*STEP
53 CALL STOER(P(1,3),BQ3,R3)
54 P(1,1)=P(1,2)-STEP*(2.*P(4,2)-P(4,3))
55 P(2,1)=P(2,2)-STEP*(2.*P(5,2)-P(5,3))
56 P(3,1)=P(3,2)-STEP
57 CALL STOER(P(1,1),BQ1,R1)
58 P(1,3)=P(1,2)+STEP*(20.*P(4,3)-3.*P(4,2)+P(4,1))/18.
59 P(2,3)=P(2,2)+STEP*(20.*P(5,3)-3.*P(5,2)+P(5,1))/18.
60 P(3,3)=P(3,2)+STEP
61 CALL STOER(P(1,3),BQ3,R3)
62 C******************INVERT SENSE IF REQUIRED
63 IF(BQ3.LE.BQ1) GOTO 2
64 STEP=-STEP
65 R3=R1
66 BQ3=BQ1
67 DO 1 I=1,5
68 ZZ=P(I,1)
69 P(I,1)=P(I,3)
70 1 P(I,3)=ZZ
71 C******************INITIALIZATION
72 2 STEP12=STEP/12.
73 VALUE=.TRUE.
74 BMIN=1.E4
75 BOLD=1.E4
76 C******************CORRECTOR (FIELD LINE TRACING)
77 N=0
78 5555 P(1,3)=P(1,2)+STEP12*(5.*P(4,3)+8.*P(4,2)-P(4,1))
79 N=N+1
80 P(2,3)=P(2,2)+STEP12*(5.*P(5,3)+8.*P(5,2)-P(5,1))
81 C******************PREDICTOR (FIELD LINE TRACING)
82 P(1,4)=P(1,3)+STEP12*(23.*P(4,3)-16.*P(4,2)+5.*P(4,1))
83 P(2,4)=P(2,3)+STEP12*(23.*P(5,3)-16.*P(5,2)+5.*P(5,1))
84 P(3,4)=P(3,3)+STEP
85 CALL STOER(P(1,4),BQ3,R3)
86 DO 1111 J=1,3
87 DO 1111 I=1,8
88 1111 P(I,J)=P(I,J+1)
89 B=SQRT(BQ3)
90 IF(B.LT.BMIN) BMIN=B
91 IF(B.LE.BOLD) THEN
92 BOLD=B
93 ROLD=1./R3
94 SP(1)=P(1,4)
95 SP(2)=P(2,4)
96 SP(3)=P(3,4)
97 GOTO 5555
98 ENDIF
99 IF(BOLD.NE.BMIN) THEN
100 VALUE=.FALSE.
101 ENDIF
102 BDELTA=(B-BOLD)/BOLD
103 IF(BDELTA.GT.BDEL) THEN
104 STEP=STEP/10.
105 GOTO 7777
106 ENDIF
107 8888 RR0=ROLD
108 BEQU=BOLD
109 BDEL=BDELTA
110 RETURN
111 END
112 C
113 C
114 SUBROUTINE SHELLG(GLAT,GLON,ALT,DIMO,FL,ICODE,B0)
115 C--------------------------------------------------------------------
116 C CALCULATES L-VALUE FOR SPECIFIED GEODAETIC COORDINATES, ALTITUDE
117 C AND GEMAGNETIC FIELD MODEL.
118 C REF: G. KLUGE, EUROPEAN SPACE OPERATIONS CENTER, INTERNAL NOTE
119 C NO. 67, 1970.
120 C G. KLUGE, COMPUTER PHYSICS COMMUNICATIONS 3, 31-35, 1972
121 C--------------------------------------------------------------------
122 C CHANGES (D. BILITZA, NOV 87):
123 C - USING CORRECT DIPOL MOMENT I.E.,DIFFERENT COMMON/MODEL/
124 C - USING IGRF EARTH MAGNETIC FIELD MODELS FROM 1945 TO 1990
125 C--------------------------------------------------------------------
126 C INPUT: ENTRY POINT SHELLG
127 C GLAT GEODETIC LATITUDE IN DEGREES (NORTH)
128 C GLON GEODETIC LONGITUDE IN DEGREES (EAST)
129 C ALT ALTITUDE IN KM ABOVE SEA LEVEL
130 C
131 C ENTRY POINT SHELLC
132 C V(3) CARTESIAN COORDINATES IN EARTH RADII (6371.2 KM)
133 C X-AXIS POINTING TO EQUATOR AT 0 LONGITUDE
134 C Y-AXIS POINTING TO EQUATOR AT 90 LONG.
135 C Z-AXIS POINTING TO NORTH POLE
136 C
137 C DIMO DIPOL MOMENT IN GAUSS (NORMALIZED TO EARTH RADIUS)
138 C
139 C COMMON
140 C X(3) NOT USED
141 C H(144) FIELD MODEL COEFFICIENTS ADJUSTED FOR SHELLG
142 C-----------------------------------------------------------------------
143 C OUTPUT: FL L-VALUE
144 C ICODE =1 NORMAL COMPLETION
145 C =2 UNPHYSICAL CONJUGATE POINT (FL MEANINGLESS)
146 C =3 SHELL PARAMETER GREATER THAN LIMIT UP TO
147 C WHICH ACCURATE CALCULATION IS REQUIRED;
148 C APPROXIMATION IS USED.
149 C B0 MAGNETIC FIELD STRENGTH IN GAUSS
150 C-----------------------------------------------------------------------
151 DIMENSION V(3),U(3,3),P(8,100),SP(3)
152 COMMON X(3),H(144)
153 COMMON/FIDB0/ SP
154 COMMON/GENER/ UMR,ERA,AQUAD,BQUAD
155 C
156 C-- RMIN, RMAX ARE BOUNDARIES FOR IDENTIFICATION OF ICODE=2 AND 3
157 C-- STEP IS STEP SIZE FOR FIELD LINE TRACING
158 C-- STEQ IS STEP SIZE FOR INTEGRATION
159 C
160 DATA RMIN,RMAX /0.05,1.01/
161 DATA STEP,STEQ /0.20,0.03/
162 BEQU=1.E10
163 C
164 C*****ENTRY POINT SHELLG TO BE USED WITH GEODETIC CO-ORDINATES
165 RLAT=GLAT*UMR
166 CT=SIN(RLAT)
167 ST=COS(RLAT)
168 D=SQRT(AQUAD-(AQUAD-BQUAD)*CT*CT)
169 X(1)=(ALT+AQUAD/D)*ST/ERA
170 X(3)=(ALT+BQUAD/D)*CT/ERA
171 RLON=GLON*UMR
172 X(2)=X(1)*SIN(RLON)
173 X(1)=X(1)*COS(RLON)
174 GOTO9
175 C
176 C*****ENTRY POINT SHELLC TO BE USED WITH CARTESIAN CO-ORDINATES
177 ENTRY SHELLC(V,FL,B0)
178 X(1)=V(1)
179 X(2)=V(2)
180 X(3)=V(3)
181 C*****CONVERT TO DIPOL-ORIENTED CO-ORDINATES
182 DATA U/ +0.3511737,-0.9148385,-0.1993679,
183 A +0.9335804,+0.3583680,+0.0000000,
184 B +0.0714471,-0.1861260,+0.9799247/
185 9 RQ=1./(X(1)*X(1)+X(2)*X(2)+X(3)*X(3))
186 R3H=SQRT(RQ*SQRT(RQ))
187 P(1,2)=(X(1)*U(1,1)+X(2)*U(2,1)+X(3)*U(3,1))*R3H
188 P(2,2)=(X(1)*U(1,2)+X(2)*U(2,2) )*R3H
189 P(3,2)=(X(1)*U(1,3)+X(2)*U(2,3)+X(3)*U(3,3))*RQ
190 C*****FIRST THREE POINTS OF FIELD LINE
191 STEP=-SIGN(STEP,P(3,2))
192 CALL STOER(P(1,2),BQ2,R2)
193 B0=SQRT(BQ2)
194 P(1,3)=P(1,2)+0.5*STEP*P(4,2)
195 P(2,3)=P(2,2)+0.5*STEP*P(5,2)
196 P(3,3)=P(3,2)+0.5*STEP
197 CALL STOER(P(1,3),BQ3,R3)
198 P(1,1)=P(1,2)-STEP*(2.*P(4,2)-P(4,3))
199 P(2,1)=P(2,2)-STEP*(2.*P(5,2)-P(5,3))
200 P(3,1)=P(3,2)-STEP
201 CALL STOER(P(1,1),BQ1,R1)
202 P(1,3)=P(1,2)+STEP*(20.*P(4,3)-3.*P(4,2)+P(4,1))/18.
203 P(2,3)=P(2,2)+STEP*(20.*P(5,3)-3.*P(5,2)+P(5,1))/18.
204 P(3,3)=P(3,2)+STEP
205 CALL STOER(P(1,3),BQ3,R3)
206 C*****INVERT SENSE IF REQUIRED
207 IF(BQ3.LE.BQ1)GOTO2
208 STEP=-STEP
209 R3=R1
210 BQ3=BQ1
211 DO 1 I=1,7
212 ZZ=P(I,1)
213 P(I,1)=P(I,3)
214 1 P(I,3)=ZZ
215 C*****SEARCH FOR LOWEST MAGNETIC FIELD STRENGTH
216 2 IF(BQ1.LT.BEQU) THEN
217 BEQU=BQ1
218 IEQU=1
219 ENDIF
220 IF(BQ2.LT.BEQU) THEN
221 BEQU=BQ2
222 IEQU=2
223 ENDIF
224 IF(BQ3.LT.BEQU) THEN
225 BEQU=BQ3
226 IEQU=3
227 ENDIF
228 C*****INITIALIZATION OF INTEGRATION LOOPS
229 STEP12=STEP/12.
230 STEP2=STEP+STEP
231 STEQ=SIGN(STEQ,STEP)
232 FI=0.
233 ICODE=1
234 ORADIK=0.
235 OTERM=0.
236 STP=R2*STEQ
237 Z=P(3,2)+STP
238 STP=STP/0.75
239 P(8,1)=STEP2*(P(1,1)*P(4,1)+P(2,1)*P(5,1))
240 P(8,2)=STEP2*(P(1,2)*P(4,2)+P(2,2)*P(5,2))
241 C*****MAIN LOOP (FIELD LINE TRACING)
242 DO 3 N=3,3333
243 C*****CORRECTOR (FIELD LINE TRACING)
244 P(1,N)=P(1,N-1)+STEP12*(5.*P(4,N)+8.*P(4,N-1)-P(4,N-2))
245 P(2,N)=P(2,N-1)+STEP12*(5.*P(5,N)+8.*P(5,N-1)-P(5,N-2))
246 C*****PREPARE EXPANSION COEFFICIENTS FOR INTERPOLATION
247 C*****OF SLOWLY VARYING QUANTITIES
248 P(8,N)=STEP2*(P(1,N)*P(4,N)+P(2,N)*P(5,N))
249 C0=P(1,N-1)**2+P(2,N-1)**2
250 C1=P(8,N-1)
251 C2=(P(8,N)-P(8,N-2))*0.25
252 C3=(P(8,N)+P(8,N-2)-C1-C1)/6.0
253 D0=P(6,N-1)
254 D1=(P(6,N)-P(6,N-2))*0.5
255 D2=(P(6,N)+P(6,N-2)-D0-D0)*0.5
256 E0=P(7,N-1)
257 E1=(P(7,N)-P(7,N-2))*0.5
258 E2=(P(7,N)+P(7,N-2)-E0-E0)*0.5
259 C*****INNER LOOP (FOR QUADRATURE)
260 4 T=(Z-P(3,N-1))/STEP
261 IF(T.GT.1.)GOTO5
262 HLI=0.5*(((C3*T+C2)*T+C1)*T+C0)
263 ZQ=Z*Z
264 R=HLI+SQRT(HLI*HLI+ZQ)
265 IF(R.LE.RMIN)GOTO30
266 RQ=R*R
267 FF=SQRT(1.+3.*ZQ/RQ)
268 RADIK=B0-((D2*T+D1)*T+D0)*R*RQ*FF
269 IF(R-RMAX)44,44,45
270 45 ICODE=2
271 RADIK=RADIK-12.*(R-RMAX)**2
272 44 IF(RADIK+RADIK.LE.ORADIK) GOTO 10
273 TERM=SQRT(RADIK)*FF*((E2*T+E1)*T+E0)/(RQ+ZQ)
274 FI=FI+STP*(OTERM+TERM)
275 ORADIK=RADIK
276 OTERM=TERM
277 STP=R*STEQ
278 Z=Z+STP
279 GOTO4
280 C*****PREDICTOR (FIELD LINE TRACING)
281 5 P(1,N+1)=P(1,N)+STEP12*(23.*P(4,N)-16.*P(4,N-1)+5.*P(4,N-2))
282 P(2,N+1)=P(2,N)+STEP12*(23.*P(5,N)-16.*P(5,N-1)+5.*P(5,N-2))
283 P(3,N+1)=P(3,N)+STEP
284 CALL STOER(P(1,N+1),BQ3,R3)
285 C*****SEARCH FOR LOWEST MAGNETIC FIELD STRENGTH
286 IF(BQ3.LT.BEQU) THEN
287 IEQU=N+1
288 BEQU=BQ3
289 ENDIF
290 3 CONTINUE
291 10 IF(IEQU.lt.2) IEQU=2
292 SP(1)=P(1,IEQU-1)
293 SP(2)=P(2,IEQU-1)
294 SP(3)=P(3,IEQU-1)
295 IF(ORADIK.LT.1E-15)GOTO11
296 FI=FI+STP/0.75*OTERM*ORADIK/(ORADIK-RADIK)
297 C
298 C-- The minimal allowable value of FI was changed from 1E-15 to 1E-12,
299 C-- because 1E-38 is the minimal allowable arg. for ALOG in our envir.
300 C-- D. Bilitza, Nov 87.
301 C
302 11 FI=0.5*ABS(FI)/SQRT(B0)+1E-12
303 C*****COMPUTE L FROM B AND I. SAME AS CARMEL IN INVAR.
304 C
305 C-- Correct dipole moment is used here. D. Bilitza, Nov 87.
306 C
307 DIMOB0=DIMO/B0
308 arg1=alog(FI)
309 arg2=alog(DIMOB0)
310 c arg = FI*FI*FI/DIMOB0
311 c if(abs(arg).gt.88.0) arg=88.0
312 XX=3*arg1-arg2
313 IF(XX.GT.23.0) GOTO 776
314 IF(XX.GT.11.7) GOTO 775
315 IF(XX.GT.+3.0) GOTO 774
316 IF(XX.GT.-3.0) GOTO 773
317 IF(XX.GT.-22.) GOTO 772
318 771 GG=3.33338E-1*XX+3.0062102E-1
319 GOTO777
320 772 GG=((((((((-8.1537735E-14*XX+8.3232531E-13)*XX+1.0066362E-9)*XX+
321 18.1048663E-8)*XX+3.2916354E-6)*XX+8.2711096E-5)*XX+1.3714667E-3)*
322 2XX+1.5017245E-2)*XX+4.3432642E-1)*XX+6.2337691E-1
323 GOTO777
324 773 GG=((((((((2.6047023E-10*XX+2.3028767E-9)*XX-2.1997983E-8)*XX-
325 15.3977642E-7)*XX-3.3408822E-6)*XX+3.8379917E-5)*XX+1.1784234E-3)*
326 2XX+1.4492441E-2)*XX+4.3352788E-1)*XX+6.228644E-1
327 GOTO777
328 774 GG=((((((((6.3271665E-10*XX-3.958306E-8)*XX+9.9766148E-07)*XX-
329 11.2531932E-5)*XX+7.9451313E-5)*XX-3.2077032E-4)*XX+2.1680398E-3)*
330 2XX+1.2817956E-2)*XX+4.3510529E-1)*XX+6.222355E-1
331 GOTO777
332 775 GG=(((((2.8212095E-8*XX-3.8049276E-6)*XX+2.170224E-4)*XX-6.7310339
333 1E-3)*XX+1.2038224E-1)*XX-1.8461796E-1)*XX+2.0007187E0
334 GOTO777
335 776 GG=XX-3.0460681E0
336 777 FL=EXP(ALOG((1.+EXP(GG))*DIMOB0)/3.0)
337 RETURN
338 C*****APPROXIMATION FOR HIGH VALUES OF L.
339 30 ICODE=3
340 T=-P(3,N-1)/STEP
341 FL=1./(ABS(((C3*T+C2)*T+C1)*T+C0)+1E-15)
342 RETURN
343 END
344 C
345 C
346 SUBROUTINE STOER(P,BQ,R)
347 C*******************************************************************
348 C* SUBROUTINE USED FOR FIELD LINE TRACING IN SHELLG *
349 C* CALLS ENTRY POINT FELDI IN GEOMAGNETIC FIELD SUBROUTINE FELDG *
350 C*******************************************************************
351 DIMENSION P(7),U(3,3)
352 COMMON XI(3),H(144)
353 C*****XM,YM,ZM ARE GEOMAGNETIC CARTESIAN INVERSE CO-ORDINATES
354 ZM=P(3)
355 FLI=P(1)*P(1)+P(2)*P(2)+1E-15
356 R=0.5*(FLI+SQRT(FLI*FLI+(ZM+ZM)**2))
357 RQ=R*R
358 WR=SQRT(R)
359 XM=P(1)*WR
360 YM=P(2)*WR
361 C*****TRANSFORM TO GEOGRAPHIC CO-ORDINATE SYSTEM
362 DATA U/ +0.3511737,-0.9148385,-0.1993679,
363 A +0.9335804,+0.3583680,+0.0000000,
364 B +0.0714471,-0.1861260,+0.9799247/
365 XI(1)=XM*U(1,1)+YM*U(1,2)+ZM*U(1,3)
366 XI(2)=XM*U(2,1)+YM*U(2,2)+ZM*U(2,3)
367 XI(3)=XM*U(3,1) +ZM*U(3,3)
368 C*****COMPUTE DERIVATIVES
369 c Changed from CALL FELDI(XI,H); XI, H are in COMMON block; results
370 c are the same; dkb Feb 1998
371 CALL FELDI
372 Q=H(1)/RQ
373 DX=H(3)+H(3)+Q*XI(1)
374 DY=H(4)+H(4)+Q*XI(2)
375 DZ=H(2)+H(2)+Q*XI(3)
376 C*****TRANSFORM BACK TO GEOMAGNETIC CO-ORDINATE SYSTEM
377 DXM=U(1,1)*DX+U(2,1)*DY+U(3,1)*DZ
378 DYM=U(1,2)*DX+U(2,2)*DY
379 DZM=U(1,3)*DX+U(2,3)*DY+U(3,3)*DZ
380 DR=(XM*DXM+YM*DYM+ZM*DZM)/R
381 C*****FORM SLOWLY VARYING EXPRESSIONS
382 P(4)=(WR*DXM-0.5*P(1)*DR)/(R*DZM)
383 P(5)=(WR*DYM-0.5*P(2)*DR)/(R*DZM)
384 DSQ=RQ*(DXM*DXM+DYM*DYM+DZM*DZM)
385 BQ=DSQ*RQ*RQ
386 P(6)=SQRT(DSQ/(RQ+3.*ZM*ZM))
387 P(7)=P(6)*(RQ+ZM*ZM)/(RQ*DZM)
388 RETURN
389 END
390 C
391 C
392 SUBROUTINE FELDG(GLAT,GLON,ALT,BNORTH,BEAST,BDOWN,BABS)
393 C-------------------------------------------------------------------
394 C CALCULATES EARTH MAGNETIC FIELD FROM SPHERICAL HARMONICS MODEL
395 C REF: G. KLUGE, EUROPEAN SPACE OPERATIONS CENTRE, INTERNAL NOTE 61,
396 C 1970.
397 C--------------------------------------------------------------------
398 C CHANGES (D. BILITZA, NOV 87):
399 C - FIELD COEFFICIENTS IN BINARY DATA FILES INSTEAD OF BLOCK DATA
400 C - CALCULATES DIPOL MOMENT
401 C--------------------------------------------------------------------
402 C INPUT: ENTRY POINT FELDG
403 C GLAT GEODETIC LATITUDE IN DEGREES (NORTH)
404 C GLON GEODETIC LONGITUDE IN DEGREES (EAST)
405 C ALT ALTITUDE IN KM ABOVE SEA LEVEL
406 C
407 C ENTRY POINT FELDC
408 C V(3) CARTESIAN COORDINATES IN EARTH RADII (6371.2 KM)
409 C X-AXIS POINTING TO EQUATOR AT 0 LONGITUDE
410 C Y-AXIS POINTING TO EQUATOR AT 90 LONG.
411 C Z-AXIS POINTING TO NORTH POLE
412 C
413 C COMMON BLANK AND ENTRY POINT FELDI ARE NEEDED WHEN USED
414 C IN CONNECTION WITH L-CALCULATION PROGRAM SHELLG.
415 C
416 C COMMON /MODEL/ AND /GENER/
417 C UMR = ATAN(1.0)*4./180. <DEGREE>*UMR=<RADIANT>
418 C ERA EARTH RADIUS FOR NORMALIZATION OF CARTESIAN
419 C COORDINATES (6371.2 KM)
420 C AQUAD, BQUAD SQUARE OF MAJOR AND MINOR HALF AXIS FOR
421 C EARTH ELLIPSOID AS RECOMMENDED BY INTERNATIONAL
422 C ASTRONOMICAL UNION (6378.160, 6356.775 KM).
423 C NMAX MAXIMUM ORDER OF SPHERICAL HARMONICS
424 C TIME YEAR (DECIMAL: 1973.5) FOR WHICH MAGNETIC
425 C FIELD IS TO BE CALCULATED
426 C G(M) NORMALIZED FIELD COEFFICIENTS (SEE FELDCOF)
427 C M=NMAX*(NMAX+2)
428 C------------------------------------------------------------------------
429 C OUTPUT: BABS MAGNETIC FIELD STRENGTH IN GAUSS
430 C BNORTH, BEAST, BDOWN COMPONENTS OF THE FIELD WITH RESPECT
431 C TO THE LOCAL GEODETIC COORDINATE SYSTEM, WITH AXIS
432 C POINTING IN THE TANGENTIAL PLANE TO THE NORTH, EAST
433 C AND DOWNWARD.
434 C-----------------------------------------------------------------------
435 DIMENSION V(3),B(3)
436 CHARACTER*12 NAME
437 COMMON XI(3),H(144)
438 COMMON/MODEL/ NAME,NMAX,TIME,G(144)
439 COMMON/GENER/ UMR,ERA,AQUAD,BQUAD
440 C
441 C-- IS RECORDS ENTRY POINT
442 C
443 C*****ENTRY POINT FELDG TO BE USED WITH GEODETIC CO-ORDINATES
444 IS=1
445 RLAT=GLAT*UMR
446 CT=SIN(RLAT)
447 ST=COS(RLAT)
448 D=SQRT(AQUAD-(AQUAD-BQUAD)*CT*CT)
449 RLON=GLON*UMR
450 CP=COS(RLON)
451 SP=SIN(RLON)
452 ZZZ=(ALT+BQUAD/D)*CT/ERA
453 RHO=(ALT+AQUAD/D)*ST/ERA
454 XXX=RHO*CP
455 YYY=RHO*SP
456 GOTO10
457 C
458 C*****ENTRY POINT FELDC TO BE USED WITH CARTESIAN CO-ORDINATES
459 ENTRY FELDC(V,B)
460 IS=2
461 XXX=V(1)
462 YYY=V(2)
463 ZZZ=V(3)
464 10 RQ=1./(XXX*XXX+YYY*YYY+ZZZ*ZZZ)
465 XI(1)=XXX*RQ
466 XI(2)=YYY*RQ
467 XI(3)=ZZZ*RQ
468 GOTO20
469 C
470 C*****ENTRY POINT FELDI USED FOR L COMPUTATION
471 ENTRY FELDI
472 IS=3
473 20 IHMAX=NMAX*NMAX+1
474 LAST=IHMAX+NMAX+NMAX
475 IMAX=NMAX+NMAX-1
476 DO 8 I=IHMAX,LAST
477 8 H(I)=G(I)
478 DO 6 K=1,3,2
479 I=IMAX
480 IH=IHMAX
481 1 IL=IH-I
482 F=2./FLOAT(I-K+2)
483 X=XI(1)*F
484 Y=XI(2)*F
485 Z=XI(3)*(F+F)
486 I=I-2
487 IF(I-1)5,4,2
488 2 DO 3 M=3,I,2
489 H(IL+M+1)=G(IL+M+1)+Z*H(IH+M+1)+X*(H(IH+M+3)-H(IH+M-1))
490 A -Y*(H(IH+M+2)+H(IH+M-2))
491 3 H(IL+M)=G(IL+M)+Z*H(IH+M)+X*(H(IH+M+2)-H(IH+M-2))
492 A +Y*(H(IH+M+3)+H(IH+M-1))
493 4 H(IL+2)=G(IL+2)+Z*H(IH+2)+X*H(IH+4)-Y*(H(IH+3)+H(IH))
494 H(IL+1)=G(IL+1)+Z*H(IH+1)+Y*H(IH+4)+X*(H(IH+3)-H(IH))
495 5 H(IL)=G(IL)+Z*H(IH)+2.*(X*H(IH+1)+Y*H(IH+2))
496 IH=IL
497 IF(I.GE.K)GOTO1
498 6 CONTINUE
499 IF(IS.EQ.3)RETURN
500 S=.5*H(1)+2.*(H(2)*XI(3)+H(3)*XI(1)+H(4)*XI(2))
501 T=(RQ+RQ)*SQRT(RQ)
502 BXXX=T*(H(3)-S*XXX)
503 BYYY=T*(H(4)-S*YYY)
504 BZZZ=T*(H(2)-S*ZZZ)
505 IF(IS.EQ.2)GOTO7
506 BABS=SQRT(BXXX*BXXX+BYYY*BYYY+BZZZ*BZZZ)
507 BEAST=BYYY*CP-BXXX*SP
508 BRHO=BYYY*SP+BXXX*CP
509 BNORTH=BZZZ*ST-BRHO*CT
510 BDOWN=-BZZZ*CT-BRHO*ST
511 RETURN
512 7 B(1)=BXXX
513 B(2)=BYYY
514 B(3)=BZZZ
515 RETURN
516 END
517 C
518 C
519 SUBROUTINE FELDCOF(YEAR,DIMO)
520 C------------------------------------------------------------------------
521 C DETERMINES COEFFICIENTS AND DIPOL MOMENT FROM IGRF MODELS
522 C
523 C INPUT: YEAR DECIMAL YEAR FOR WHICH GEOMAGNETIC FIELD IS TO
524 C BE CALCULATED (e.g.:1995.5 for day 185 of 1995)
525 C OUTPUT: DIMO GEOMAGNETIC DIPOL MOMENT IN GAUSS (NORMALIZED
526 C TO EARTH'S RADIUS) AT THE TIME (YEAR)
527 C D. BILITZA, NSSDC, GSFC, CODE 633, GREENBELT, MD 20771,
528 C (301)286-9536 NOV 1987.
529 C -corrected for 2000 update - dkb- 5/31/2000
530 C ### updated to IGRF-2000 version -dkb- 5/31/2000
531 C ### updated to IGRF-2005 version -dkb- 3/24/2000
532 C-----------------------------------------------------------------------
533 CHARACTER*12 FILMOD, FIL1, FIL2
534 C ### FILMOD, DTEMOD arrays +1
535 DIMENSION GH1(144),GH2(120),GHA(144),FILMOD(14),DTEMOD(14)
536 DOUBLE PRECISION X,F0,F
537 COMMON/MODEL/ FIL1,NMAX,TIME,GH1
538 COMMON/GENER/ UMR,ERAD,AQUAD,BQUAD
539 C ### changed to conform with IGRF 45-95, also FILMOD, DTEMOD arrays +1
540 DATA FILMOD /'dgrf45.dat', 'dgrf50.dat',
541 1 'dgrf55.dat', 'dgrf60.dat', 'dgrf65.dat',
542 2 'dgrf70.dat', 'dgrf75.dat', 'dgrf80.dat',
543 3 'dgrf85.dat', 'dgrf90.dat', 'dgrf95.dat',
544 4 'dgrf00.dat','igrf05.dat','igrf05s.dat'/
545 DATA DTEMOD / 1945., 1950., 1955., 1960., 1965., 1970.,
546 1 1975., 1980., 1985., 1990., 1995., 2000.,2005.,2010./
547 C
548 C ### numye is number of 5-year priods represented by IGRF
549 C
550 NUMYE=13
551 C
552 C IS=0 FOR SCHMIDT NORMALIZATION IS=1 GAUSS NORMALIZATION
553 C IU IS INPUT UNIT NUMBER FOR IGRF COEFFICIENT SETS
554 C
555 IU = 10
556 IS = 0
557 C-- DETERMINE IGRF-YEARS FOR INPUT-YEAR
558 TIME = YEAR
559 IYEA = INT(YEAR/5.)*5
560 L = (IYEA - 1945)/5 + 1
561 IF(L.LT.1) L=1
562 IF(L.GT.NUMYE) L=NUMYE
563 DTE1 = DTEMOD(L)
564 FIL1 = FILMOD(L)
565 DTE2 = DTEMOD(L+1)
566 FIL2 = FILMOD(L+1)
567 C-- GET IGRF COEFFICIENTS FOR THE BOUNDARY YEARS
568 CALL GETSHC (IU, FIL1, NMAX1, ERAD, GH1, IER)
569 IF (IER .NE. 0) STOP
570 CALL GETSHC (IU, FIL2, NMAX2, ERAD, GH2, IER)
571 IF (IER .NE. 0) STOP
572 C-- DETERMINE IGRF COEFFICIENTS FOR YEAR
573 IF (L .LE. NUMYE-1) THEN
574 CALL INTERSHC (YEAR, DTE1, NMAX1, GH1, DTE2,
575 1 NMAX2, GH2, NMAX, GHA)
576 ELSE
577 CALL EXTRASHC (YEAR, DTE1, NMAX1, GH1, NMAX2,
578 1 GH2, NMAX, GHA)
579 ENDIF
580 C-- DETERMINE MAGNETIC DIPOL MOMENT AND COEFFIECIENTS G
581 F0=0.D0
582 DO 1234 J=1,3
583 F = GHA(J) * 1.D-5
584 F0 = F0 + F * F
585 1234 CONTINUE
586 DIMO = DSQRT(F0)
587
588 GH1(1) = 0.0
589 I=2
590 F0=1.D-5
591 IF(IS.EQ.0) F0=-F0
592 SQRT2=SQRT(2.)
593
594 DO 9 N=1,NMAX
595 X = N
596 F0 = F0 * X * X / (4.D0 * X - 2.D0)
597 IF(IS.EQ.0) F0 = F0 * (2.D0 * X - 1.D0) / X
598 F = F0 * 0.5D0
599 IF(IS.EQ.0) F = F * SQRT2
600 GH1(I) = GHA(I-1) * F0
601 I = I+1
602 DO 9 M=1,N
603 F = F * (X + M) / (X - M + 1.D0)
604 IF(IS.EQ.0) F = F * DSQRT((X - M + 1.D0) / (X + M))
605 GH1(I) = GHA(I-1) * F
606 GH1(I+1) = GHA(I) * F
607 I=I+2
608 9 CONTINUE
609 RETURN
610 END
611 C
612 C
613 SUBROUTINE GETSHC (IU, FSPEC, NMAX, ERAD, GH, IER)
614
615 C ===============================================================
616 C
617 C Version 1.01
618 C
619 C Reads spherical harmonic coefficients from the specified
620 C file into an array.
621 C
622 C Input:
623 C IU - Logical unit number
624 C FSPEC - File specification
625 C
626 C Output:
627 C NMAX - Maximum degree and order of model
628 C ERAD - Earth's radius associated with the spherical
629 C harmonic coefficients, in the same units as
630 C elevation
631 C GH - Schmidt quasi-normal internal spherical
632 C harmonic coefficients
633 C IER - Error number: = 0, no error
634 C = -2, records out of order
635 C = FORTRAN run-time error number
636 C
637 C A. Zunde
638 C USGS, MS 964, Box 25046 Federal Center, Denver, CO 80225
639 C
640 C ===============================================================
641
642 CHARACTER FSPEC*(*)
643 DIMENSION GH(*)
644
645 C ---------------------------------------------------------------
646 C Open coefficient file. Read past first header record.
647 C Read degree and order of model and Earth's radius.
648 C ---------------------------------------------------------------
649 OPEN (IU, FILE=FSPEC, STATUS='OLD', IOSTAT=IER, ERR=999)
650 READ (IU, *, IOSTAT=IER, ERR=999)
651 READ (IU, *, IOSTAT=IER, ERR=999) NMAX, ERAD
652 C ---------------------------------------------------------------
653 C Read the coefficient file, arranged as follows:
654 C
655 C N M G H
656 C ----------------------
657 C / 1 0 GH(1) -
658 C / 1 1 GH(2) GH(3)
659 C / 2 0 GH(4) -
660 C / 2 1 GH(5) GH(6)
661 C NMAX*(NMAX+3)/2 / 2 2 GH(7) GH(8)
662 C records \ 3 0 GH(9) -
663 C \ . . . .
664 C \ . . . .
665 C NMAX*(NMAX+2) \ . . . .
666 C elements in GH \ NMAX NMAX . .
667 C
668 C N and M are, respectively, the degree and order of the
669 C coefficient.
670 C ---------------------------------------------------------------
671
672 I = 0
673 DO 2211 NN = 1, NMAX
674 DO 2233 MM = 0, NN
675 READ (IU, *, IOSTAT=IER, ERR=999) N, M, G, H
676 IF (NN .NE. N .OR. MM .NE. M) THEN
677 IER = -2
678 GOTO 999
679 ENDIF
680 I = I + 1
681 GH(I) = G
682 IF (M .NE. 0) THEN
683 I = I + 1
684 GH(I) = H
685 ENDIF
686 2233 CONTINUE
687 2211 CONTINUE
688
689 999 CLOSE (IU)
690
691 RETURN
692 END
693 C
694 C
695 SUBROUTINE INTERSHC (DATE, DTE1, NMAX1, GH1, DTE2,
696 1 NMAX2, GH2, NMAX, GH)
697
698 C ===============================================================
699 C
700 C Version 1.01
701 C
702 C Interpolates linearly, in time, between two spherical
703 C harmonic models.
704 C
705 C Input:
706 C DATE - Date of resulting model (in decimal year)
707 C DTE1 - Date of earlier model
708 C NMAX1 - Maximum degree and order of earlier model
709 C GH1 - Schmidt quasi-normal internal spherical
710 C harmonic coefficients of earlier model
711 C DTE2 - Date of later model
712 C NMAX2 - Maximum degree and order of later model
713 C GH2 - Schmidt quasi-normal internal spherical
714 C harmonic coefficients of later model
715 C
716 C Output:
717 C GH - Coefficients of resulting model
718 C NMAX - Maximum degree and order of resulting model
719 C
720 C A. Zunde
721 C USGS, MS 964, Box 25046 Federal Center, Denver, CO 80225
722 C
723 C ===============================================================
724
725 DIMENSION GH1(*), GH2(*), GH(*)
726
727 C ---------------------------------------------------------------
728 C The coefficients (GH) of the resulting model, at date
729 C DATE, are computed by linearly interpolating between the
730 C coefficients of the earlier model (GH1), at date DTE1,
731 C and those of the later model (GH2), at date DTE2. If one
732 C model is smaller than the other, the interpolation is
733 C performed with the missing coefficients assumed to be 0.
734 C ---------------------------------------------------------------
735
736 FACTOR = (DATE - DTE1) / (DTE2 - DTE1)
737
738 IF (NMAX1 .EQ. NMAX2) THEN
739 K = NMAX1 * (NMAX1 + 2)
740 NMAX = NMAX1
741 ELSE IF (NMAX1 .GT. NMAX2) THEN
742 K = NMAX2 * (NMAX2 + 2)
743 L = NMAX1 * (NMAX1 + 2)
744 DO 1122 I = K + 1, L
745 1122 GH(I) = GH1(I) + FACTOR * (-GH1(I))
746 NMAX = NMAX1
747 ELSE
748 K = NMAX1 * (NMAX1 + 2)
749 L = NMAX2 * (NMAX2 + 2)
750 DO 1133 I = K + 1, L
751 1133 GH(I) = FACTOR * GH2(I)
752 NMAX = NMAX2
753 ENDIF
754
755 DO 1144 I = 1, K
756 1144 GH(I) = GH1(I) + FACTOR * (GH2(I) - GH1(I))
757
758 RETURN
759 END
760 C
761 C
762 SUBROUTINE EXTRASHC (DATE, DTE1, NMAX1, GH1, NMAX2,
763 1 GH2, NMAX, GH)
764
765 C ===============================================================
766 C
767 C Version 1.01
768 C
769 C Extrapolates linearly a spherical harmonic model with a
770 C rate-of-change model.
771 C
772 C Input:
773 C DATE - Date of resulting model (in decimal year)
774 C DTE1 - Date of base model
775 C NMAX1 - Maximum degree and order of base model
776 C GH1 - Schmidt quasi-normal internal spherical
777 C harmonic coefficients of base model
778 C NMAX2 - Maximum degree and order of rate-of-change
779 C model
780 C GH2 - Schmidt quasi-normal internal spherical
781 C harmonic coefficients of rate-of-change model
782 C
783 C Output:
784 C GH - Coefficients of resulting model
785 C NMAX - Maximum degree and order of resulting model
786 C
787 C A. Zunde
788 C USGS, MS 964, Box 25046 Federal Center, Denver, CO 80225
789 C
790 C ===============================================================
791
792 DIMENSION GH1(*), GH2(*), GH(*)
793
794 C ---------------------------------------------------------------
795 C The coefficients (GH) of the resulting model, at date
796 C DATE, are computed by linearly extrapolating the coef-
797 C ficients of the base model (GH1), at date DTE1, using
798 C those of the rate-of-change model (GH2), at date DTE2. If
799 C one model is smaller than the other, the extrapolation is
800 C performed with the missing coefficients assumed to be 0.
801 C ---------------------------------------------------------------
802
803 FACTOR = (DATE - DTE1)
804
805 IF (NMAX1 .EQ. NMAX2) THEN
806 K = NMAX1 * (NMAX1 + 2)
807 NMAX = NMAX1
808 ELSE IF (NMAX1 .GT. NMAX2) THEN
809 K = NMAX2 * (NMAX2 + 2)
810 L = NMAX1 * (NMAX1 + 2)
811 DO 1155 I = K + 1, L
812 1155 GH(I) = GH1(I)
813 NMAX = NMAX1
814 ELSE
815 K = NMAX1 * (NMAX1 + 2)
816 L = NMAX2 * (NMAX2 + 2)
817 DO 1166 I = K + 1, L
818 1166 GH(I) = FACTOR * GH2(I)
819 NMAX = NMAX2
820 ENDIF
821
822 DO 1177 I = 1, K
823 1177 GH(I) = GH1(I) + FACTOR * GH2(I)
824
825 RETURN
826 END
827 C
828 C
829 SUBROUTINE INITIZE
830 C----------------------------------------------------------------
831 C Initializes the parameters in COMMON/GENER/
832 C
833 C UMR = ATAN(1.0)*4./180. <DEGREE>*UMR=<RADIANT>
834 C ERA EARTH RADIUS FOR NORMALIZATION OF CARTESIAN
835 C COORDINATES (6371.2 KM)
836 C EREQU MAJOR HALF AXIS FOR EARTH ELLIPSOID (6378.160 KM)
837 C ERPOL MINOR HALF AXIS FOR EARTH ELLIPSOID (6356.775 KM)
838 C AQUAD SQUARE OF MAJOR HALF AXIS FOR EARTH ELLIPSOID
839 C BQUAD SQUARE OF MINOR HALF AXIS FOR EARTH ELLIPSOID
840 C
841 C ERA, EREQU and ERPOL as recommended by the INTERNATIONAL
842 C ASTRONOMICAL UNION .
843 C-----------------------------------------------------------------
844 COMMON/GENER/ UMR,ERA,AQUAD,BQUAD
845 ERA=6371.2
846 EREQU=6378.16
847 ERPOL=6356.775
848 AQUAD=EREQU*EREQU
849 BQUAD=ERPOL*ERPOL
850 UMR=ATAN(1.0)*4./180.
851 RETURN
852 END

  ViewVC Help
Powered by ViewVC 1.1.23