DESCRIPTION OF PROGRAM TJI95

This software package is self-contained and capable of being compiled
and executed on a variety of platforms ranging from a personal computer
to large scale \"super computers\". The software is written in FORTRAN
77. The software is designed to efficiently compute the trajectory of an
energetic charged particle of a specified momentum per unit charge
(rigidity) through a model magnetic field. For cosmic ray access to the
earth, the geocenter becomes the origin of the coordinate system. All
calculations are done in the r, , coordinate system (a right-handed,
orthogonal coordinate system). The magnetic field subroutine included in
this software is designed for efficient evaluation of the IGRF95 model
of the earth\'s magnetic field.

In its usual mode of computing the path of cosmic ray trajectories in a
model of the earth\'s magnetic field, we utilize this program to
determine the path of a cosmic ray (a positively charged particle) from
interplanetary space arriving at the earth at a specified position and
direction. To accomplish this, a negatively charged particle is
\'launched\' from the \'top\' of the atmosphere at a specified position
(latitude and longitude) in a specific direction (zenith and azimuth),
and its path traced though the model magnetic field until it either (1)
reaches a specified radial distance, (2) reenters the atmosphere, or (3)
fails to reach either condition by a specified number of iterative
steps. If the negative test particle path penetrates the specified outer
boundary (reaches interplanetary space) the direction of the particle
velocity vector at the boundary crossing is specified as asymptotic
latitude and longitude (in corresponding geocentric coordinates). If the
charged particle re-enters the atmosphere, then the re-entrant
coordinates (geocentric latitude and longitude) are given. In this
version of the software, an oxygen nuclei (^16^O) is used as the test
particle. Since rigidity is a canonical coordinate, the path of any
charged particle having the specified rigidity will be the same.

The software can be adapted to trace the path of any particle of a
specified rigidity from a specified position and direction through any
model magnetic field as long as the magnetic field is expressed as
vectors in the r, , coordinate system.

> The \'input\' to this program is a data line that specifies the
> initial position (latitude, longitude and altitude above the earth\'s
> surface), direction (azimuth and zenith), and rigidity (momentum per
> unit charge) along with control parameters to do \'N\' trajectory
> calculations at specified rigidity increments beginning at an initial
> specified rigidity.

The \'output\' of the software is in summary files with one line for
each trajectory calculation.

> The traditional summary output is called TAPE7 and is in an 80-column
> \"card image\" format. This contains the initial conditions (the
> geodetic latitude, longitude, rigidity, zenith, and azimuth), the
> final results (asymptotic latitude and longitude), fate of the
> trajectory, the number of steps in the trajectory calculation, plus a
> magnetic field identifier.
>
> There is also a second output summary, traditionally called TAPE8,
> which is in a line printer 132 column format. This contains more
> detail; the initial conditions (the geodetic and geocentric latitude,
> the longitude, rigidity, zenith, and azimuth), the final results
> (asymptotic latitude and longitude, path length (in earth radii),
> trajectory transit time, time at altitudes under 100 km, number of
> maximum and minimum in radial distance along the trajectory, the
> trajectory fate, the number of steps in the calculation, plus a
> magnetic field identifier.

This structured FORTRAN 77 software assembly consists of a main program
and four associated subroutines. The software has extensive internal
comments to aid the user in understanding the program. The program and
subroutines are:

Program TJI95 Main program; primary purpose is control

Subroutine GDGC Conversion from geodetic to geocentric coordinates

Subroutine SINGLTJ Calculates a particle trajectory

Subroutine FGRAD Evaluates the **V X B** force vectors on the particle

Subroutine MAGNEW95 Evaluates the vector magnetic field at position r, ,

[Program Organization:]{.underline}

Each subroutine has a separate unique function. Critical and often used
variables are defined in labeled common blocks. The important
\"working\" variables are in the common block WRKVLU. The trigonometric
sines and cosines are in common block WRKTSC. Definitions associated
with the shape of the ellipsoid representing the surface of the earth
are in common block GEOID. We have found that some \"super computers\"
do not allow mixing of real and integer variables in the same common
block; therefore the are two additional common blocks associated with
subroutine SINGLTJ. These are common block SNGLR (real variables) and
common block SNGLI (integer variables)

[Accuracy and Precision]{.underline}

It is recommended that the REAL\*8 precision always be used. The primary
limitation affecting the results is the accuracy of the magnetic field
expansion. For reasonably simple trajectories the results should be
repeatable, independent of the computer platform used. For long complex
trajectories, default compiler options (round off or truncate, and the
precision of intrinsic function) begin to affect the result. The
calculation procedure includes automatic error checking. The particle
acceleration terms are monitored in subroutine SINGLTJ. When significant
increases in the force on the particle are noted the step size is
reduced and the calculations are continued at smaller step intervals.
The quantity BETA ( = v/c) should be invariant throughout the
calculation and is monitored. Changes of BETA exceeding 1 part in 10^5^
results in an automatic restart and the trajectory is recalculated at
smaller step size increments.

On some \"main fames\" the intrinsic functions are automatically derived
in REAL\*8 precision; on some other systems the intrinsic functions are
evaluated in a REAL\*4 mode unless the double precision argument is
specified. In this version intended for the Desktop Computer, all
intrinsic functions are specified in the double precision mode; however,
we have left a single precision statement \"commented out\" immediately
before each double precision statement. For simple trajectories the user
probably cannot note any difference; however, for long complex
trajectories, the differences between the use of single precision
intrinsic functions and double precision intrinsic functions will become
apparent. If a specific rigidity and direction is used for comparison
and the position of each trajectory step is monitored, the effect of the
small differences between REAL\*4 and REAL\*8 accumulate and eventually
the trajectory path will differ is it is a long complex trajectory.

In the interest of computational speed, the magnetic field calculation
routine drops the evaluation of the high order terms when they make an
\"insignificant\" contribution to the total magnetic field. Again for
long complex trajectories, these small differences accumulate and the
trajectory paths may diverge when different criteria are used for
dropping magnetic field expansion.

[User Defined Parameters:]{.underline}

Two variables are intended to be user defined. These are FSTEP and
LIMIT. Default values have been set in the program, LIMIT = 600,000, and
FSTEP = 4 x 10^8^.

LIMIT is the number of Runge-Kutta steps allowed before a trajectory is
declared failed.

FSTEP is the total number of Runge-Kutta steps allowed before the run is
terminated.

Simple high rigidity trajectories often require only several hundred
steps. Simple trajectories above the upper cutoff rigidity often can be
completed in a few thousand steps. Most cosmic ray trajectories will
complete in about 10,000 steps. Some quasi-trapped periodic orbits may
require more than 100,000 steps. Trapped orbits require an infinite
number of steps. Very low rigidity trajectories initiated at high polar
latitudes will exhibit the quasi-trapped behavior and probably fail to
reach a solution. (The step size criteria is based on the time to travel
about one percent of a gyro-distance. Therefore trajectories with many
loops require many steps to complete.)

Assuming the user wants to operate in a \"batch mode\" some job control
parameters are needed. This is the quantity FSTEP. Some estimate of the
computer speed is necessary. For desktop personal computers this can
range from a few hundred steps per second on old obsolete 486 chips to
the order of 50,000 steps per second obtainable with current Pentium III
chips operating at approximately 1 GHz clock cycle time. We have found a
very significant difference in the program computational speed on the
same computer that can be attributed to the efficiency of the object
code generated by the compiler. In our testing on desktop platforms we
have found that the executable code generated by the COMPAQ Visual
Fortran operates efficiently on a Microsoft Windows operating system.
The worst performing executable code (derived from an old, no longer
sold system) ran about five times slower on the same test set of
trajectory calculation. It is assumed that workstations will have
trajectory computational speeds of the order or at least 10,000 steps
per second. The default FSTEP setting will allow a batch run of the
order of 10 hours if the program executes at 10,000 Runge-Kutta steps
per second.

[Program Operation]{.underline}

This program operates in the r, , coordinate system. The variables Y(1),
Y(2), and Y(3) are the position vectors in the r, , coordinate system
and the variables Y(4), Y(5), and Y(6) are the velocity vectors.

The program initially defines the physical constants used in the
calculation and control parameters. It then enters a control loop
beginning with reading a data line to determine the initial position and
direction, the specified starting rigidity and how many trajectories to
calculate at specified increments.

For each control line read, a call to subroutine GDGC converts the
initial geodetic coordinates (map makers coordinates on the earth\'s
surface) to geocentric r, , coordinates. Then the trajectory
calculations are done by subroutine SINGLTJ

The control loop continues (read in control line, convert coordinates,
do trajectory calculations) until a negative (or zero) value of rigidity
is read in. When this occurs the program terminates.

[Labeled Common arguments:]{.underline}

Block name: /WRKVLU/

Arguments in block F(6), Y(6), ERAD, EOMC, VEL, BR, BT, BP

F(6) Array of force and acceleration vectors

Y(6) Array of position and velocity vectors

ERAD Average radius of the earth in kilometers

EOMC Electronic charge divided by relativistic mass energy equivalent
(mixed units)

VEL Particle velocity in earth radii per second

BP Value of the B() magnetic field vector (in units of Gauss)

BR Value of the B(r) magnetic field vector (in units of Gauss)

BT Value of the B() magnetic field vector (in units of Gauss)

Block name: /WRKTSC/

Arguments in block TSY2, TCY2, TSY3, TCY3

TSY2 Sine of the Y(2) coordinate (theta coordinate)

TCY2 Cosine of the Y(2) coordinate (theta coordinate)

TSY3 Sine of the Y(3) coordinate (phi coordinate)

TCY3 Cosine of the Y(3) coordinate (phi coordinate)

Block name: /TRIG/

Arguments in block PI, RAD, PI02

PI Value of pi

RAD Value of degrees in a radian

PI02 Value of pi/2.0

Block name: /GEOID/

Arguments in block ERADPL, ERECSQ

ERADPL Polar radius of the earth in kilometers

ERECSQ Eccentricity of ellipsoid squared

Block name: /SNGLR/

Arguments in block SALT, DISOUT, GCLATD, GDLATD, GLOND, GDAZD, GDZED,

RY1, RY2, RY3, RHT, TSTEP

SALT Start altitude of trajectory above surface of geoid

DISOUT Radial distance (in earth radii) for termination of calculation

GCLATD Geocentric latitude in degree

GDLATD Geodetic latitude in degrees

GLOND East longitude in degrees

GDAZD Geodetic azimuth in degrees

GDZED Geodetic zenith in degrees

RY1 Original start position Y(1), (radial component in the r, ,
coordinate system)

RY2 Original start position Y(2), (theta component in the r, ,
coordinate system)

RY3 Original start position Y(3), (phi component in the r, , coordinate
system)

RHT Height above geoid where trajectory re-enters the atmosphere

TSTEP Total number of steps in this run.

Block name: /SNGLI/

Arguments in block LIMIT, NTRAJC, IERRPT

LIMIT Maximum number of steps before \'failed\' trajectory

NTRAJC Number of trajectories calculated in this run

IERRPT Integer control for printing diagnostics (normally = 0)

[Subroutines called:]{.underline}

GCGC (TCD, TSD)

TCD Cosine of the rotation angle

TSD Sine of the rotation angle

SINGLTJ (PC, IRSLT, INDXPC, Y1GC, Y2GC, Y3GC)

PC Particle rigidity (in units of GV)

IRSLT Integer result of trajectory calculation

(+1 = Allowed, 0 = Failed, -1 = Re-entrant)

INDXPC Integer value of PC in units of MV

Y1GC Y(1) position in geocentric coordinates

Y2GC Y(2) position in geocentric coordinates

Y3GC Y(3) position in geocentric coordinates

[Dimensioned variables:]{.underline} all in labeled common

F(6) Array of force and acceleration vectors

Y(6) Array of position and velocity vectors

[Data files]{.underline}: none

[Output files:]{.underline}

TAPE7 from subroutine SINGTJ (80 character summary)

TAPE8 from subroutine SINGTJ (132 character line printer summary)

TAPE16 (diagnostic output; if desired, set IERRPT to \> 0)

[Listing of all variables]{.underline} in program TJI95

B Magnitude of the magnetic field (in Gauss) \[in COMMON /WRKVLU/\]

BP Value of the B() magnetic field vector (in Gauss) \[in COMMON
/WRKVLU/\]

BR Value of the B(r) magnetic field vector (in Gauss) \[in COMMON
/WRKVLU/\]

BT Value of the B() magnetic field vector (in Gauss) \[in COMMON
/WRKVLU/\]

CF Character variable \"F\"

DELPC Increments of rigidity spacing search in control line

DISOUT Distance (in earth radii) for trajectory termination \[in COMMON
/SNGLR/\]

ERAD Average radius of the earth in kilometers \[in COMMOM /WRKVLU/\]

ERADPL Polar radius of the earth in km \[in COMMON /GEOID/\]

ERECSQ Eccentricity of ellipsoid squared \[in COMMON /GEOID/\]

F(6) Array of \"F\" values (velocity and acceleration in program
coordinates)

\[in COMMON /WRKVLU/\]

GCLATD Geocentric latitude in degrees (+ = north, south = -) \[in COMMON
/SNGLR/\]

GDAZ Geodetic azimuth in radians (measured clockwise from north)

GDAZD Geodetic azimuth in degrees (measured clockwise from north)

GDLATD Geodetic latitude in degrees (+ = north, south = -) \[in COMMON
/SNGLR/\]

GDZE Geodetic zenith in radians (0 = vertical)

GDZED Geodetic zenith in degrees (0 = vertical)

GLOND Geodetic East longitude in degrees (from Greenwich meridian)

\[in COMMON /SNGLR/\]

IDELPC Integer value of rigidity change increment in MV (attempt to
avoid round off)

IERRPT Integer control for printing diagnostics (normally set to 0)

INDEX Arbitrary index number of input control line (optional)

INDO Integer control of number of trajectories to calculate

INDXPC Integer value of rigidity in MV increments (attempt to avoid
round off)

IOSTAT Integer system argument of status of read

IRSLT Internal result of particle trajectory (+1 = allowed, 0 = failed,
-1 = re-entrant)

ISALT Integer value of start altitude (in km) above geoid surface

LIMIT Limit of number of steps before trajectory is declared \"Failed\"

LSTEP Number of times step size control has been reduced to overcome
trajectory error

NDO Integer control read in (number of trajectories to compute from this
control line)

NTRAJC Number of trajectories in this computer run

PC Rigidity of particle (in units of GV)

PI Real value of the quantity Pi (\~3.1415926535) \[in COMMON /TRIG/\]

PIO2 Real value of Pi divided by 2.0 \[in COMMON /TRIG/\]

RAD Real value of one radian (\~57.29578 degrees) \[in COMMON /TRIG/\]

RHT Height above geoid where a trajectory re-entered the atmosphere

RY1 Real value of the starting position of the r coordinate in r, ,
coordinates

\[in COMMON /SNGLR/\]\]

RY2 Real value of the starting position of the theta coordinate in r, ,
coordinates

\[in COMMON /SNGLR/\]

RY3 Real value of the starting position of the phi coordinate in r, ,
coordinates

\[in COMMON /SNGLR/\]

SALT Real value of the starting altitude above the surface of the geoid

\[in COMMON /SNGLR/\]

TCD Trigonometric cosine of the rotation angle from geodetic to
geocentric

TCGDAZ Trigonometric cosine of the geodetic azimuth

(Measured clockwise from north) \[in COMMON /SNGLR/\]

TCGDZE Trigonometric cosine of the geodetic zenith (measured clockwise
from north)

TCY2 Trigonometric cosine of the vector theta angle in r, , coordinates

\[in COMMON /SNGLR/\]

TCY3 Trigonometric cosine of the vector phi angle in r, , coordinates

\[in COMMON /SNGLR/\]

TSTEP Number steps executed in this run

TSY2 Trigonometric sine of the vector theta () angle in r, , coordinates

\[in COMMON /SNGLR/\]

TSY3 Value of the trigonometric sine of the vector phi () angle in r, ,
coordinates

\[in COMMON /SNGLR/\]

VEL Particle velocity in earth radii per second \[in COMMON /WRKVLU/\]

Y(6) Array of \"Y\" values (position and velocity in r, , coordinates)

\[in COMMON /WRKVLU/\]

Y1GC Starting position r component in geocentric coordinates

Y1GD Starting position r component in geodetic coordinates

Y2GC Starting position component in geocentric coordinates

Y2GD Starting position component in geodetic coordinates

Y3GC Starting position component in geocentric coordinates

Y3GD Starting position component in geodetic coordinates

Subroutine GDGC (TCD, TSD)

This subroutine calculates the angle between geodetic and geocentric
coordinates. The arguments TCD and TSD are the trigonometric cosine and
sine of the rotation angle from a normal from the surface of the geoid
(geodetic coordinates) and a radial from the center of the earth
(geocentric coordinates). See Appendix B of NSSDC 72-12)

[Arguments in call statement]{.underline}

TCD Cosine of the rotation angle

TSD Sine of the rotation angle

[Labeled Common arguments:]{.underline}

Block name: /WRKVLU/

Arguments in block F(6), Y(6), ERAD, EOMC, VEL, BR, BT, BP

F(6) Array of force and acceleration vectors

Y(6) Array of position and velocity vectors

ERAD Average radius of the earth in kilometers

EOMC Electronic charge divided by relativistic mass energy equivalent
(mixed units)

VEL Particle velocity in earth radii per second

BP Value of the B() magnetic field vector (in units of Gauss)

BR Value of the B(r) magnetic field vector (in units of Gauss)

BT Value of the B() magnetic field vector (in units of Gauss)

Block name: /WRKTSC/

Arguments in block TSY2, TCY2, TSY3, TCY3

TSY2 Sine of the Y(2) coordinate (theta coordinate)

TCY2 Cosine of the Y(2) coordinate (theta coordinate)

TSY3 Sine of the Y(3) coordinate (phi coordinate)

TCY3 Cosine of the Y(3) coordinate (phi coordinate)

Block name: /TRIG/

Arguments in block PI, RAD, PI02

PI Value of pi

RAD Value of degrees in a radian

PI02 Value of pi/2.0

Block name: /GEOID/

Arguments in block ERADPL, ERECSQ

ERADPL Polar radius of the earth in kilometers

ERECSQ Eccentricity of ellipsoid squared

Block name: /SNGLR/

Arguments in block SALT, DISOUT, GCLATD, GDLATD, GLOND, GDAZD, GDZED,

RY1, RY2, RY3, RHT, TSTEP

SALT Start altitude of trajectory above surface of geoid

DISOUT Radial distance (in earth radii) for termination of calculation

GCLATD Geocentric latitude in degrees

GDLATD Geodetic latitude in degrees

GLOND East longitude in degrees

GDAZD Geodetic azimuth in degrees

GDZED Geodetic zenith in degrees

RY1 Original start position Y(1)

RY2 Original start position Y(2)

RY3 Original start position Y(3)

RHT Height above geoid where trajectory re-enters the atmosphere

TSTEP Total number of steps in this run.

Block name: /SNGLI/

Arguments in block LIMIT, NTRAJC, IERRPT

LIMIT Maximum number of steps before \'failed\' trajectory

NTRAJC Number of trajectories calculated in this run

IERRPT Integer control for printing diagnostics (normally = 0)

[Dimensioned variables]{.underline}: all in labeled common

F(6) Array of force and acceleration vectors

Y(6) Array of position and velocity vectors

[Subroutines called]{.underline}: none

[Data files]{.underline}: none

[Output files]{.underline}: none

[Operation]{.underline}:

The shape of the earth used is not a sphere, but an ellipsoid having a
specified polar radius, equatorial radius, and eccentricity. When this
subroutine is called, it defines the shape of an oblate earth from the
polar and equatorial radius, and calculates vectors from a normal on the
surface of the ellipsoid to the specified position in geodetic
coordinates, at the specified latitude, and determines the vector
rotation angle between geodetic coordinates and geocentric coordinates.
The sine and cosine of this rotation angle are passed to the calling
program. Geodetic latitude is a measure of latitude in a coordinate
system normal to the surface of the earth. At a position on or above the
surface of the ellipsoid, there is a slight difference between a
direction normal to the surface of the ellipsoid and a direction to the
geocentric. This difference is latitude dependent. (It is zero at the
equator or poles and can be as large as approximately 1/2 of a degree at
mid latitudes.) The vector rotation angle allows for direction
specification in both geodetic (map) coordinates and geocentric
coordinates. This small correction for the direction may be
insignificant for some applications, but may be significant for
precision calculation in a specific direction at high rigidities.

[Data checking]{.underline}: none. The data to describe the shape of the
earth are included in the subroutine.

[Listing of all variables]{.underline} used in subroutine GDGC of
program TJI95

B Magnitude of the magnetic field (in Gauss) \[in COMMON /WRKVLU/\]

BP Value of the B() magnetic field vector (in Gauss) \[in COMMON
/WRKVLU/\]

BR Value of the B(r) magnetic field vector (in Gauss) \[in COMMON
/WRKVLU/\]

BT Value of the B() magnetic field vector (in Gauss) \[in COMMON
/WRKVLU/\]

DISOUT Distance (in earth radii) for trajectory termination \[in COMMON
/SNGLR/\]

DISTKM Starting position geocentric radial distance from geocenter

EOMC Electronic charge divided by relativistic mass energy equivalent
(mixed units)

\[in COMMON /WRKVLU/\]

ERAD Average radius of the earth in kilometers \[in COMMOM /WRKVLU/\]

ERADPL Polar radius of the earth in km \[in COMMON /GEOID/\]

ERECSQ Eccentricity of ellipsoid squared \[in COMMON /GEOID/\]

ERPLSQ Polar radius of earth (in km) squared

F(6) Array of \"F\" values (velocity and acceleration in program
coordinates)

\[in COMMON /WRKVLU/\]

GCLATD Geocentric latitude in degrees (+ = north, south = -) \[in COMMON
/SNGLR/\]

GDAZD Geodetic azimuth in degrees (measured clockwise from north)

GDCLT Geodetic co-latitude (in radians)

GDLATD Geodetic latitude in degrees (+ = north, south = -) \[in COMMON
/SNGLR/\]

GDZED Geodetic zenith in degrees (0 = vertical)

GLOND Geodetic East longitude in degrees (from Greenwich meridian)

\[in COMMON /SNGLR/\]

ONE Intermediate term in computations (see NSSDC ALLMAG description)

PI Real value of the quantity Pi (\~3.1415926535) \[in COMMON /TRIG/\]

PIO2 Real value of Pi divided by 2.0 \[in COMMON /TRIG/\]

RAD Real value of one radian (\~57.29578 degrees) \[in COMMON /TRIG/\]

RHO Intermediate term in computations (see NSSDC ALLMAG description)

RHT Height above geoid where a trajectory is declared re-entrant

\[in COMMON /SNGLR/\]

RY1 Real value of the starting position of the r coordinate in r, ,
coordinates

\[in COMMON /SNGLR/\]

RY2 Real value of the starting position of the theta coordinate in r, ,
coordinates

\[in COMMON /SNGLR/\]

RY3 Real value of the starting position of the phi coordinate in r, ,
coordinates

\[in COMMON /SNGLR/\]

SALT Real value of the starting altitude above the surface of the geoid

\[in COMMON /SNGLR/\]

TCD Trigonometric cosine of the rotation angle from geodetic to
geocentric

TCGDCLT Trigonometric cosine of the geocentric co-latitude

TCY2 Trigonometric cosine of the vector theta angle in r, , coordinates

\[in COMMON /SNGLR/\]

TCY3 Trigonometric cosine of the vector phi angle in r, , coordinates

\[in COMMON /SNGLR/\]

THREE Intermediate term in computations (see NSSDC ALLMAG description)

TSD Trigonometric sine of the rotation angle from geodetic to geocentric

TSGDCLT Trigonometric sine of the geocentric co-latitude

TSTEP Number steps executed in this run

TSY2 Trigonometric sine of the vector theta () angle in r, , coordinates

\[in COMMON /SNGLR/\]

TSY3 Trigonometric sine of the vector phi () angle in r, , coordinates

\[in COMMON /SNGLR/\]

TWO Intermediate term in computations (see NSSDC ALLMAG description)

VEL Particle velocity in earth radii per second \[in COMMON /WRKVLU/\]

Y(6) Array of \"Y\" values (position and velocity in r, , coordinates)

> \[in COMMON /WRKVLU/\]

Reference publication:

[ALLMAG, GCAZLMG, LINTRA: Commuter Programs For Geomagnetic Field And
Field-Line Calculations]{.underline}, E.G. Stassinopulous and G.D. Mead,
NSSDC 72-12, February 1972, NASA, GFSC, Greenbelt MD.

Subroutine SINGLTJ (PC, IRSLT, INDXPC, Y1GC, Y2GC, Y3GC)

This subroutine does the actual trajectory tracing. When called it
initially defines control parameters and constants used in the particle
tracing and initializes the Runge-Kutta variables to zero. It sets up
the initial position and direction, and defines the relativistic
parameters relating to the particle total energy and speed.

In this version of the subroutine, an oxygen nuclei (^16^O) is used as
the test particle. By definition a ^16^O nuclei has a mass of 16 Atomic
Mass Units (AMU) and an atomic charge of 8. The mass-energy conversion
for one AMU is 0.93114 GeV. If it were desired to modify the program for
some other nuclei, such as a proton that has an atomic charge of 1 and
atomic mass of 1.0081415 AMU, then the rest mass energy for atomic
nuclei must be adjusted.

After the initial definitions, the subroutine then chooses an initial
starting step length (a relatively small value) and starts the
Runge-Kutta process of tracing the particle trajectory. After each step
it goes though an error checking and detection process. If the checks
are satisfactory, it determines the particle location with respect to
the atmosphere and the outer boundary.

> If the charged particle is between the atmosphere and the outer
> boundary, it adjusts the size of the next step and continues the
> trajectory tracing until the LIMIT on the number of steps is reached.
>
> If the charged particle is entering the atmosphere, it terminates the
> calculation.
>
> If the charged particle is less than 100 km above the earth\'s
> surface, it maintains a running sum of the time at low altitudes.
>
> If the charged particle is approaching the outer boundary, it adjusts
> the step size so it penetrates this boundary at small step lengths.
>
> If the charged particle has penetrated the outer boundary at a small
> step, it computes the final coordinates.

When the particle has reached a solution (allowed or re-entrant) or
reached the step limit, it writes out the result and returns to the
calling program.

[Arguments in call statement]{.underline}

PC, IRSLT, INDXPC, Y1GC, Y2GC, Y3GC

PC Particle rigidity

IRSLT Integer result of trajectory (+1 = allowed, 0 = failed, -1 =
re-entrant)

INDXPC Integer value of PC in MV units

Y1GC Y(1) position in geocentric coordinates

Y2GC Y(2) position in geocentric coordinates

Y3GC Y(3) position in geocentric coordinates

[Labeled Common arguments]{.underline}:

Block name: /WRKVLU/

Arguments in block F(6), Y(6), ERAD, EOMC, VEL, BR, BT, BP

F(6) Array of force and acceleration vectors

Y(6) Array of position and velocity vectors

ERAD Average radius of the earth in kilometers

EOMC Electronic charge divided by relativistic mass energy equivalent
(mixed units)

VEL Particle velocity in earth radii per second

BP Value of the B() magnetic field vector (in units of Gauss)

BR Value of the B(r) magnetic field vector (in units of Gauss)

BT Value of the B() magnetic field vector (in units of Gauss

Block name: /WRKTSC/

Arguments in block TSY2, TCY2, TSY3, TCY3

TSY2 Sine of the Y(2) coordinate (theta coordinate)

TCY2 Cosine of the Y(2) coordinate (theta coordinate)

TSY3 Sine of the Y(3) coordinate (phi coordinate)

TCY3 Cosine of the Y(3) coordinate (phi coordinate)

Block name: /TRIG/

Arguments in block PI, RAD, PI02

PI Value of pi

RAD Value of degrees in a radian

PI02 Value of pi/2.0

Block name: /GEOID/

Arguments in block ERADPL, ERECSQ

ERADPL polar radius of the earth in kilometers

ERECSQ Eccentricity of ellipsoid squared

Block name: /SNGLR/

Arguments in block SALT, DISOUT, GCLATD, GDLATD, GLOND, GDAZD, GDZED,

RY1, RY2, RY3, RHT, TSTEP

SALT Start altitude of trajectory above surface of geoid

DISOUT Radial distance (in earth radii) for termination of calculation

GCLATD Geocentric latitude in degrees

GDLATD Geodetic latitude in degrees

GLOND East longitude in degrees

GDAZD Geodetic azimuth in degrees

GDZED Geodetic zenith in degrees

RY1 Original start position Y(1), (radial component in the r, ,
coordinate system)

RY2 Original start position Y(2) , (theta component in the r, ,
coordinate system)

RY3 Original start position Y(3) , (phi component in the r, , coordinate
system)

RHT Height above geoid where trajectory re-enters the atmosphere

TSTEP Total number of steps in this run.

Block name: /SNGLI/

Arguments in block LIMIT, NTRAJC, IERRPT

LIMIT Maximum number of steps before \'failed\' trajectory

NTRAJC Number of trajectories calculated in this run

IERRPT Integer control for printing diagnostics (normally = 0)

[Dimension variables]{.underline}: (not in labeled common)

P(6), Q(6), R(6), S(6), YB(6), FOLD(6), YOLD(6)

P(6) Runge-Kutta variable

Q(6) Runge-Kutta variable

R(6) Runge-Kutta variable

S(6) Runge-Kutta variable

YB(6) Runge-Kutta variable

FOLD(6) \"F\" vectors of previous step

YOLD(6) \"Y\" vectors of previous step

[Subroutines called]{.underline}: FGRAD

[Data files]{.underline}: none

[Output files]{.underline}:

TAPE7 (80 character summary)

TAPE8 (132 character line printer summary)

TAPE16 (diagnostic output; if desired set IERRPT to \> 0)

[Program Operation]{.underline}:

This program operates in the r, , coordinate system. The variables Y(1),
Y(2), and Y(3) are the position vectors in the r, , coordinate system
and the variables Y(4), Y(5), and Y(6) are the velocity vectors.

When this subroutine is called, it initially defines control parameters
and constants used in the particle path tracing, and initializes the
Runge-Kutta variables to zero. It obtains the particle\'s height with
respect to the surface of an oblate earth. It sets up the initial
position vectors, Y(1), Y(2) and Y(3), and based on the particle
rigidity, sets up velocity vectors, Y(4), Y(5), and Y(6). It then
defines the relativistic parameters TENG (total energy), EOMC (charge
per relativistic mass/energy equivalent), and GMA (the relativistic
parameter of total energy over the rest mass energy). It defines scalar
quantities relating to the particle, BETA (the particle speed with
respect to light), PVEL (the particle speed in earth radii per second),
and HMAX (a maximum step length allowed for this particle rigidity).

Next it defines an initial starting step length (a relatively small
value) and starts the Runge-Kutta process of tracing the particle
trajectory. Comment cards specifically indicate the Runge-Kutta
iteration process, which is the coding between FORTRAN statement numbers
130 and 170. The calls to subroutine FGRAD evaluate the **V X B** force
on the particle during this step. The logic is very similar to that
documented in Ralston and Wilf (1960). After each Runge-Kutta iteration
step there is an extensive error checking and detection process.

The error checking process begins with a check on the particle speed
(BETA), which should remain invariant throughout the trajectory. If the
difference between the initial particle speed (BETA) and it\'s current
speed (RCKBETA) is greater than EDIF, then the trajectory tracing
process is re-initialized (including the NSTEP variable) and the
trajectory re-started at a smaller step size selection criteria. Up to
five re-starts are allowed before the specific trajectory is declared
impossible to calculate, evaluated as \"failed\", and the path length
made negative in order to distinguish it from successful trajectories.
In order to attempt to reach a solution the EDIF variable is widened by
a factor of two after each successive trajectory failure.

After the error check, then the acceleration of the particle is compared
with previous values. We have found that computational errors are most
likely to occur when there are rapid changes in the acceleration. If the
average change in acceleration exceeds a factor of five, or if any
component of the acceleration exceeds a factor of three, then the step
length for the next Runge-Kutta step is reduced.

Along the particle path the software checks the particle location with
respect to the atmosphere and the outer boundary. If the charged
particle is less than 100 km above the earth\'s surface, it maintains a
running sum of the time at low altitudes. If the charged particle is
entering the atmosphere, it terminates the calculation.

The next check determines if the particle has penetrated the outer
termination boundary. The step length can be relatively large at extreme
distances from the earth. If the outer boundary has been penetrated at a
large step size, the trajectory is \"backed up\" and the step size
reduced until it penetrates the boundary at a small step size. This
results in a more precise determination of the penetration location and
can significantly affect the computed asymptotic direction.

If there are no errors and the charged particle is between the
atmosphere and the outer boundary, the software adjusts the size of the
next step appropriate for the magnitude of the magnetic field (the step
size is normally about one percent of the gyro-distance) and continues
the trajectory tracing. The basic step length algorithm is:

> H = ((2.0\*\*33.333\*PC)/(B\*\*C))/100.0\
> \
> where \"H\" is time in seconds, \"PC\" is the particle rigidity in GV,
> \"B\" is the magnitude of the magnetic field in Gauss and \"C\" is the
> speed of light in km/sec.
>
> (A handy formula to remember is the gyro-radius is 33 km per GV per
> Gauss)

The software initially starts at a trajectory calculation at a small
step size and the step size is permitted to grow at a maximum of about
20 percent each step. If the particle trajectory starts to loop to a
lower altitude, then the step size is reduced to compensate for the
increasing magnitude of the magnetic field. When a loop in the
trajectory develops, the acceleration forces increase and step size
adjustments are made in the case of a significant increase in the
acceleration forces.

During the trajectory tracing, the software notes the number of maximum
and minima the trajectory experiences. This information is useful in
ascertaining the complexity of the trajectory.

When the particle has reached a solution (allowed or re-entrant) or
reached the step limit, the subroutine will write out the result and
return to the calling program. The fate is coded in the variable IFATE:
(0 = Allowed, 1 = Failed, 2 = Re-entrant, 3 = Failed, but max alt \> 6.6
earth radii)

If the trajectory is allowed (penetrates the outer boundary), then the
velocity vectors are transformed into asymptotic latitude and longitude.
Asymptotic latitude and longitude are the geocentric coordinates the
velocity vector would have at infinity. If the trajectory re-enters the
atmosphere, then the position coordinates are transformed to geocentric
latitude and longitude.

The \'output\' of the software is in summary files, on one line for each
trajectory calculation.

The output files are:

> TAPE7 is in an 80-column \"card image\" format. This contains the
> initial conditions (the geodetic latitude, longitude, rigidity,
> zenith, and azimuth), the final results (asymptotic latitude and
> longitude, fate and number of steps), and a magnetic field identifier.
>
> TAPE8 is in a line printer 132-column format. This contains more
> detail; the initial conditions (the geodetic latitude and the
> geocentric latitude, the longitude, rigidity, zenith, and azimuth),
> the final results (asymptotic latitude and longitude, path length,
> trajectory time, time at altitudes under 100 km, number of maximum and
> minimum in trajectory, fate, and number of steps), and a magnetic
> field identifier.
>
> TAPE16 is a diagnostic output, including a record of restarts due to
> BETA checks\
> or trajectory failures.

See the examples section for samples of the output.

[Possible Additions for Trajectory Plotting]{.underline}

If it is desired to plot a trajectory, the position variables Y(1),
Y(2), and Y(3) must be stored after each Runge-Kutta step in a suitable
array. The task of adding such a modification should be straightforward
and is left to the individual program user.

[Listing of all variables]{.underline} in subroutine SINGLTJ

ACCER Magnitude of current value of the particle acceleration

ACCOLD Magnitude of last value of the particle acceleration

AFOLD Absolute value of FOLD

AHLT Variable to control step size at high latitude

ANUC Atomic number of number of nucleons in atom

ATRG1 Intermediate value to compute asymptotic latitude (avoid 0/0
problem)

ATRG2 Intermediate value to compute asymptotic latitude (avoid 0/0
problem)

AZD Azimuth angle in degrees (measured clockwise from north)

B Magnitude of the magnetic field (in Gauss) \[in COMMON /WRKVLU/\]

BETA Particle speed as fraction of light speed ( = v/c where c is speed
of light)

BETAST Control variable for reducing step size if error has occurred

BP Value of the B() magnetic field vector (in Gauss) \[in COMMON
/WRKVLU/\]

BR Value of the B(r) magnetic field vector (in Gauss) \[in COMMON
/WRKVLU/\]

BT Value of the B() magnetic field vector (in Gauss) \[in COMMON
/WRKVLU/\]

CF Character variable \"F\"

CNAME Character variable (up to 6 characters) identifying magnetic field
used

CR Character variable \"R\"

DELACC Change in particle acceleration from previous step

DISCK Step length control variable to approach boundary at 10 %
increments

DISOUT Distance (in earth radii) for trajectory termination \[in COMMON
/SNGLR/\]

DISTR Radial distance (in earth radii) from current distance to
termination boundary

EDIF Variation in allowed before error declared

EMCSQ Mass energy equivalent for a AMU (0.931141 GeV)

EOMC Electronic charge divided by relativistic mass energy equivalent
(mixed units)

\[in COMMON /WRKVLU/\]

ERAD Average radius of the earth in kilometers \[in COMMOM /WRKVLU/\]

ERADPL Polar radius of the earth in km \[in COMMON /GEOID/\]

ERECSQ Eccentricity of ellipsoid squared \[in COMMON /GEOID/\]

F(6) Array of \"F\" values (velocity and acceleration in program
coordinates)

\[in COMMON /WRKVLU/\]

FOLD(6) Array of \"F\" values (velocity and acceleration) from previous
step

FASLAT Asymptotic latitude (in degrees)

FASLON Asymptotic longitude (in degrees east of the Greenwich meridian)

GCLATD Geocentric latitude in degrees (+ = north, south = -) \[in COMMON
/SNGLR/\]

GDAZD Geodetic azimuth in degrees (measured clockwise from north)

GDLATD Geodetic latitude in degrees (+ = north, south = -) \[in COMMON
/SNGLR/\]

GDZED Geodetic zenith in degrees (0 = vertical)

GLOND Geodetic East longitude in degrees (from Greenwich meridian)

\[in COMMON /SNGLR/\]

GMA Relativistic factor (total energy/rest energy)

GRNDKM Altitude above surface of earth at this latitude (in km)

H Runge-Kutta step size (in seconds)

HB Preliminary value of step size for = 1

HCK Control to limit step size growth to 20%

HCNG Change of step size from previous step

HMAX Maximum value of step size allowed

HOLD Value of previous Runge-Kutta step size

HSNEK Control to approach 90% of distance to boundary

HSTART Starting step size value (deliberately made small)

I Index variable in do loops

IAZ Integer value of azimuth (measured counter clockwise from north)

ICK Index for checking acceleration growth

IERRPT Integer control for printing diagnostics (normally set to 0)

IFATE Integer fate of particle trajectory ( 0 = Allowed, 1 = Failed, 2 =
Re-entrant

3 = Failed, but max alt \> 6.6 earth radii)

INDXPC Index of particle rigidity in MV

IRT Integer control for writing results (+1 = Allowed, 0 = Failed, -1 =
Re-entrant)

IRSLT Internal result of particle trajectory (+1 = Allowed, 0 = Failed,
-1 = Re-entrant)

ISALT Integer value of start altitude (in km) above earth surface

IZE Integer value of zenith angle (in degrees)

KBF Number of failed attempts to trace this trajectory

LIMIT Limit of number of steps before trajectory declared \"Failed\"

LSTEP Number of times the step size control reduced to overcome
trajectory error

NMAX Number of maxima in complex trajectory path

NMIN Number of minima in complex trajectory path

NSTEP Number of steps in current trajectory

NSTEPT Temporary variable that can be used to print out first 1000 steps

NTRAJC Number of trajectories in this computer run

P(6) Array of intermediate values used in Runge-Kutta integration

PATH Total distance of trajectory path from start to termination

PC Rigidity of particle (in units of GV)

PI Real value of the quantity Pi (\~3.1415926535) \[in COMMON /TRIG/\]

PIO2 Real value of Pi divided by 2.0 \[in COMMON /TRIG/\]

PSALT Current particle distance from ground (used for re-entrant
calculations)

PTCY2 Absolute value of cosine Y(2) (used in control of polar step size)

PVEL Particle velocity (in earth radii per second)

R(6) Array of intermediate values used in Runge-Kutta integration

R100KM Y(1) distance of 100 km altitude at this latitude

R120KM Y(1) distance of 120 km altitude at this latitude

RAD Real value of one radian (\~57.29578 degrees) \[in COMMON /TRIG/\]

RC1O6 Constant in Runge-Kutta integration (1.0/6.0)

RCKBETA Current value of particle after this step

RENLAT Latitude of re-entrant particle intersection with atmosphere

RENLON Longitude of re-entrant particle intersection with atmosphere

RFA Ratio of acceleration magnitude between current step and last step

RFCK Ratio of acceleration component between current step and last step

RHT Height above geoid where a trajectory re-entered the atmosphere

RY1 Real value of the starting position of the r coordinate in r, ,
coordinates

\[in COMMON /SNGLR/\]

RY2 Real value of the starting position of the theta coordinate in r, ,
coordinates

\[in COMMON /SNGLR/\]

RY3 Real value of the starting position of the phi coordinate in r, ,
coordinates

\[in COMMON /SNGLR/\]

S(6) Array of intermediate values used in Runge-Kutta integration

SALT Real value of the starting altitude above the surface of the geoid

\[in COMMON /SNGLR/\]

SR2 Runge-Kutta constant (square root of 2.0)

TAU Time (in seconds) for a trajectory transit from start to termination

TBETA Difference between current value of and starting value of

TCY2 Trigonometric cosine of the vector theta () angle in r, ,
coordinates

\[in COMMON /SNGLR/\]

TCY3 Trigonometric cosine of the vector phi () angle in r, , coordinates

\[in COMMON /SNGLR/\]

TENG Total energy of particle (kinetic energy plus rest mass energy)

TMS2O2 Runge-Kutta constant (2.0 - SR2/2.0)

TPS2O2 Runge-Kutta constant ((2.0 + SR2/2.0)

TSTEP Number steps executed in this run

TSY2 Trigonometric sine of the vector theta () angle in r, , coordinates

\[in COMMON /SNGLR/\]

TSY2SQ Square of TSY2

TSY3 Value of the trigonometric sine of the vector phi () angle in r, ,
coordinates

\[in COMMON /SNGLR/\]

TU100 Time (in seconds) the particle is under 100 km altitude

VEL Particle velocity in earth radii per second \[in COMMON /WRKVLU/\]

Y(6) Array of \"Y\" values (position and velocity in r, , coordinates)

\[in COMMON /WRKVLU/\]

YB(6) Array of intermediate values used in Runge-Kutta integration

Y10 Y(1) radial coordinate for re-entrant distance at this latitude

Y1GC Starting position r component in geocentric coordinates

Y2GC Starting position component in geocentric coordinates

Y3GC Starting position component in geocentric coordinates

YDA5 Intermediate value for computing asymptotic latitude

YMAX Maximum radial distance attained by trajectory

ZCHARGE Atomic charge number

ZED Zenith angle in degrees

Reference publication:

Ralston, A, and Wilf, S.H., [Mathematical Models for Digital
Computers]{.underline}, JohnWiley and Sons, New York, 1960.

Subroutine FGRAD

This subroutine calculates the **V X B** force on the charged particle.

[Arguments in call statement]{.underline}: none. All in labeled common

[Labeled Common arguments]{.underline}:

Block name: /WRKVLU/

Arguments in block F(6), Y(6), ERAD, EOMC, VEL, BR, BT, BP

F(6) Array of force and acceleration vectors

Y(6) Array of position and velocity vectors

ERAD Average radius of the earth in kilometers

EOMC Electronic charge divided by relativistic mass energy equivalent
(mixed units)

VEL Particle velocity in earth radii per second

BP Value of the B() magnetic field vector (in units of Gauss)

BR Value of the B(r) magnetic field vector (in units of Gauss)

BT Value of the B() magnetic field vector (in units of Gauss)

Block name: /WRKTSC/

Arguments in block TSY2, TCY2, TSY3, TCY3

TSY2 Sine of the Y(2) coordinate (theta coordinate)

TCY2 Cosine of the Y(2) coordinate (theta coordinate)

TSY3 Sine of the Y(3) coordinate (phi coordinate)

TCY3 Cosine of the Y(3) coordinate (phi coordinate)

[Subroutines called]{.underline}:

MAGNEW95

(All arguments are in labeled common /WRKVLU/ and /WRKTSC/)

[Data files]{.underline}: none

[Output files]{.underline}: none

[Program Operation]{.underline}:

When this subroutine is called, the force vectors, ((F(1), F(2), F(3))
are defined. The sine and cosine of the Y(2) coordinates are determined.
The magnetic field vectors (BR, BT, BP) at the particle position (Y(1),
Y(2), Y(3)) are obtained by the call to subroutine MAGNEW95. Then the
acceleration vectors (F(4), F(5), F(6)) are calculated.

[Listing of all variables]{.underline} in subroutine FGRAD in program
TJI95

B Magnitude of the magnetic field (in Gauss) \[in COMMON /WRKVLU/\]

BP Value of the B() magnetic field vector (in Gauss) \[in COMMON
/WRKVLU/\]

BR Value of the B(r) magnetic field vector (in Gauss) \[in COMMON
/WRKVLU/\]

BT Value of the B() magnetic field vector (in Gauss) \[in COMMON
/WRKVLU/\]

EOMC Electronic charge divided by relativistic mass energy equivalent
(mixed units)

\[in COMMON /WRKVLU/\]

ERAD Average radius of the earth in kilometers \[in COMMOM /WRKVLU/\]

F(6) Array of \"F\" values (velocity and acceleration in program
coordinates)

\[in COMMON /WRKVLU/\]

SQY6 Intermediate term (Y6)\*Y(6)/Y(1)

TAY2 Intermediate term TSY2/TSC2

TCY2 Trigonometric cosine of the vector theta angle in r, , coordinates

\[in COMMON /SNGLR/\]

TCY3 Trigonometric cosine of the vector phi angle in r, , coordinates

\[in COMMON /SNGLR/\]

TSY2 Trigonometric sine of the vector theta () angle in r, , coordinates

\[in COMMON /SNGLR/\]

TSY3 Trigonometric sine of the vector phi () angle in r, , coordinates

\[in COMMON /SNGLR/\]

VEL Particle velocity in earth radii per second \[in COMMON /WRKVLU/\]

Y(6) Array of \"Y\" values (position and velocity in r, , coordinates)

\[in COMMON /WRKVLU/\]

Y5OY1 Intermediate term Y(5)/Y(1)

Subroutine MAGNEW95

This is the magnetic field evaluation program containing the IGRF95
magnetic field model. The normalized magnetic field coefficients have
been pre-processed and loaded in as data statements. This is a serial
computation designed to repeatedly evaluate the same magnetic field
model at different positions in space at optimum efficiency. Each term
of the magnetic field expansion is the derivative of the pervious term.

[Arguments in call statement]{.underline}: none. (All in labeled common)

[Labeled Common arguments]{.underline}:

Block name: /WRKVLU/

Arguments in block F(6), Y(6), ERAD, EOMC, VEL, BR, BT, BP

F(6) Array of force and acceleration vectors

Y(6) Array of position and velocity vectors

ERAD Average radius of the earth in kilometers

EOMC Electronic charge divided by relativistic mass energy equivalent
(mixed units)

VEL Particle velocity in earth radii per second

BP Value of the B() magnetic field vector (in units of Gauss)

BR Value of the B(r) magnetic field vector (in units of Gauss)

BT Value of the B() magnetic field vector (in units of Gauss)

Block name: /WRKTSC/

Arguments in block TSY2, TCY2, TSY3, TCY3

TSY2 Sine of the Y(2) coordinate (theta coordinate)

TCY2 Cosine of the Y(2) coordinate (theta coordinate)

TSY3 Sine of the Y(3) coordinate (phi coordinate)

TCY3 Cosine of the Y(3) coordinate (phi coordinate)

[Dimensioned variables:]{.underline} (not in labeled common)

G(11,11), BM(11)

G(11,11) Normalized coefficients ordered for fast serial computation

BM(11) Values to determine if the expansion should terminated at order N

[Program Operation]{.underline}:

This subroutine is designed for efficient serial computation of the
earth\'s main magnetic field. This procedure expands the terms into
FORTRAN coding (resulting in pages and pages of FORTRAN code) and then
evaluates the normalized field coefficients. The result of this serial
expansion is approximately an order of magnitude speed increase over the
recursion method which is much more compact in program size but requires
the expansion of the Legender polynomials each time the subroutine is
called.

When this subroutine is called it checks the value of the variable
JDATA. If JDATA is not 77, then it loads in the data coefficients;
otherwise it proceeds directly to the magnetic field evaluation. The
array G(11,11) contains the pre-processed and normalized magnetic field
model coefficients ordered for fast serial computation. The array BM(11)
contains a set of values that determine if the expansion should be
terminated above a specified order because the additional contribution
of the magnetic field would not be significant. (This is an additional
technique to speed up the computation.)

At the beginning of the magnetic field determination, the sine and
cosine of Y(3) (the phi coordinate) are calculated. (This has not been
done since the completion of the last Runge-Kutta step.) It then
determines the value of the AR variable; AR is reciprocal of the radial
distance, (Y(1)) ,in earth radii. Each order of the field expansion
requires an evaluation of AR^N^ where N is the order of the field
expansion. (The N = 1 term, a description of a monopole magnetic field,
is zero.) The first magnetic field evaluation designated by N = 2 is the
dipole field component. Each subsequent order of expansion evaluates the
contribution of the next order and adds this to the contribution of the
previous orders. The final computation converts the magnetic field
vectors to units of Gauss. ( One Gauss - 10^5^ Nt). A more detailed
description of the process is given in NSSDC DATA USERS NOTE 68-11.

[Listing of all variables]{.underline} in subroutine MAGNEW95

AOR Temporary value of radial distance (in earth radii) to the N^th^
power

AR 1.0/Radial distance (in earth radii) \[AR = 1.0/Y(1)\]

B Magnitude of the magnetic field (in Gauss) \[in COMMON /WRKVLU/\]

BERR Variable for dropping terms if field magnitude is less than
significant value

(set to 0.001 for improved accuracy)

BM(11) Array of check values for dropping terms of magnetic field
expansion

BP Value of the B() magnetic field vector (in Gauss) \[in COMMON
/WRKVLU/\]

BR Value of the B(r) magnetic field vector (in Gauss) \[in COMMON
/WRKVLU/\]

BT Value of the B() magnetic field vector (in Gauss) \[in COMMON
/WRKVLU/\]

DP101 Derivative of the polynomial (10,1) term

DP102 Derivative of the polynomial (10,2) term

DP103 Derivative of the polynomial (10,3) term

DP104 Derivative of the polynomial (10,4) term

DP105 Derivative of the polynomial (10,5) term

DP106 Derivative of the polynomial (10,6) term

DP107 Derivative of the polynomial (10,7) term

DP108 Derivative of the polynomial (10,8) term

DP109 Derivative of the polynomial (10,9) term

DP111 Derivative of the polynomial (11,1) term

DP1110 Derivative of the polynomial (11,10) term

DP1111 Derivative of the polynomial (11,11) term

DP112 Derivative of the polynomial (11,2 term

DP113 Derivative of the polynomial (11,3 term

DP114 Derivative of the polynomial (11,4) term

DP115 Derivative of the polynomial (11,5) term

DP116 Derivative of the polynomial (11,6) term

DP117 Derivative of the polynomial (11,7) term

DP118 Derivative of the polynomial (11,8) term

DP119 Derivative of the polynomial (11,9) term

DP21 Derivative of the polynomial (2,1) term

DP22 Derivative of the polynomial (2,2) term

DP31 Derivative of the polynomial (3,1) term

DP32 Derivative of the polynomial (3,2) term

DP33 Derivative of the polynomial (3,3) term

DP41 Derivative of the polynomial (4,1) term

DP42 Derivative of the polynomial (4,2) term

DP43 Derivative of the polynomial (4,3) term

DP44 Derivative of the polynomial (4,4) term

DP51 Derivative of the polynomial (5,1) term

DP52 Derivative of the polynomial (5,2) term

DP53 Derivative of the polynomial (5,3) term

DP54 Derivative of the polynomial (5,4) term

DP55 Derivative of the polynomial (5,5) term

DP61 Derivative of the polynomial (6,1) term

DP62 Derivative of the polynomial (6,2) term

DP63 Derivative of the polynomial (6,3) term

DP64 Derivative of the polynomial (6,4) term

DP65 Derivative of the polynomial (6,5) term

DP66 Derivative of the polynomial (6,6) term

DP71 Derivative of the polynomial (7,1) term

DP72 Derivative of the polynomial (7,2) term

DP73 Derivative of the polynomial (7,3) term

DP74 Derivative of the polynomial (7,4) term

DP75 Derivative of the polynomial (7,5) term

DP76 Derivative of the polynomial (7,6) term

DP77 Derivative of the polynomial (7,7) term

DP81 Derivative of the polynomial (8,1) term

DP82 Derivative of the polynomial (8,2) term

DP83 Derivative of the polynomial (8,3) term

DP84 Derivative of the polynomial (8,4) term

DP85 Derivative of the polynomial (8,5) term

DP86 Derivative of the polynomial (8,6) term

DP87 Derivative of the polynomial (8,7) term

DP88 Derivative of the polynomial (8,8) term

DP91 Derivative of the polynomial (9,1) term

DP92 Derivative of the polynomial (9,2) term

DP93 Derivative of the polynomial (9,3) term

DP94 Derivative of the polynomial (9,4) term

DP95 Derivative of the polynomial (9,5) term

DP96 Derivative of the polynomial (9,6) term

DP97 Derivative of the polynomial (9,7) term

DP98 Derivative of the polynomial (9,8) term

DP99 Derivative of the polynomial (9,9) term

DP1010 Derivative of the polynomial (10,10) term

EOMC Electronic charge divided by relativistic mass energy equivalent
(mixed units)

(not used in this subroutine) \[in COMMON /WRKVLU/\]

ERAD Average radius of the earth in kilometers (not used in this
subroutine)

\[in COMMOM /WRKVLU/\]

ERR Intermediate term used to evaluate when higher order terms can be
dropped

F(6) Array of \"F\" values (velocity and acceleration in program
coordinates)

\[in COMMON /WRKVLU/\]

G(11,11) Array of normalized magnetic field coefficients ordered for
fast computation

GMSUM Check value to assure that proper magnetic field coefficients are
loaded

GSUM Check data for testing that proper magnetic field coefficients are
loaded

JDATA Integer test value for loading in data statement for magnetic
field coefficients

JMAG Integer order of magnetic field expansion (N+1)

L Intermediate index for data checking

M Intermediate index for data checking

MGNMAX Integer value of maximum order of magnetic field expansion

P101 Polynomial (10,1) term

P1010 Polynomial (10,10) term

P102 Polynomial (10,2) term

P103 Polynomial (10,3) term

P104 Polynomial (10,4) term

P105 Polynomial (10,5) term

P106 Polynomial (10,6) term

P107 Polynomial (10,7) term

P108 Polynomial (10,8) term

P109 Polynomial (10,9) term

P111 Polynomial (11,1) term

P1110 Polynomial (11,10) term

P1111 Polynomial (11,11) term

P112 Polynomial (11,2 term

P113 Polynomial (11,3 term

P114 Polynomial (11,4) term

P115 Polynomial (11,5) term

P116 Polynomial (11,6) term

P117 Polynomial (11,7) term

P118 Polynomial (11,8) term

P119 Polynomial (11,9) term

P21 Polynomial (2,1) term

P22 Polynomial (2,2) term

P31 Polynomial (3,1) term

P32 Polynomial (3,2) term

P33 Polynomial (3,3) term

P41 Polynomial (4,1) term

P42 Polynomial (4,2) term

P43 Polynomial (4,3) term

P44 Polynomial (4,4) term

P51 Polynomial (5,1) term

P52 Polynomial (5,2) term

P53 Polynomial (5,3) term

P54 Polynomial (5,4) term

P55 Polynomial (5,5) term

P61 Polynomial (6,1) term

P62 Polynomial (6,2) term

P63 Polynomial (6,3) term

P64 Polynomial (6,4) term

P65 Polynomial (6,5) term

P66 Polynomial (6,6) term

P71 Polynomial (7,1) term

P72 Polynomial (7,2) term

P73 Polynomial (7,3) term

P74 Polynomial (7,4) term

P75 Polynomial (7,5) term

P76 Polynomial (7,6) term

P77 Polynomial (7,7) term

P81 Polynomial (8,1) term

P82 Polynomial (8,2) term

P83 Polynomial (8,3) term

P84 Polynomial (8,4) term

P85 Polynomial (8,5) term

P86 Polynomial (8,6) term

P87 Polynomial (8,7) term

P88 Polynomial (8,8) term

P91 Polynomial (9,1) term

P92 Polynomial (9,2) term

P93 Polynomial (9,3) term

P94 Polynomial (9,4) term

P95 Polynomial (9,5) term

P96 Polynomial (9,6) term

P97 Polynomial (9,7) term

P98 Polynomial (9,8) term

P99 Polynomial (9,9) term

RC10 Intermediate value (N=10) in magnetic field expansion

RC11 Intermediate value (N=11) in magnetic field expansion

RC2 Intermediate value (N=2) in magnetic field expansion

RC3 Intermediate value (N=3) in magnetic field expansion

RC4 Intermediate value (N=4) in magnetic field expansion

RC5 Intermediate value (N=5) in magnetic field expansion

RC6 Intermediate value (N=6) in magnetic field expansion

RC7 Intermediate value (N=7) in magnetic field expansion

RC8 Intermediate value (N=8) in magnetic field expansion

RC9 Intermediate value (N=9) in magnetic field expansion

TCP10 Trigonometric cosine function for the P10 term

TCP11 Trigonometric cosine function for the P11 term

TCP2 Trigonometric cosine function for the P2 term

TCP3 Trigonometric cosine function for the P3 term

TCP4 Trigonometric cosine function for the P4 term

TCP5 Trigonometric cosine function for the P5 term

TCP6 Trigonometric cosine function for the P6 term

TCP7 Trigonometric cosine function for the P7 term

TCP8 Trigonometric cosine function for the P8 term

TCP9 Trigonometric cosine function for the P9 term

TCY2 Trigonometric cosine of the vector theta angle in r, , coordinates

\[in COMMON /SNGLR/\]

TCY3 Trigonometric cosine of the vector phi angle in r, , coordinates

\[in COMMON /SNGLR/\]

TSP10 Trigonometric sine function for the P10 term

TSP11 Trigonometric sine function for the P11 term

TSP2 Trigonometric sine function for the P2 term

TSP3 Trigonometric sine function for the P3 term

TSP4 Trigonometric sine function for the P4 term

TSP5 Trigonometric sine function for the P5 term

TSP6 Trigonometric sine function for the P6 term

TSP7 Trigonometric sine function for the P7 term

TSP8 Trigonometric sine function for the P8 term

TSP9 Trigonometric sine function for the P9 term

TSY2 Trigonometric sine of the vector theta angle in r, , coordinates

\[in COMMON /SNGLR/\]

TSY3 Trigonometric sine of the vector phi angle in r, , coordinates

\[in COMMON /SNGLR/\]

VEL Particle velocity in earth radii per second (not used in this
subroutine)

> \[in COMMON /WRKVLU/\]

Y(6) Array of \"Y\" values (position and velocity in r, , coordinates)

> \[in COMMON /WRKVLU/\]

Reference publication:

[Computation of the Main Geomagnetic Field for Spherical Harmonic
Expansion,]{.underline} Data Users Note NSSDC 68-11, February 1968,
NSDC, NASA, GFSC, Greenbelt MD.