DESCRIPTION OF PROGRAM TJALLMAG

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 the NASA NSSDC ALLMAG program which has the option of
including a number of magnetic field models. This version has the
control (MODEL = 14) which specifies 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 TJALLMAG Main program; primary purpose is control

Subroutine GDGC Conversion from geodetic to geocentric coordinates

Subroutine SINGLTJ Calculates a particle trajectory

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

Subroutine ALLMAG Calculates 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.

[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 TJALLMAG

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 r, , coordinates

Y1GD Starting position, r component in geodetic r, , coordinates

Y2GC Starting position, component in geocentric r, , coordinates

Y2GD Starting position, component in geodetic r, , coordinates

Y3GC Starting position, component in geocentric r, , coordinates

Y3GD Starting position, component in geodetic r, , 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 TJALLMAG

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}: FGRADA

[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 FGRADA 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 FGRADA

This subroutine calculates the **V X B** force on the charged particle.
This is a slight modification of the standard FGRAD subroutine adapted
to call subroutine ALLMAG to calculate the magnetic field vectors. The
version of ALLMAG included in this software is an update from the
earlier version in the reference.

[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}:

ALLMAG (MODEL,TM,RKM,TSY2,TCY2,TSY3,TCY3,BR,BT,BP,B)

MODEL Integer designating the magnetic field model (14 = IGRF95)

TM Decimal value designating the epoch of the magnetic field

year (including decimal fraction of year)

RKM Radial distance of particle from geocenter in km.

TSY2 Sine of the theta () coordinate

TCY2 Cosine of the theta () coordinate

TSY3 Sine of the phi () coordinate

TCY3 Cosine of the phi () coordinate

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)

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

B Magnitude of the magnetic field (in Gauss)

[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. Then for compatibility with the standard ALLMAG call
statement the epoch of the magnetic field and the magnetic field model
are specified, the sine and cosine of the Y(2) and Y(3) coordinates (the
and coordinates) are calculated, and the radial distance (in km)
determined. The call to subroutine ALLMAG obtains the magnetic field
vectors (BR, BT, BP) at the particle position. Then the acceleration
vectors (F(4), F(5), F(6)) are calculated.

[Listing of all variables]{.underline} in subroutine FGRADA in program
TJALLMAG

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/\]

MODEL Integer designating the magnetic field model used in subroutine
ALLMAG

(needed for compatibility with standard call to subroutine ALLMAG)

RKM Radial distance of particle from geocenter in km.

(needed for compatibility with standard call to subroutine ALLMAG)

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/\]

TM Decimal value designating the epoch of the magnetic field

year (including decimal fraction of year)

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)

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 ALLMAG (MODEL,TM,RKM,ST,CT,SPH,CPH,BR,BT,BP,B)

This magnetic field subroutine is the NASA NSSDC ALLMAG software which
has the option of using a number of magnetic field models. This version
sets the variable MODEL to 14, which specifies the IGRF95 model of the
earth\'s magnetic field. The variable TM specifies the epoch of the
magnetic field expansion. *The FORTRAN coding in this subroutine is
unchanged from that received from NASA.* The use of this NASA version of
ALLMAG results in a slight deviation from the convention used in the
FORTRAN software coding used elsewhere in this software assembly. In the
previous subroutines, variables beginning with the character \"C\" are
used exclusively to denote character variables. In this subroutine
variables beginning with the character \"C\" are used, primarily to
denote cosine and co-latitude, and also intermediate terms in the
computation. In this subroutine the array \"G\"(13,13) has an
equivalence set of variables and can be addressed as either real or
integer.

When this subroutine is called, it first determines the model to be used
and the epoch of the model to be evaluated. It then derives the magnetic
field coefficients and normalizes the selected model coefficient for
fast computation. Then it loads the normalized and ordered coefficients
into the G arrays for magnetic field evaluation and initiates a serial
computation where each term is the derivative of the pervious term to
evaluate the magnetic field model at the specified position.

This technique is about 30 percent slower in execution speed than the
streamlined model specific code used in the TJI95 program, but is much
more versatile in the ability to select different magnetic field models
for any specified epoch.

[Arguments in call statement]{.underline}:

MODEL Integer designating the magnetic field model (14 = IGRF95)

TM Decimal value designating the epoch of the magnetic field

year (including decimal fraction of year)

RKM Radial distance of particle from geocenter in km.

ST Sine of the theta () coordinate

CT Cosine of the theta () coordinate

SPH Sine of the phi () coordinate

CPH Cosine of the phi () coordinate

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)

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

B Magnitude of the magnetic field (in Gauss)

[Labeled Common arguments]{.underline}: (designed to be used with other
software)

Block name: /TRAJAC/

Arguments in block CONSTEM, T, FILENAM

Block name: /DIPOLE/

Arguments in block WLONG, COLAT, EM

[Dimensioned variables:]{.underline}

TO(14), NMX(14), ISUM(14,3), G(13,13) LSUM(14,3),

TO(14) Array of the epoch years for the various magnetic field models

NMX(14) The order of the field expansion for each field model

ISUM(14,3) Check of the coefficients in data statement

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

LSUM(14,3) Check sums for the various models

Note that the following variables are REAL\*4:

GG(13,13,14), GGT(13,13,14), GGTT(13,13,14), SHMIT(13,13)

Note that the following variables are INTEGER\*4:

G1(13,13), GT1(13,13), GTT1(13,13)

G2(13,13), GT2(13,13), GTT2(13,13)

G3(13,13), GT3(13,13), GTT3(13,13)

G4(13,13), GT4(13,13), GTT4(13,13)

G5(13,13), GT5(13,13), GTT5(13,13)

G6(13,13), GT6(13,13), GTT6(13,13)

G7(13,13), GT7(13,13), GTT7(13,13)

G8(13,13), GT8(13,13), GTT8(13,13)

G9(13,13), GT9(13,13), GTT9(13,13)

G10(13,13), GT10(13,13), GTT10(13,13)

G11(13,13), GT11(13,13), GTT11(13,13)

G12(13,13), GT12(13,13), GTT12(13,13)

G13(13,13), GT13(13,13), GTT13(13,13)

G14(13,13), GT14(13,13), GTT14(13,13)

LG(13,13,14), LGT(13,13,14), LGTT(13,13,14)

Word of caution. We have found that on some of the old SGI machines, the
FORTRAN compile cannot accept the equivalence statement used in this
subroutine.

[Program Operation]{.underline}:

This subroutine is designed for 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
much faster than 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 first determines the model to be used
and the epoch of the model to be evaluated. It then derives the magnetic
field coefficients from the specified data statement and adjusts the
epoch year by the time derivative in the data statements and generates
normalized, ordered coefficients. (A byproduct of the magnetic field
evaluation is the determination of the dipole magnitude, the position of
the north magnetic dipole axis, and the offset of the dipole position
from the geocenter.) Then it loads the normalized and ordered
coefficients into the G arrays for magnetic field evaluation. It then
determines the value of the AR variable; AR is reciprocal of the radial
distance 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).

See the 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, for more details.

[Magnetic Field Models]{.underline}

[Model \#]{.underline} [Identification]{.underline}

1 HENDRICKS & CAIN GSFC 99-TERM 9/65 EPOCH 1960

2 CAIN ET AL. GSFC 120-TERM 12/6 EPOCH 1960

3 CAIN & LANGEL POGO 143-TERM 10/68 EPOCH 1960

4 CAIN + SWEENEY POGO 120-TERM 8/69 EPOCH 1960

5 IGRF 1965.0 80-TERM 10/68 EPOCH 1965

6 LEATON, MALIN & EVANS, 1965 80-TERM EPOCH 1965

7 HURWITZ (US COAST & GEODETIC) 168-TERM EPOCH 1970

8 IGRF 1980 168-TERM EPOCH 1980

9 IGRF 1975 80-TERM EPOCH 1975

10 BARRACLOUGH 168-TERM EPOCH 1975

11 AWC 168-TERM EPOCH 1975

12 IGRF 1985 168-TERM EPOCH 1985

13 IGRF 1990 168-TERM EPOCH 1990

14 IGRF 1995 168-TERM EPOCH 1995

[Listing of variables]{.underline} in subroutine ALLMAG

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

AR 1.0/Radial distance (in earth radii)

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/\]

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

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

C12 Intermediate value (N=12) in magnetic field expansion

C13 Intermediate value (N=13) in magnetic field expansion

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

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

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

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

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

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

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

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

COLAT Co-latitude of the north dipole axis

CONSTEM Magnitude of the dipole term

CP10 Trigonometric cosine function for the P10 term

CP11 Trigonometric cosine function for the P11 term

CP12 Trigonometric cosine function for the P12 term

CP13 Trigonometric cosine function for the P13 term

CP2 Trigonometric cosine function for the P2 term

CP3 Trigonometric cosine function for the P3 term

CP4 Trigonometric cosine function for the P4 term

CP5 Trigonometric cosine function for the P5 term

CP6 Trigonometric cosine function for the P6 term

CP7 Trigonometric cosine function for the P7 term

CP8 Trigonometric cosine function for the P8 term

CP9 Trigonometric cosine function for the P9 term

CP Cosine of the phi coordinate

CT Cosine of the theta coordinate

DP101 Derivative of the polynomial (10,1) term

DP1010 Derivative of the polynomial (10,10) 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

FILENAM Externally provided identification

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

I Index in do loops

JJ Integer control parameter for loading coefficients

K Index in do loops

L Intermediate index for data checking

M Intermediate index for data checking

MODOLD Integer control variable for model changing

N Index in specifying degree and 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

SHMIT Array of coefficients from model data statements

SP10 Trigonometric sine function for the P10 term

SP11 Trigonometric sine function for the P11 term

SP2 Trigonometric sine function for the P2 term

SP3 Trigonometric sine function for the P3 term

SP4 Trigonometric sine function for the P4 term

SP5 Trigonometric sine function for the P5 term

SP6 Trigonometric sine function for the P6 term

SP7 Trigonometric sine function for the P7 term

SP8 Trigonometric sine function for the P8 term

SP9 Trigonometric sine function for the P9 term

SPH Sine of the phi coordinate

ST Sine of the theta coordinate

T Difference in decimal years between specified epoch and model
reference year

TO Reference year for the specified magnetic field model

TM Time (decimal year) for the selected magnetic field model

to be evaluated.

TMOLD Year of the previous magnetic field model used

WLONG West longitude of dipole axis

Reference publications:

[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.

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