-
Notifications
You must be signed in to change notification settings - Fork 2
/
apex_routines.f
2918 lines (2740 loc) · 149 KB
/
apex_routines.f
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
subroutine apex_to_geo(date, aLat, aLon, Alt, gLat, gLon, sLat, sLon)
C Input:
C
C DATE = Year and fraction (1990.0 = 1990 January 1, 0 UT)
C aLat = Apex latitude in degrees
C aLon = Apex longitude in degrees
C ALT = Altitude in km
C
C Output:
C
C gLat = Geographic Latitude
C gLon = Geographin Longitude
PARAMETER (RTOD=5.72957795130823E1,DTOR=1.745329251994330E-2)
PARAMETER (RE=6371.2,REQ=6378.160)
parameter (dr=0.001*req)
COMMON/DIPOLE/COLAT,ELON,VP,CTP,STP
real lShell, rStart, lonStart, ang_save, ang, cte, ste, tal
real stfcpa, stfspa, dif, dif_save, magpot
real xMag, yMag, zMag, bMag, GeoLat, GeoLon, GeoAlt
real aLatTest, aLonTest
integer i
CALL COFRM(DATE)
if (sLat < -90.0) then
lShell = 1.0/(cos(aLat*DTOR)**2) ! L value
rStart = lShell * Req
cte = 0.0
ste = sqrt(1.0-cte*cte)
tal = tan(alon*dtor)
ang_save = 1000.0
dif_save = 1000.0
mid = alon
if (mid < 0.0) mid = mid + 360.0
do i=mid*10.0-100,mid*10.0+100
ang = real(i)/10.0*dtor
stfcpa = ste*ctp*cos(ang)-cte*stp
stfspa = sin(ang)*ste
dif = abs(tal*stfcpa - stfspa)
if (dif < dif_save) then
dif_save = dif
ang_save = ang
endif
enddo
lonStart = ang_save + elon*dtor
if (lonStart < 0.0) lonStart = lonStart + 6.2831855
rotate = - cos(ang_save) * colat*dtor
ang = ang_save
stfcpa = ste*ctp*cos(ang)-cte*stp
stfspa = sin(ang)*ste
c write(*,*) 'xlon:',lonStart,ang_save, alon, cte, STFCPA, STFSPA, ctp
GeoAlt = rStart - Req
GeoLon = lonStart
GeoLat = 0.0
SGN = SIGN(1.,alat)
CALL GD2CART (GeoLAT,GeoLON,geoALT,XXX,YYY,ZZZ)
r = sqrt(xxx*xxx + yyy*yyy + zzz*zzz)
do while (r > req+Alt) ! Alt = AltMinIono
CALL FELDG(2,xxx,yyy,zzz,xMAG,yMAG,zMAG,bMag)
xmag = xmag/bmag
ymag = ymag/bmag
zmag = zmag/bmag
xxx = xxx + xmag * dr*sgn
yyy = yyy + ymag * dr*sgn
zzz = zzz + zmag * dr*sgn
r = sqrt(xxx*xxx + yyy*yyy + zzz*zzz)
c write(*,*) xxx,yyy,zzz,r, xmag,ymag,zmag
enddo
geolat = asin(zzz/r)
geolon = asin(xxx/sqrt(xxx*xxx+yyy*yyy))
if (yyy < 0.0) geolon = 6.2831855 - geolon
geoalt = r-req
gLat = GeoLat * rtod
gLon = GeoLon * rtod
aLatTest = 1000.0
aLonTest = 1000.0
gLatGuess = gLat
gLonGuess = gLon
iCount = 0
do while ((abs(aLatTest-aLat) > 0.01 .or. abs(aLonTest-aLon) > 0.1) .and.
! iCount < 20)
call APEX(DATE,gLatGuess,gLonGuess,Alt,lShell,aLatTest,aLonTest,
! bmag,xmag,ymag,zmag,MagPot)
c write(*,*) aLat, aLatTest, aLon, aLonTest, gLatGuess, gLonGuess, iCount
gLatGuess = gLatGuess + (aLat-aLatTest)/4.0
if (gLatGuess > 90.0) then
iCount = 20
gLatGuess = 88.0
gLonGuess = gLonGuess + 30.0
elseif (gLatGuess < -90.0) then
iCount = 20
gLatGuess = -88.0
gLonGuess = gLonGuess + 30.0
else
dLon = aLon-aLonTest
if (dLon > 300) dLon = dLon - 360.0
if (dLon < -300) dLon = dLon + 360.0
gLonGuess = gLonGuess + dLon/4.0
endif
if (gLonGuess > 360.0) gLonGuess = gLonGuess - 360.0
if (gLonGuess < 0.0) gLonGuess = gLonGuess + 360.0
iCount = iCount + 1
enddo
else
GeoLat = sLat
GeoLon = sLon
iCount = 20
gLat = sLat
endif
c ACTUAL
c write(*,*) "iCount : ",iCount, gLatGuess, gLonGuess, gLat, gLon
if (iCount >= 20) then
CALL GD2CART (aLat,aLon,ALT,aXXX,aYYY,aZZZ)
if (abs(gLatGuess) < 85.0) then
gLatGuess = gLat
else
gLat = gLatGuess
endif
iLonBest = -100
diff = 10000.0
do iLon = 0, 360, 10
gLonGuess = real(iLon)
call APEX(DATE,gLatGuess,gLonGuess,Alt,lShell,aLatTest,aLonTest,
! bmag,xmag,ymag,zmag,MagPot)
CALL GD2CART (aLatTest,aLonTest,ALT,tXXX,tYYY,tZZZ)
dist = sqrt((aXXX-tXXX)**2 + (aYYY-tYYY)**2 + (aZZZ-tZZZ)**2)
if (dist < diff) then
diff = dist
iLonBest = iLon
endif
enddo
gLon = real(iLonBest)
dLon = 20.0
LatFac = 110.0
iCount = 1
do while (diff > 1.0 .and. iCount < 20)
dLat = diff / LatFac
LatFac = LatFac * 0.99
do iLat = -5,5
gLatGuess = gLat + real(iLat)/10 * dLat
gLonGuess = gLon
call APEX(DATE,gLatGuess,gLonGuess,Alt,lShell,aLatTest,aLonTest,
! bmag,xmag,ymag,zmag,MagPot)
CALL GD2CART (aLatTest,aLonTest,ALT,tXXX,tYYY,tZZZ)
dist = sqrt((aXXX-tXXX)**2 + (aYYY-tYYY)**2 + (aZZZ-tZZZ)**2)
if (dist < diff) then
diff = dist
gLat = gLatGuess
if (gLat < -90.0) gLat = -180.0 - gLat
if (gLat > 90.0) gLat = 180.0 - gLat
c write(*,*) aLat, aLatTest, aLon, aLonTest, gLat, gLon, dLat, dLon, diff
endif
enddo
do iLon = -5,5
gLatGuess = gLat
gLonGuess = gLon + real(iLon)/10 * dLat
call APEX(DATE,gLatGuess,gLonGuess,Alt,lShell,aLatTest,aLonTest,
! bmag,xmag,ymag,zmag,MagPot)
CALL GD2CART (aLatTest,aLonTest,ALT,tXXX,tYYY,tZZZ)
dist = sqrt((aXXX-tXXX)**2 + (aYYY-tYYY)**2 + (aZZZ-tZZZ)**2)
if (dist < diff) then
diff = dist
gLon = gLonGuess
if (gLon < 0.0) gLon = gLon + 360.0
if (gLon > 360.0) gLon = gLon - 360.0
c write(*,*) aLat, aLatTest, aLon, aLonTest, gLat, gLon, dLat, dLon, diff
endif
enddo
dLon = dLon * exp((real(iCount)-5)/10)
iCount = iCount + 1
enddo
else
gLat = gLatGuess
gLon = gLonGuess
endif
if (gLon > 360) gLon = gLon - 360
if (gLon < 0) gLon = gLon + 360
end
C FILE NAME: apex.f
SUBROUTINE APEX (DATE,DLAT,DLON,ALT,
+ A,ALAT,ALON,BMAG,XMAG,YMAG,ZMAG,V)
C Calculate apex radius, latitude, longitude; and magnetic field and
C scalar magnetic potential.
C
C INPUTS:
C DATE = Year and fraction (1990.0 = 1990 January 1, 0 UT)
C DLAT = Geodetic latitude in degrees
C DLON = Geodetic longitude in degrees
C ALT = Altitude in km
C
C RETURNS:
C A = (Apex height + REQ)/REQ, where REQ = equatorial Earth radius.
C A is analogous to the L value in invariant coordinates.
C ALAT = Apex latitude in degrees (negative in S. magnetic hemisphere)
C ALON = Apex longitude (geomagnetic longitude of apex) in degrees
C BMAG = geomagnetic field magnitude (nT)
C XMAG = geomagnetic field component (nT): north
C YMAG = geomagnetic field component (nT): east
C ZMAG = geomagnetic field component (nT): downward
C V = geomagnetic potential (T-m)
C
C COMMON BLOCKS:
C COMMON /DIPOLE/ COLAT,ELON,VP,CTP,STP
C
C DIPOLE has IGRF variables obtained from routines in magfld.f:
C COLAT = Geocentric colatitude of geomagnetic dipole north pole (deg)
C ELON = East longitude of geomagnetic dipole north pole (deg)
C VP = Magnitude (T-m) of dipole component of magnetic potential at
C geomagnetic pole and geocentric radius of 6371.0088 km
C CTP = cosine of COLAT
C STP = sine of COLAT
C
C------------------------------------------------------------------------------
C HISTORY:
C Aug 1994: First version completed on the 22nd by A.D. Richmond.
C May 1999: Revise DS calculation in LINAPX to avoid divide by zero.
C Apr 2004: - Change definition of earth's equatorial radius (REQ)
C from the IAU-1966 spheroid (6378.160 km) to the WGS-1984
C spheroid (6378.137 km); see description below.
C - Revise comments toward a consistent format so they are
C easy to read.
C - Replace computed GO TO in ITRACE with IF blocks.
C - Refine FNDAPX to insure |Bdown/Btot| < 1.E-6 at apex
C Nov 2009: Change definition of earth's mean radius (RE) from 6371.2
C to the WGS84 value (6371.0088), by J.T. Emmert, NRL
C
C------------------------------------------------------------------------------
C Reference Spheroid Change March 2004
C
C Apex geomagnetic coordinates are based on the International Reference
C Geomagnetic Field (IGRF) which involves the earth's shape when converting
C geographic latitude and altitude to geocentric coordinates. For this
C purpose, the earth is assumed to be an ellipsoid which is fatter at the
C equator than the poles. The World Geodetic System 1984 spheroid
C (WGS-1984) is recommended in the recent release of IGRF-9 because it is
C used to position current satellite magnetic data (EOS Volume 84 Number 46
C November 18 2003). This differs from previous IGRF releases which favored
C the International Astronomical Union 1966 spheroid (IAU-1966) so the Apex
C program conversion from geographic to geocentric coordinates in subroutine
C CONVRT of file magfld.f has been revised from the IAU-1966 spheroid to the
C WGS-1984 spheroid.
C
C The spheroid used to prepare earlier IGRF releases is not always known but
C changing spheroids now produces differences at the earth's surface less
C than 1 nT, significantly less than other error sources: viz., 9 nT RMS
C error for older measurements reporting 1 nT accuracy, 200 nT in the
C vicinity of magnetized rocks, or as much as 1000 nT during and after a
C geomagnetic storm (www.ngdc.noaa.gov/IAGA/vmod/index.html).
C
C The elliptical shape is characterized in subroutine CONVRT by eccentricity
C (e) which is related to the the earth's equatorial radius (a) and polar
C radius (b) by
C
C e**2 = 1 - (b/a)**2 (1)
C
C This term is part of an eighth order Lagrange expansion formula (Astron.
C J. Vol. 66, p. 15-16, 1961) designed to give eight digit conversion
C accuracy. The following table summarizes the relevant spheroids:
C
C a b e**2 Source
C ----------- ----------- ----------- --------------
C - - 0.006722670 Astron J. 1961
C 6378.160 km 6356.775 km 0.006701642 IAU-1966
C 6378.137 km 6356.752 km 0.006694478 WGS-1984
C
C The previous formulation in CONVRT used the oblateness factor (epsilon),
C a surrogate for eccentricity, where
C
C e**2 = 2*epsilon - epsilon**2
C
C with oblateness revised from the 1961 paper's suggested 1/297 to 1/298.25
C in accordance with the IAU-1966 spheroid. Now CONVRT is reformulated to
C use equation 1 with the WGS-1984 spheroid's major axis (a) and minor axis
C (b) for which epsilon is approximately 1/298.2528.
C
C In addition to earth's equatorial radius and polar radius, the reference
C radius (Re) of 6371.2 km is explicit in the IGRF formula for magnetic
C potential and implicit in the derived magnetic field coefficients. The
C reference radius has not changed in IGRF releases.
C
C------------------------------------------------------------------------------
PARAMETER (RE = 6371.0088, DTOR = .01745329251994330)
COMMON /DIPOLE/ COLAT,ELON,VP,CTP,STP
CALL COFRM (DATE)
CALL DYPOL (CLATP,POLON,VPOL)
COLAT = CLATP
CTP = COS(CLATP*DTOR)
STP = SQRT(1. - CTP*CTP)
ELON = POLON
VP = VPOL
CALL LINAPX (DLAT,DLON,ALT, A,ALAT,ALON,XMAG,YMAG,ZMAG,BMAG)
XMAG = XMAG*1.E5 ! convert from gauss to nT
YMAG = YMAG*1.E5
ZMAG = ZMAG*1.E5
BMAG = BMAG*1.E5
CALL GD2CART (DLAT,DLON,ALT,X,Y,Z)
CALL FELDG (3, X/RE,Y/RE,Z/RE, BX,BY,BZ,V)
RETURN
END
SUBROUTINE LINAPX (GDLAT,GLON,ALT, A,ALAT,ALON,XMAG,YMAG,ZMAG,F)
C Transform geographic coordinates to Apex coordinates.
C
C INPUTS:
C GDLAT = Latitude (degrees, positive northward)
C GLON = Longitude (degrees, positive eastward)
C ALT = Height of starting point (km above mean sea level)
C
C OUTPUTS:
C A = (Apex height + REQ)/REQ, where REQ = equatorial Earth radius.
C A is analogous to the L value in invariant coordinates.
C ALAT = Apex Lat. (deg)
C ALON = Apex Lon. (deg)
C XMAG = Geomagnetic field component (gauss): north
C YMAG = Geomagnetic field component (gauss): east
C ZMAG = Geomagnetic field component (gauss): down
C F = Geomagnetic field magnitude (gauss)
C
C Trace the geomagnetic field line from the given location to find the
C apex of the field line. Before starting iterations to trace along
C the field line: (1) Establish a step size (DS, arc length in km)
C based on the geomagnetic dipole latitude; (2) determine the step
C direction from the sign of the vertical component of the geomagnetic
C field; and (3) convert to geocentric cartesian coordinates. Each
C iteration increments a step count (NSTP) and calls ITRACE to move
C along the the field line until reaching the iteration count limit
C (MAXS) or passing the apex (IAPX=2) and then calling FNDAPX to
C determine the apex location from the last three step locations
C (YAPX); however, if reaching the iteration limit, apex coordinates
C are calculated by DIPAPX which assumes a simplified dipole field.
C
C COMMON BLOCKS:
C COMMON /APXIN/ YAPX(3,3)
C COMMON /DIPOLE/ COLAT,ELON,VP,CTP,STP
C COMMON /FLDCOMD/ BX, BY, BZ, BB
C COMMON /ITRA/ NSTP, Y(3), YOLD(3), SGN, DS
C
C APXIN has step locations determined in ITRACE:
C YAPX = Matrix of cartesian coordinates (loaded columnwise) of the
C three points about the apex. Set in subroutine ITRACE.
C
C DIPOLE has IGRF variables obtained from routines in magfld.f:
C COLAT = Geocentric colatitude of geomagnetic dipole north pole (deg)
C ELON = East longitude of geomagnetic dipole north pole (deg)
C VP = Magnitude (T-m) of dipole component of magnetic potential at
C geomagnetic pole and geocentric radius of 6371.0088 km
C CTP = cosine of COLAT
C STP = sine of COLAT
C
C FLDCOMD has geomagnetic field at current trace point:
C BX = X component (Gauss)
C BY = Y component (Gauss)
C BZ = Z component (Gauss)
C BB = Magnitude (Gauss)
C
C ITRA has field line tracing variables determined in LINAPX:
C NSTP = Step count.
C Y = Array containing current tracing point cartesian coordinates.
C YOLD = Array containing previous tracing point cartesian coordinates.
C SGN = Determines direction of trace.
C DS = Step size (arc length in km).
C
C REFERENCES:
C Stassinopoulos E. G. , Mead Gilbert D., X-841-72-17 (1971) GSFC,
C Greenbelt, Maryland
C
C EXTERNALS:
C GD2CART = Convert geodetic to geocentric cartesian coordinates (in magfld.f)
C CONVRT = Convert geodetic to geocentric cylindrical or geocentric spherical
C and back (in magfld.f).
C FELDG = Obtain IGRF magnetic field components (in magfld.f).
C ITRACE = Follow a geomagnetic field line
C DIPAPX = Compute apex coordinates assuming a geomagnetic dipole field
C FNDAPX = Compute apex coordinates from the last three traced field line points
C
C------------------------------------------------------------------------------
C HISTORY:
C Oct 1973: Initial version completed on the 29th by Wally Clark, NOAA
C ERL Lab.
C Feb 1988: Revised on the 1st by Harsh Anand Passi, NCAR.
C Aug 1994: Revision by A. D. Richmond, NCAR.
C Nov 2009: Change definition of earth's mean radius (RE) from 6371.2
C to the WGS84 value (6371.0088), by J.T. Emmert, NRL.
PARAMETER (MAXS = 200, RTOD = 57.2957795130823, RE =6371.0088,
+ DTOR = .01745329251994330, REQ=6378.137)
COMMON /FLDCOMD/ BX, BY, BZ, BB
COMMON /APXIN/ YAPX(3,3)
COMMON /DIPOLE/ COLAT,ELON,VP,CTP,STP
COMMON /ITRA/ NSTP, Y(3), YP(3), SGN, DS
C Set step size based on the geomagnetic dipole latitude of the starting point
CALL CONVRT (2,GDLAT,ALT,GCLAT,R)
SINGML = CTP*SIN(GCLAT*DTOR) + STP*COS(GCLAT*DTOR)*
+ COS((GLON-ELON)*DTOR)
C May 1999: avoid possible divide by zero (when SINGML = 1.): the old version
C limited DS to its value at 60 deg GM latitude with: DS = .06*R/(1.-SINGML*SINGML) - 370.
C IF (DS .GT. 1186.) DS = 1186.
CGML2 = AMAX1 (0.25,1.-SINGML*SINGML)
DS = .06*R/CGML2 - 370.
C Initialize YAPX array
DO 4 J=1,3
DO 4 I=1,3
4 YAPX(I,J) = 0.
C Convert from geodetic to earth centered cartesian coordinates
CALL GD2CART (GDLAT,GLON,ALT,Y(1),Y(2),Y(3))
NSTP = 0
C Get magnetic field components to determine the direction for
C tracing the field line
CALL FELDG (1,GDLAT,GLON,ALT,XMAG,YMAG,ZMAG,F)
SGN = SIGN (1.,-ZMAG)
C Use cartesian coordinates to get magnetic field components
C (from which gradients steer the tracing)
10 CALL FELDG (2, Y(1)/RE,Y(2)/RE,Y(3)/RE, BX,BY,BZ,BB)
NSTP = NSTP + 1
IF (NSTP .LT. MAXS) THEN
CALL ITRACE (IAPX) ! trace along field line
IF (IAPX .EQ. 1) GO TO 10
CALL FNDAPX (ALT,ZMAG,A,ALAT,ALON) ! (IAPX=2) => passed max radius; find its coordinates
ELSE
RHO = SQRT (Y(1)*Y(1) + Y(2)*Y(2)) ! too many steps; get apex from dipole approximation
CALL CONVRT (3,XLAT,HT,RHO,Y(3))
XLON = RTOD*ATAN2 (Y(2),Y(1))
CALL FELDG (1,XLAT,XLON,HT,BNRTH,BEAST,BDOWN,BABS)
CALL DIPAPX (XLAT,XLON,HT,BNRTH,BEAST,BDOWN,A,ALON)
ALAT = -SGN*RTOD*ACOS (SQRT(1./A))
ENDIF
RETURN
END
SUBROUTINE ITRACE (IAPX)
C Follow a geomagnetic field line until passing its apex
C
C INPUTS:
C (all are in common blocks)
C OUTPUTS:
C IAPX = 2 (when apex passed) or 1 (not)
C
C This uses the 4-point Adams formula after initialization.
C First 7 iterations advance point by 3 steps.
C
C COMMON BLOCKS:
C COMMON /APXIN/ YAPX(3,3)
C COMMON /FLDCOMD/ BX, BY, BZ, BB
C COMMON /ITRA/ NSTP, Y(3), YOLD(3), SGN, DS
C
C APXIN has step locations determined in ITRACE:
C YAPX = Matrix of cartesian coordinates (loaded columnwise) of the
C three points about the apex. Set in subroutine ITRACE.
C
C FLDCOMD has geomagnetic field at current trace point:
C BX = X component (Gauss)
C BY = Y component (Gauss)
C BZ = Z component (Gauss)
C BB = Magnitude (Gauss)
C
C ITRA has field line tracing variables determined in LINAPX:
C NSTP = Step count.
C Y = Array containing current tracing point cartesian coordinates.
C YOLD = Array containing previous tracing point cartesian coordinates.
C SGN = Determines direction of trace.
C DS = Step size (arc length in km).
C
C REFERENCES:
C Stassinopoulos E. G. , Mead Gilbert D., X-841-72-17 (1971) GSFC,
C Greenbelt, Maryland
C------------------------------------------------------------------------------
C HISTORY:
C Oct 1973: Initial version completed on the 29th by W. Clark, NOAA ERL
C Laboratory.
C Feb 1988: Revised by H. Passi, NCAR.
C Apr 2004: Replace computed GO TO with IF blocks because some compilers
C are threatening to remove this old feature
C
COMMON /APXIN/ YAPX(3,3)
COMMON /FLDCOMD/ BX, BY, BZ, BB
COMMON /ITRA/ NSTP, Y(3), YOLD(3), SGN, DS
DIMENSION YP(3,4)
SAVE
C Statement function
RDUS(D,E,F) = SQRT (D**2 + E**2 + F**2)
IAPX = 1
C Cartesian component magnetic field (partial) derivitives steer the trace
YP(1,4) = SGN*BX/BB
YP(2,4) = SGN*BY/BB
YP(3,4) = SGN*BZ/BB
IF (NSTP .LE. 7) THEN
DO 10 I=1,3
IF (NSTP .EQ. 1) THEN
D2 = DS/2.
D6 = DS/6.
D12 = DS/12.
D24 = DS/24.
YP(I,1) = YP(I,4)
YOLD(I) = Y(I)
YAPX(I,1) = Y(I)
Y(I) = YOLD(I) + DS*YP(I,1)
ELSE IF (NSTP .EQ. 2) THEN
YP(I,2) = YP(I,4)
Y(I) = YOLD(I) + D2*(YP(I,2)+YP(I,1))
ELSE IF (NSTP .EQ. 3) THEN
Y(I) = YOLD(I) + D6*(2.*YP(I,4)+YP(I,2)+3.*YP(I,1))
ELSE IF (NSTP .EQ. 4) THEN
YP(I,2) = YP(I,4)
YAPX(I,2) = Y(I)
YOLD(I) = Y(I)
Y(I) = YOLD(I) + D2*(3.*YP(I,2)-YP(I,1))
ELSE IF (NSTP .EQ. 5) THEN
Y(I) = YOLD(I) + D12*(5.*YP(I,4)+8.*YP(I,2)-YP(I,1))
ELSE IF (NSTP .EQ. 6) THEN
YP(I,3) = YP(I,4)
YOLD(I) = Y(I)
YAPX(I,3) = Y(I)
Y(I) = YOLD(I) + D12*(23.*YP(I,3)-16.*YP(I,2)+5.*YP(I,1))
ELSE IF (NSTP .EQ. 7) THEN
YAPX(I,1) = YAPX(I, 2)
YAPX(I,2) = YAPX(I, 3)
Y(I) = YOLD(I) + D24*(9.*YP(I,4) + 19.*YP(I,3) -
+ 5.*YP(I,2) + YP(I,1))
YAPX(I,3) = Y(I)
ENDIF
10 CONTINUE
IF (NSTP .EQ. 6 .OR. NSTP .EQ. 7) THEN ! signal if apex passed
RC = RDUS (YAPX(1,3), YAPX(2,3), YAPX(3,3))
RP = RDUS (YAPX(1,2), YAPX(2,2), YAPX(3,2))
IF (RC .LT. RP) IAPX = 2
ENDIF
ELSE ! NSTP > 7
DO 30 I=1,3
YAPX(I,1) = YAPX(I,2)
YAPX(I,2) = Y(I)
YOLD(I) = Y(I)
Y(I) = YOLD(I) + D24*(55.*YP(I,4) - 59.*YP(I,3) +
+ 37.*YP(I,2) - 9.*YP(I,1))
YAPX(I,3) = Y(I)
DO 20 J=1,3
20 YP(I,J) = YP(I,J+1)
30 CONTINUE
RC = RDUS ( Y(1), Y(2), Y(3))
RP = RDUS (YOLD(1), YOLD(2), YOLD(3))
IF (RC .LT. RP) IAPX = 2
ENDIF
RETURN
END
SUBROUTINE FNDAPX (ALT,ZMAG,A,ALAT,ALON)
C Find apex coordinates once tracing (in subroutine ITRACE) has
C signalled that the apex has been passed.
C INPUTS:
C ALT = Altitude of starting point
C ZMAG = Downward component of geomagnetic field at starting point
C OUTPUT
C A = Apex radius, defined as (Apex height + Req)/Req, where
C Req = equatorial Earth radius.
C A is analogous to the L value in invariant coordinates.
C ALAT = Apex Lat. (deg)
C ALON = Apex Lon. (deg)
C
C COMMON BLOCKS:
C COMMON /APXIN/ YAPX(3,3)
C COMMON /DIPOLE/ COLAT,ELON,VP,CTP,STP
C
C APXIN has step locations determined in ITRACE:
C YAPX = Matrix of cartesian coordinates (loaded columnwise) of the
C three points about the apex. Set in subroutine ITRACE.
C
C DIPOLE has IGRF variables obtained from routines in magfld.f:
C COLAT = Geocentric colatitude of geomagnetic dipole north pole (deg)
C ELON = East longitude of geomagnetic dipole north pole (deg)
C VP = Magnitude (T-m) of dipole component of magnetic potential at
C geomagnetic pole and geocentric radius of 6371.0088 km
C CTP = cosine of COLAT
C STP = sine of COLAT
C
C EXTERNALS:
C FINT = Second degree interpolation routine
C------------------------------------------------------------------------------
C HISTORY:
C Oct 1973: Initial version completed on the 23rd by Clark, W., NOAA
C Boulder.
C Aug 1994: Revision on the 3rd by A.D. Richmond, NCAR
C Apr 2004: Repair problem noted by Dan Weimer where the apex location
C produced by FINT may still have a non-zero vertical magnetic
C field component.
PARAMETER (RTOD = 57.2957795130823,
+ DTOR = .01745329251994330, REQ=6378.137)
COMMON /APXIN/ YAPX(3,3)
COMMON /DIPOLE/ COLAT,ELON,VP,CTP,STP
DIMENSION BD(3), Y(3)
C Get geodetic height and vertical (downward) component of the magnetic
C field at last three points found by ITRACE
DO 10 I=1,3
RHO = SQRT (YAPX(1,I)**2 + YAPX(2,I)**2)
CALL CONVRT (3,GDLT,HT, RHO,YAPX(3,I))
GDLN = RTOD*ATAN2 (YAPX(2,I),YAPX(1,I))
10 CALL FELDG (1,GDLT,GDLN,HT, BN,BE,BD(I),BMAG)
C Interpolate to where Bdown=0 to find cartesian coordinates at dip equator
NITR = 0
20 Y(1) = FINT (BD(1),BD(2),BD(3),YAPX(1,1),YAPX(1,2),YAPX(1,3), 0.)
Y(2) = FINT (BD(1),BD(2),BD(3),YAPX(2,1),YAPX(2,2),YAPX(2,3), 0.)
Y(3) = FINT (BD(1),BD(2),BD(3),YAPX(3,1),YAPX(3,2),YAPX(3,3), 0.)
C Insure negligible Bdown or
C
C |Bdown/Btot| < 2.E-6
C
C For instance, Bdown must be less than 0.1 nT at low altitudes where
C Btot ~ 50000 nT. This ratio can be exceeded when interpolation is
C not accurate; i.e., when the middle of the three points interpolated
C is too far from the dip equator. The three points were initially
C defined with equal spacing by ITRACE, so replacing point 2 with the
C most recently fit location will reduce the interpolation span.
RHO = SQRT (Y(1)**2 + Y(2)**2)
GDLN = RTOD*ATAN2 (Y(2),Y(1))
CALL CONVRT (3,GDLT,HTA, RHO,Y(3))
CALL FELDG (1,GDLT,GDLN,HTA, BNA,BEA,BDA,BA)
ABDOB = ABS(BDA/BA)
IF (ABDOB .GT. 2.E-6) THEN
IF (NITR .LT. 4) THEN ! 4 was chosen because tests rarely required 2 iterations
NITR = NITR + 1
YAPX(1,2) = Y(1)
YAPX(2,2) = Y(2)
YAPX(3,2) = Y(3)
BD(2) = BDA
GO TO 20
ELSE
WRITE (0,'(''APEX: Imprecise fit of apex: |Bdown/B| ='',1PE7.1
+ )') ABDOB
ENDIF
ENDIF
C Ensure altitude of the Apex is at least the initial altitude when
C defining the Apex radius then use it to define the Apex latitude whose
C hemisphere (sign) is inferred from the sign of the dip angle at the
C starting point
A = (REQ + AMAX1(ALT,HTA)) / REQ
IF (A .LT. 1.) THEN
WRITE (0,'(''APEX: A can not be less than 1; A, REQ, HTA: '',1P3
+E15.7)') A,REQ,HTA
CALL EXIT (1)
ENDIF
RASQ = ACOS (SQRT(1./A))*RTOD
ALAT = SIGN (RASQ,ZMAG)
C ALON is the dipole longitude of the apex and is defined using
C spherical coordinates where
C GP = geographic pole.
C GM = geomagnetic pole (colatitude COLAT, east longitude ELON).
C XLON = longitude of apex.
C TE = colatitude of apex.
C ANG = longitude angle from GM to apex.
C TP = colatitude of GM.
C TF = arc length between GM and apex.
C PA = ALON be geomagnetic longitude, i.e., Pi minus angle measured
C counterclockwise from arc GM-apex to arc GM-GP.
C then, spherical-trigonometry formulas for the functions of the angles
C are as shown below. Notation uses C=cos, S=sin and STFCPA = sin(TF) * cos(PA),
C STFSPA = sin(TF) * sin(PA)
XLON = ATAN2 (Y(2),Y(1))
ANG = XLON-ELON*DTOR
CANG = COS (ANG)
SANG = SIN (ANG)
R = SQRT (Y(1)**2 + Y(2)**2 + Y(3)**2)
CTE = Y(3)/R
STE = SQRT (1.-CTE*CTE)
STFCPA = STE*CTP*CANG - CTE*STP
STFSPA = SANG*STE
ALON = ATAN2 (STFSPA,STFCPA)*RTOD
RETURN
END
SUBROUTINE DIPAPX (GDLAT,GDLON,ALT,BNORTH,BEAST,BDOWN, A,ALON)
C Compute A, ALON from local magnetic field using dipole and spherical
C approximation.
C
C INPUTS:
C GDLAT = geodetic latitude, degrees
C GDLON = geodetic longitude, degrees
C ALT = altitude, km
C BNORTH = geodetic northward magnetic field component (any units)
C BEAST = eastward magnetic field component
C BDOWN = geodetic downward magnetic field component
C OUTPUTS:
C A = apex radius, 1 + h_A/R_eq
C ALON = apex longitude, degrees
C
C Use spherical coordinates and define:
C GP = geographic pole.
C GM = geomagnetic pole (colatitude COLAT, east longitude ELON).
C G = point at GDLAT,GDLON.
C E = point on sphere below apex of dipolar field line passing
C through G.
C TD = dipole colatitude of point G, found by applying dipole
C formula for dip angle to actual dip angle.
C B = Pi plus local declination angle. B is in the direction
C from G to E.
C TG = colatitude of G.
C ANG = longitude angle from GM to G.
C TE = colatitude of E.
C TP = colatitude of GM.
C A = longitude angle from G to E.
C APANG = A + ANG
C PA = geomagnetic longitude, i.e., Pi minus angle measured
C counterclockwise from arc GM-E to arc GM-GP.
C TF = arc length between GM and E.
C Then, using notation C=cos, S=sin, COT=cot, spherical-trigonometry
C formulas for the functions of the angles are as shown below. Note:
C STFCPA = sin(TF) * cos(PA)
C STFSPA = sin(TF) * sin(PA)
C
C COMMON BLOCKS:
C COMMON /DIPOLE/ COLAT,ELON,VP,CTP,STP
C
C DIPOLE has IGRF variables obtained from routines in magfld.f:
C COLAT = Geocentric colatitude of geomagnetic dipole north pole (deg)
C ELON = East longitude of geomagnetic dipole north pole (deg)
C VP = Magnitude (T-m) of dipole component of magnetic potential at
C geomagnetic pole and geocentric radius of 6371.0088 km
C CTP = cosine of COLAT
C STP = sine of COLAT
C------------------------------------------------------------------------------
C HISTORY:
C May 1994: Completed on the 1st by A. D. Richmond
C Nov 2009: Change definition of earth's mean radius (RE) from 6371.2
C to the WGS84 value (6371.0088), by J.T. Emmert, NRL.
PARAMETER (RTOD = 57.2957795130823, RE =6371.0088,
+ DTOR = .01745329251994330, REQ=6378.137)
COMMON /DIPOLE/ COLAT,ELON,VP,CTP,STP
BHOR = SQRT(BNORTH*BNORTH + BEAST*BEAST)
IF (BHOR .EQ. 0.) THEN
ALON = 0.
A = 1.E34
RETURN
ENDIF
COTTD = BDOWN*.5/BHOR
STD = 1./SQRT(1.+COTTD*COTTD)
CTD = COTTD*STD
SB = -BEAST /BHOR
CB = -BNORTH/BHOR
CTG = SIN (GDLAT*DTOR)
STG = COS (GDLAT*DTOR)
ANG = (GDLON-ELON)*DTOR
SANG = SIN(ANG)
CANG = COS(ANG)
CTE = CTG*STD + STG*CTD*CB
STE = SQRT(1. - CTE*CTE)
SA = SB*CTD/STE
CA = (STD*STG - CTD*CTG*CB)/STE
CAPANG = CA*CANG - SA*SANG
SAPANG = CA*SANG + SA*CANG
STFCPA = STE*CTP*CAPANG - CTE*STP
STFSPA = SAPANG*STE
ALON = ATAN2 (STFSPA,STFCPA)*RTOD
R = ALT + RE
HA = ALT + R*COTTD*COTTD
A = 1. + HA/REQ
RETURN
END
FUNCTION FINT (X1,X2,X3,Y1,Y2,Y3, XFIT)
C Second degree interpolation used by FNDAPX
C INPUTS:
C X1 = point 1 ordinate value
C X2 = point 2 ordinate value
C X3 = point 3 ordinate value
C Y1 = point 1 abscissa value
C Y2 = point 2 abscissa value
C Y3 = point 3 abscissa value
C XFIT = ordinate value to fit
C RETURNS:
C YFIT = abscissa value corresponding to XFIT
C
C MODIFICATIONS:
C Apr 2004: Change from subroutine to function, rename variables and
C add intermediates which are otherwise calculated twice
X12 = X1-X2
X13 = X1-X3
X23 = X2-X3
XF1 = XFIT-X1
XF2 = XFIT-X2
XF3 = XFIT-X3
FINT = (Y1*X23*XF2*XF3 - Y2*X13*XF1*XF3 + Y3*X12*XF1*XF2) /
+ (X12*X13*X23)
RETURN
END
SUBROUTINE COFRM (DATE)
C Define the International Geomagnetic Reference Field (IGRF) as a
C scalar potential field using a truncated series expansion with
C Schmidt semi-normalized associated Legendre functions of degree n and
C order m. The polynomial coefficients are a function of time and are
C interpolated between five year epochs or extrapolated at a constant
C rate after the last epoch.
C
C INPUTS:
C DATE = yyyy.fraction (UT)
C OUTPUTS (in comnon block MAGCOF):
C NMAX = Maximum order of spherical harmonic coefficients used
C GB = Coefficients for magnetic field calculation
C GV = Coefficients for magnetic potential calculation
C ICHG = Flag indicating when GB,GV have been changed in COFRM
C
C It is fatal to supply a DATE before the first epoch. A warning is
C issued to Fortran unit 0 (stderr) if DATE is later than the
C recommended limit, five years after the last epoch.
C
C HISTORY (blame):
C Apr 1983: Written by Vincent B. Wickwar (Utah State Univ.) including
C secular variation acceleration rate set to zero in case the IGRF
C definition includes such second time derivitives. The maximum degree
C (n) defined was 10.
C
C Jun 1986: Updated coefficients adding Definitive Geomagnetic Reference
C Field (DGRF) for 1980 and IGRF for 1985 (EOS Volume 7 Number 24). The
C designation DGRF means coefficients will not change in the future
C whereas IGRF coefficients are interim pending incorporation of new
C magnetometer data. Common block MAG was replaced by MAGCOF, thus
C removing variables not used in subroutine FELDG. (Roy Barnes)
C
C Apr 1992 (Barnes): Added DGRF 1985 and IGRF 1990 as given in EOS
C Volume 73 Number 16 April 21 1992. Other changes were made so future
C updates should: (1) Increment NDGY; (2) Append to EPOCH the next IGRF
C year; (3) Append the next DGRF coefficients to G1DIM and H1DIM; and (4)
C replace the IGRF initial values (G0, GT) and rates of change indices
C (H0, HT).
C
C Apr 1994 (Art Richmond): Computation of GV added, for finding magnetic
C potential.
C
C Aug 1995 (Barnes): Added DGRF for 1990 and IGRF for 1995, which were
C obtained by anonymous ftp to geomag.gsfc.nasa.gov (cd pub, mget table*)
C as per instructions from Bob Langel (langel@geomag.gsfc.nasa.gov) with
C problems reported to baldwin@geomag.gsfc.nasa.gov.
C
C Oct 1995 (Barnes): Correct error in IGRF-95 G 7 6 and H 8 7 (see email
C in folder). Also found bug whereby coefficients were not being updated
C in FELDG when IENTY did not change so ICHG was added to flag date
C changes. Also, a vestigial switch (IS) was removed from COFRM; it was
C always zero and involved 3 branch if statements in the main polynomial
C construction loop now numbered 200.
C
C Feb 1999 (Barnes): Explicitly initialize GV(1) in COFRM to avoid the
C possibility of compiler or loader options initializing memory to
C something else (e.g., indefinite). Also simplify the algebra in COFRM
C with no effect on results.
C
C Mar 1999 (Barnes): Removed three branch if's from FELDG and changed
C statement labels to ascending order.
C
C Jun 1999 (Barnes): Corrected RTOD definition in GD2CART.
C
C May 2000 (Barnes): Replace IGRF 1995, add IGRF 2000, and extend the
C earlier DGRF's back to 1900. The coefficients came from an NGDC web
C page. Related documentation is in $APXROOT/docs/igrf.2000.* where
C $APXROOT, defined by 'source envapex', is traditionally ~bozo/apex).
C
C Mar 2004 (Barnes): Replace 1995 and 2000 coefficients; now both are
C DGRF. Coefficients for 2000 are degree 13 with precision increased to
C tenths nT and accommodating this instigated changes: (1) degree (NMAX)
C is now a function of epoch (NMXE) to curtail irrelevant looping over
C unused high order terms (n > 10 in epochs before 2000) when calculating
C GB; (2) expand coefficients data statement layout for G1D and H1D,
C formerly G1DIM and H1DIM; (3) omit secular variation acceleration terms
C which were always zero; (4) increase array dimensions in common block
C MAGCOF and associated arrays G and H in FELDG; (5) change earth's shape
C in CONVRT from the IAU-1966 to the WGS-1984 spheroid; (6) eliminate
C reference to 'definitive' in variables in COFRM which were not always
C definitive; (7) change G to GB in COFRM s.t. arrays GB and GV in common
C block MAGCOF are consistently named in all subroutines; (8) remove
C unused constants in all five subroutines. See EOS Volume 84 Number 46
C November 18 2003, www.ngdc.noaa.gov/IAGA/vmod/igrf.html or local files
C $APXROOT/docs/igrf.2004.*
C
C Sept. 2005 (Maute): update with IGRF10 from
C http://www.ngdc.noaa.gov/IAGA/vmod/igrf.html use script
C ~maute/apex.d/apex_update/igrf2f Note that the downloaded file the start
C column of the year in the first line has to be before the start of each
C number in the same column
C
C Jan. 2010 (Maute) update with IGRF11 (same instructions as Sep. 2005
C comment