Computer program for transforming spherical harmonic coefficients to geomagnetic dipole coordinate

Mead, G. D.

English

NASA Technical Reports Server

General Disclaimer 
One or more of the Following Statements may affect this Document 
 
 This document has been reproduced from the best copy furnished by the 
organizational source. It is being released in the interest of making available as 
much information as possible. 
 
 This document may contain data, which exceeds the sheet parameters. It was 
furnished in this condition by the organizational source and is the best copy 
available. 
 
 This document may contain tone-on-tone or color graphs, charts and/or pictures, 
which have been reproduced in black and white. 
 
 This document is paginated as submitted by the original source. 
 
 Portions of this document are not fully legible due to the historical nature of some 
of the material. However, it is the best reproduction available from the original 
submission. 
 
 
 
 
 
 
 
Produced by the NASA Center for Aerospace Information (CASI) 
https://ntrs.nasa.gov/search.jsp?R=19700019385 2020-03-23T19:35:09+00:00Z
1I
--^	 x-moo-use
NASA TM F63121
THE INTERNATIONAL GEOMAGNETIC
REFERENCE FIELD 1965.0
IN DIPOLE COORDINATES
GILBERT D. MEAD-
MAY 1970
_`	 - Fx
	 .1J	 ..18-.	 T
GODDARD SPACE FLIGHT CENTER
GREENBELT, -
 MARYLAND
	
N70 
	 5"	 t	 U1CCE8flON NUMB[A!
	 -[	 ITNRU)O t
J ICODEI63cta^ 	 13
/NABA CR OR TM% OA AD NUMf[RI
ICAT[60RT)
THE INTERNATIONAL GEOMAGNETIC REFERENCE FIELD 1965.0
IN DIPOLE COORDINATES
Gilbert D. Mead
Laboratory for Space Physics
NASA Goddard Space Flight Center
Greenbelt, Md. 20771
At the working group session of the meeting on Quantitative :Iagnetcspheric
Models held in Boulder March 18-20, 1970, the need was recognized for a ref-
erence model of the earth's internal field with input and output in geomagnetic
dipole coordinates (Chapman, 1963) instead of geographic coordinates. Such a
model could be used much more readily in conjunction with various models of
external field sources (e.g., ring current, magnetopause, or tail fields)
whose orientation depends upon the direction of the magnetic dipole axis
rather than the rotation axis. 'Model field and field-line-tracing calculations
could then be performed in a physically more meaningful coordinate system whose
z-axis is aligned with the magnetic dipole; simple longitudinal transformations
would suffice to link this system with the solar magnetic coordinate system,
which is also aligned with the dipole. Deviations from a dipole field would
be much more readily identifiable.
A program has been written to transform the spherical harmonic coefficients
I
of any model field into any arbitrarily-tilted coordinate system, using a method
similar to Stern's (1965). For each value of n in the spherical harmonic
expansion, the scalar potential V was calculated at 2 n + 1 randomly-selected
points and set equal to the potential in the tilted system. The resulting
equations were solved for the 2 n + 1 unknown coefficients. The process was
then repeated for the time derivatives of the coefficients representing the
secular change.
* To be published (without Fortran listings) in the Journal of Geophysical
Research, Vol. 75, 1970.
2The resulting coefficients for the International Geomagnetic Reference
Field 190'5.0 ( I.AGA Commission 2 Working Group 4, 1969) are given in Table 1.
The normalization of the Legendre polynomials and all definitions are the
IAGA
same as in the^ reference except for the orientation of the coordinate system.
The geographic coordinates of the new pole are:
eo = 11.4350 geocentric colatitude
X0 = -69.7610 east longitude
as determined from the three n = 1 coefficients of the IGRF field at
epoch 1965.0. The input to a subroutine using these coefficients is
geomagnetic dipole colatitude (6d) with respect to this pole and east
geomagnetic dipole longitude (Ad) with respect to a meridian passing through
the north and south dipole poles and the south geographic pole. The XX
(northward) and YY(eastward) components of the field intensity obtained
as output are then measured with respect to the north dipole pole rather than
the geographic pole. The Z (downward) component is independent of the
orientation of the coordinate system.
It is often necessary to transform positions and/or field components
from one system to the other. The geometry is shown in Figure 1. P is
the geographic pole, D is the dipole pole, G is Greenwich, and A is a
point on the reference sphere (or the intersection of a vector through A
f
upon the celestial sphere) at which the field is to be evaluated; 8 g
 and
X  are the geographic colatitude and longitude of A; and d is the rotational
c
angle at A used to transform the x and Y field components; 6 is positive
if D appears to be west of P, negative if east. In a pure dipole field,
6 would be the negative of the declination angle everywhere.
3I
110° -Xd
00	 P
DIPOLE
MERIDIAN	 D	 -X0 X
9
A d	 09
GREENWICH 0d
	 SMERIDIAN
A
G
/--- DIPOLE
EQUATOR
GEOGRAPHIC
EQUATOR
Figure 1. Earth's surface (or, alternatively, surface of celestial
sphere), showing spherical triangle DPA needed for coordinate transformations.
The dipole tilt is exaggerated for greater clarity.
4From the solution of the spherical triangle EPA,
cosed = coseo coseg + sineo sine  cos(Xg-X0)
sinAd = + 1 - cos29-
cosXd = (cos8o cosed - coseg) / (sine0 sine 
r
sinAd = sineg
cos8
	 (coseo
sins = sineo
sin(ag-X0) / sine 
- cosed coseg) / (sinAd sine 
3in(Xg-Xo) / sine 
and
X  = Xg cos8 - Yg sins
Yd
 = X  sins + Y  cos8
where the sine and cosine of each angle determine the proper quadrant.
The inverse transformations are given by
coseg = coseo cosed - sineo sined cosxd
cos(?,g-ao) = (cosed - cose0 coseg) / (sineo sine g)
sin(Xg-Xo) = sine  sinAd / sineg
with similar expressions for the rotation of the Xd and Yd field components. 	
9
The dipole latitude will always differ from geographic latitude by less
than 11.50 . Except near the poles, dipole longitude will generally be
about 700 greater than geographic longitude and the angle 8 will usually
be less than ::15 0 . In a pure dipole field, Yd = 0 everywhere, and in the
IGRF field, Yd
 is small and apprc^ahes zero rapidly at large distances.
5The coefficients were checked by calculating the field components
at a number of points and at several time periods, both in geographic
coordinates, using the original coefficients, and in the corresponding
dipole coordinates, using the present coefficients. The calculated
intensities agreed with each other to within about 4 gammas at the
earth's surface and 0.2 gammas at a geocentric distance of 3 RE.
This is consistent with the fact that each of the transformed coeffi-
cients was rounded to the nearest gamma (derivatives to the nearest
tenth gamma/yr). Field directions corresponded to within 0.005 degrees.
A Fortran program to compute field values in geomagnetic dipole
coordinates from the IGRF 1965.0 field, as well as a program to transform
positions and field components from geographic to dipole coordinates and
vice versa, can be obtained from the author or from the U.S. National
Space Science Data Center, NASA Goddard Space Flight Center, Greenbelt,
Md. 20771.
REFERENCES
Chapman, Sydney, Geomagnetic nomenclature, J. Geophys. Res., 68, 1174, 1963.
IAGA Commission 2 Working Group 4 (Analysis of the Geomagnetic Field),
International Geomagnetic Reference Field 1965.0, J. Geo h . Res.,
.L4, 1+1 07-44o8, 1969.
Stern, David P., Classification of Magnetic Shells, J. Geophys.: Res.,
EO, 3629-3634, 1965.
r
I
6r
Table 1.	 IGRF 1965.0 Dipole Coefficients:
e° = 11.4350, ^.°69.761°
n	 m
m
gn
h m
n
m
9n
h m
n
1	 0 -3095; Y 16.0 v/Yr
1	 1 0 0 2.0 7.4
2	 0 -618 -18.8
2	 1 2997 2255 20.9 -1.4
2	 2 -1875 481 8.8 12.2
3	 0 906 -3.4
3	 1 -1238 -1758 -6.8 -8.o
3	 2 -1052 1170 1.3 4.6
3	 3 -546 -485 -0.5 8.0
4	 o 837 -0.3
4	 1 -496 962 1.1 -1.3
4	 2 15 176 1.2
-3.5
4	 3 311 -39 -0.3 -0.7
4	 4 -317 -312 -4.7 1.9
5	 0 -140 -0.5
5	 1 90 344 -3.9 2.1
5	 2 -322 -49 -2.0 0.3
5	 3 53 170 0.2 1.4
5	 4 -138 103 1.4 -0.7
5	 5 7 46 0.8 -o.4
6	 o 48 0.3
6	 1 -17 7 0.0 -o.4
6	 2 6o -55 -1.3 -0.3
6	 3 185 78 -0.7 -2.3
6	 4 -151 -56 -0.9 1.7
6	 5 22 -29 -0.2 -0.2
6	 6 -48 -95 -o.8 0.3
7	 0 69 o.6
7	 1 -4o -39 1.0 -1.0
7	 2 8 61 0.0 -0.1
7	 3 -14 -9 0.5 0.1
7	 4 15 36 -0.2 0.0
7	 5 3 4 -0.1 0.3
7	 6 27 -3 -0.2 -0.5
7	 7 5 10 0.2 -0.5
8	 o to -0.1
8	 1 -1 12 -0.3 0.7
8	 2 12 -2 -0.3 0.1
8	 3 -11 1 0.0 -0.1
8	 4 -14 4 -o.1 -0.2
8	 5 4 7 0.3 -0.1
8	 6 -26 10 o.6 0.2
8	 7 4 11 0.2 0.1
8	 8 -10 9 0.2 o.4
9
4
SULRUUiINE TSFUkM(GLAT,GLON,ULAT,DLON,SINU,COSD,I)
PROGRAM TO TRANSFORM GEOGRAPHIC COORDINATES TO GEOMAGNETIC DIPOLE
COORDINATES (1 = 1), OR VICE VERSA (1-2).
	
SEE MEAD(1970).
INPUT (UR, IF 1=2, OUTPUT) ARE GEOCENTRIC LATITUDE AND EAST
LONGITUDE IN DEGREES.	 OUTPUT (OR, IF 1 n 2, INPUT) ARE DIPOLE
LATITUDE AND LONGITUDE W/R/T THE DIPOLE MERIDIAN PASSING THROUGH
THE SUUTII GEOGRAPHIC POLE. SIND AND COSU ARE USED TO ROTATE THE
X AND Y FIELD COMPONENTS FROM GEOGRAPHIC TO GEOMAGNETIC OR V/V.
IF (RAU.EQ.57.29578) GO TO (1,2), 1
RAD = 57.29578
COLAT = 11.435
WLUN = b'J.7b1
STU = SIN(COLAT/RAD)
CTO = COS(CULAT/RAD)
GO TU (1, 2 ), 1
CTG = SIN(GLAT/RAD)
STG = COS(GLAT/RAD)
CTD = CTO*CTG + STO*STG*COS((GLON+IJLON)/RAD)
STU = S(1RT(1. -CTU* *2 )
ULAT = RAU * ATAN(CTU/STD)
CLU = (CTU*CTD-CTG) / (STO*STU)
SLGLO - SIN((GLON+WLON)/RAU)
SLU = STG * SLGLU / STU
ULUN = 18U. - RAU*ATAN2(SLD,-CLU)
GO TO 3
CTU = SIN(OLAT/RAD)
STU = CUS(ULAT/RAD)
CTG = CTO*CTU - STO*STU*COS(ULUN/RAD)
STG = SQRT(1.-CTG**2)
GLAT = RAU * ATAN(CTG/STG)
CLGLO - (CTU-CTU*CTG) / (STO*STG)
SLGLO - STU * SIN(D_ON/RAD) / STG
GLUM = RAD*ATAN2(SLGL0,CLGL0) - WLON
IF(GLON.LT.O.) GLUN = GLON+3b0.
CUSD = (CTO-CTD*CTG) / (STD*STG)
SINU = STU * SLGLO / STU
U = RAU * ATAN2(SIND,CUSU)
RETURN
END
SUBROUTINE
	 UIPFLD(TM,RKM,ST,CT,SL,CL,BR,BT,BP,B)
C FIELD MODEL
	 IS	 IGRF	 19b5.0	 IN	 DIPOLE	 COORDINATES	 (MEAD,1970).
C TM	 IS	 TIME	 IN	 YEARS	 (E.G.	 19119.2)
	
FOR	 WHICH	 THE	 FIELD	 IS	 DESIRED.
C INPUT
	 ARE
	 GEOCENTRIC
	
DISTANCE
	
114	 KM	 (RKM)	 AND	 SINE	 AND	 COSINE	 OF
C DIPOLE
	 COLATITUDE	 (ST,CT)	 AND	 EAST	 DIPOLE	 LONGITUDE W/R/T	 THE
C UIPULE MERIDIAN	 PASSING	 THROUGH	 THE	 SOUTH	 GEOGRAPHIC	 PULE	 (SL,CL).
C OUTPUT	 BT AND
	 BP	 (GAMHAS)	 ARE ALSU WITH RESPECT	 TO DIPOLE	 AXIS.
DIMENSION	 F(9,9),FT(9,'J),SHMIT(9,9),CONST(9,9),SP(9),CP(9),
1	 P(9,J),UP(9,J),G(9,'J)
DATA	 SlifIIT(1,1);P(1,1),CP(1),UP(1,1),SP(1)	 /	 -1.,2*1.,2*0./
DATA	 TZERU,NMAX	 /1965.0,9/
` DATA	 F/U.,-30953.,-016.,906.,837.,-140
	
,48.,69.,10.,0.,0.,2997.,
1-1238.,-496.,90.,-17.,-4U.,-1.,2255.,481.,-1875.,-1052.,15.,-322.,
2	 6U.,8.,12.,-1758.,117U.,-485.,-54b.,311.,53.,185.,-14.,-11.,962.,
3	 110.,-39.,-312.,-317.,-138.,-151.,15.,-14.,344.,-49.,170.,103.,
4	 40.,7.,22.,3.,4.,7.,-55.,78.,-56.,-29.,-95.,-48.,27.,-2b.,-39.,
5
DATA
	 FT/U.,16.,-18.8,-3.4,-.3,-.5,.3,.6,-.1,1.4,2.,20.9,-6.8,1.1,
2
3 .3, 1.4, -. 7, -. 4, . 8, -. 2, -. 1, . 3, -. 4, -. 3, -2.3, 1.1, -. 2, . 3, •. 8, -. T, . 6,
4
	
	 .2,.7,.1,-.i,-.2,-.1,.2,.1,.4,.2/
IF(SHM IT( 2,1) .EQ. -1.) GO TO 3
DO 2 N-2, NMAX
SHMIT(N,1) n SIIMIT(N-1,1) * (2*N-3) / (N-1)
DO T M-1, N
IF(M. EQ. 1) GO TO i
JJ - 1
IF(M.EQ.2) JJ - 2
SIIMIT(N,M) - SIIt4IT(N,M-1) * SQRT(FLOAT((N-M+1)*JJ)/(N+M-2))
F(M-1,N) - F(M-1,N) * SHMIT(N,M)
FT(M-1,N)	 FT(4-1,N) * SHMIT(N,M)
1 F(N,M) = F(N,M) * SHMIT(N,M)
FT(N,M) = FT(N,M) * SHMIT(N,M)
2 CONST(N,M) - FLOAT((N-2)**2-(M-1)**2) / ((2*N-3)*(2*N-5))
3 IF(TM.EQ.TMLAS+) GO TO 5
THLAST n TM
T n TM-TLERO
DO 4 N n 2,NMAX
DO 4 M-1, N
G(N,M) - F(N,M) + T*FT(N,M)
IF(M.EQ.1) GO TO 4
G(M-1,N) - F(M-1,N) + T*FT(M-1,N)
4 CONTINUE
C	 CALCULATION FOR EACH POINT USUALLY BEGINS HERE.
5 CONTINUE
SP(2) - SL
CP(2 ) - CL
DO 6 M-3, NMAX
SP(M) - SP(2)*CP(M-1) + CP(2)*SP(M-1)
6 CP(14) - CP(2)*CP(M-1) - SP(2)*SP(M-1)
AOR - 6371.16 / RKM
AR - AOR**2
BR - 0.
BT - 0.
BP - 0.
DO 11 N-2, NMAX
AR = AOR*AR
DO it M=l,N
IF(M.EQ.N) GO TO 7
P(N,M) - CT*P(N-1,M) - CONST(N,M)*P(N-2,M)
DP(N,M) = CT*DP(N-1,M) - ST*P(N-1,M) - CONST(N,M)*DP(N-2,M)
GO TO 8
7 P(N,N) = ST*P(N-1,N-1)
DP(N,N) - ST*DP(N-l,N-1) + CT*P(N-1,N-1)
8 PAR - P(N,M)*AR
IF (M.EQ.1) GO TO 9
TEMP - G(N,M)*CP(M) + G(M-1,N)*SP(M)
BP - SP - (M-1) * PAR * (G(N,M)*SP(M)-G(M-1,N)*CP(M))
GO TO 10
9 TEMP - G(N,M)
10 BT - ST + AR*DP(N,M)*TEMP
11 BR - BR - PAR*N*TEMP
BP - BP/ST
B = SQRT(BR**2+BT**2+BP**2)
RETURN
END


The international geomagnetic reference field 1965.0 in dipole coordinates

Abstract

Similar works

Full text

Available Versions

NASA Technical Reports Server