milmag.f
80 lines
| 2.4 KiB
| text/x-fortran
|
FortranFixedLexer
| r0 | C $Id: milmag.f 7206 2020-09-09 15:33:40Z brideout $ | |||
| C | ||||
| SUBROUTINE MILMAG(TM,RKM,ST,CT,SPH,CPH,BR,BT,BP,B) | ||||
| C | ||||
| C MILMAG evaluates the geomagnetic field at a point specified by | ||||
| C its geocentric coordinates. | ||||
| C | ||||
| C Modified by B. Rideout - Sep 9, 2020 | ||||
| C This method is now simply a thin wrapper around igrf13.f, | ||||
| C method igrf13syn. See igrf13.f for details | ||||
| C | ||||
| C Modified by B. Rideout - Aug. 31, 2005 | ||||
| C This method is now simply a thin wrapper around Geopack-2005 code, | ||||
| C method igrf_geo. See Geopack-2005.f for details | ||||
| C | ||||
| C Modified by B. Rideout - Dec. 26, 2002 | ||||
| C This method is now simply a thin wrapper around geo-cgm code, | ||||
| C method igrf. See geo-cgm.f for details | ||||
| C | ||||
| C Input: | ||||
| C TM - time in years for desired field (e.g. 1971.25) | ||||
| C RKM - geocentric distance (km) | ||||
| C ST,CT - sin and cos of geocentric colatitude | ||||
| C SPH,CPH - sin and cos of east longitude | ||||
| C | ||||
| C Output: | ||||
| C BR,BT,BP - geocentric field components (gauss) | ||||
| C B - magnitude of field (gauss) | ||||
| C .. | ||||
| C .. Scalar Arguments .. | ||||
| DOUBLE PRECISION B,BP,BR,BT,CPH,CT,RKM,SPH,ST,TM | ||||
| C .. | ||||
| C .. Local Scalars .. | ||||
| DOUBLE PRECISION GCLAT,GLON | ||||
| DOUBLE PRECISION R,T,F,RBR,RBT,RBP,RB | ||||
| INTEGER IYEAR | ||||
| C .. | ||||
| C .. Intrinsic Functions .. | ||||
| INTRINSIC ASIN,ANINT,SQRT | ||||
| C .. external functions | ||||
| EXTERNAL igrf13syn | ||||
| C convert from sin/cos to GCLAT, GLON | ||||
| IF (ST .GE. 0.0 .AND. CT .GE. 0.0) THEN | ||||
| GCLAT = ASIN(ST)*57.2958 | ||||
| ELSE IF (ST .GE. 0.0 .AND. CT .LE. 0.0) THEN | ||||
| GCLAT = 180 - ASIN(ST)*57.2958 | ||||
| ELSE IF (ST .LE. 0.0 .AND. CT .LE. 0.0) THEN | ||||
| GCLAT = 180 - ASIN(ST)*57.2958 | ||||
| ELSE | ||||
| GCLAT = ASIN(ST)*57.2958 + 360.0 | ||||
| END IF | ||||
| C leave GCLAT as degrees | ||||
| IF (SPH .GE. 0.0 .AND. CPH .GE. 0.0) THEN | ||||
| GLON = ASIN(SPH)*57.2958 | ||||
| ELSE IF (SPH .GE. 0.0 .AND. CPH .LE. 0.0) THEN | ||||
| GLON = 180 - ASIN(SPH)*57.2958 | ||||
| ELSE IF (SPH .LE. 0.0 .AND. CPH .LE. 0.0) THEN | ||||
| GLON = 180 - ASIN(SPH)*57.2958 | ||||
| ELSE | ||||
| GLON = ASIN(SPH)*57.2958 + 360.0 | ||||
| END IF | ||||
| C leave GLON as degrees | ||||
| CALL igrf13syn(0,TM,2,RKM,GCLAT,GLON,RBT,RBP,RBR,RB) | ||||
| C Set outputs as double percision | ||||
| BR=DBLE(RBR)*(-1E-5) | ||||
| BT=DBLE(RBT)*(-1E-5) | ||||
| BP=DBLE(RBP)*1E-5 | ||||
| B = SQRT(BT*BT + BP*BP + BR*BR) | ||||
| RETURN | ||||
| C | ||||
| END | ||||
