m.e.ash   june 1966    subroutine bodred
control and data constants are initialized and read for
earth, moon, planets, asteroids, satellites, artificial space
probes, earth rotation, moon rotation and planet rotation.

--- harmonic+shape storage ---
for nmlst2 with positive nplnt, gravitational
potential harmonics is assumed.
if a -nplnt nmlst2 is input, this is assumed
to be planet shape. model depends on nshp.
nshp=0,(default) spherical harmonic expansion
nshp=1,          fourier series expansion
nshp=2,          altitude grid (local model)
nshp=3,          external model expected to be supplied on obslibs
(nothing for nshp .gt.3 yet)
note: see later comments on shape for full
explanation.

ovely(lovely)equivalenced starting with j2(lj2)
for the next 252 real*8(integer*2)elements.
*** historical reasons for 252 elements. doesn't hurt. ***
fa thru fbpp variables are the corresponding
fourier coefficients for the two dimensional fourier
series representation (note:  read p as primed and
pp as double primed). lfa thru lfbpp are the
corresponding l-vectors. see later comments on
shape for full explanation.

imnk,idyk,iyrk = month,day,year converted to jdk (k=0,1,2)
if iyrk.gt.0
year is 0 to 99 in 1900s
year is greater than or equal to 100 in 2000s

sclzon = 0 unscaled zonal harmonics are input (default)
sclzon = 1 scaled zonal harmonics are input (change to unscaled
before return to calling program)
scltes = 0 unscaled tesseral harmonics are input (change to scaled
before return to calling program)
scltes = 1 scaled tesseral harmonics are input (default)

NPLNT = body number
1 = Mercury
2 = Venus
3 = Earth-Moon barycenter
4 = Mars
5 = Jupiter
6 = Saturn
7 = Uranus
8 = Neptune
9 = Pluto
10 = Moon
11,...,30 = other possible natural planets,asteroids or satellites
31,32,... = artificial space probes
-1,-2,-3,-4,...,-9,-10,-11,... = rotation for planet iabs(NPLNT)
0 = indicates end of sequence of &nmlst2 namelists,
subroutine input gets out of loop calling subroutine
bodred. in other parts of program, planet number 0 is
associated with the sun.
--- However, if SMALL=t, then NPLNT is an asteroid number,
starting with 1 for Ceres, 2 for Juno, etc.  Asteroids may be
input either as ordinary minor planets with NPLNT=11 to 30 or
as limited asteroids with any NPLNT number.

NCENTR = central body of numerical integration
-1,0 = Sun
1 = Mercury
2 = Venus
3 = Earth (not Earth-Moon barycenter)
4,5,...same as above

if NPLNT between 1 and 30 and NCENTR=0, individual body integration
with sun as central body is done in subroutine FN
if NPLNT between 1 and 30 and NCENTR.lt.0, individual body integration
with sun as central body is done in subroutine SBFN
if NPLNT.gt.30 or NCENTR.gt.0, individual body integration with
central body NCENTR is done in subroutine SBFN
if NPLNT=10, moon integration is always done in MONFN
if NPLNT= 3, earth-moon barycenter integration is done in FN or SBFN
depending on whether ICT(40)=0 or 1
There is no indivdual body integration if JD0=0 or if the specific
body is one of the NPLBDY(.) bodies in the n-body /BDCTRL/
labeled common

(mass of probe NPLNT)/(mass of central body NCENTR)=0 if
NPLNT.gt.30
if NCENTR.le.0, (mass of body NPLNT)/(mass of sun)=MASS(NPLNT) if
NPLNT between 1 and 30 (MASS vector is in /PARAM/
labeled common)
if NCENTR.gt.0, (mass of body NPLNT)/(mass of planet-satellite system
NCENTR) =MASS(NPLNT) if NPLNT between 1 and 30
Note that when we integrate the motion of Jupiter, for instance, we
are actually integrating for the motion of the center of mass of the
Jovian planet-satellite system relative to the sun and that mass(5) is
the total mass of the Jovian planet-satellite system divided by the
mass of the sun

JD0 =initial epoch of NPLNT integration is midnight beginning of
day of day with julian day number   JD0  (if negative there is
checkpoint restart at the epoch   -JD0  unless
JD0 .le.-3000000. if JD0=-1 there is checkpoint restart
at record just before end of file on ephemeris tape.)
(if JD0=0 or JD0=-3000000, there is no
numerical integration, comparison of theory and observation
assumes motion is already on tape, but if initial conditions
are adjusted,   JD0  is set to positive value on second record
of tape so there will be reintegration on next least squares
iteration, whereas if no initial conditions are adjusted,
JD0  is set to zero in the compar link so there will be no
reintegration on subsequent least squares iterations)
(  JD0 .le.-3000000 is signal to compar link on first
iteration to override end time on tape by input   JD2 )

Initial time is midnight ephemeris time on day with Julian day
number JD0 (Julian date at noon on given day), except that
integrations in SBFN (NPLNT.gt.30 or NCENTR.ne.0) can have
initial epoch at day JD0 and nonzero fraction of a day
INT1*(2**INT2). INT2 negative and less than 30 in absolute
value. INT1 positive and less than 2**iabs(INT2).  Initial
epoch time can be input as IHR,IMIN,SEC ephemeris time or as
FRACT0 fraction of day, rather than INT1,INT2, but it will be
converted to INT1 with INT2=-24 (accurate to .01 second)

JD1 =julian day number of first data record on output tape
JD2  julian day number of last data record on output tape
(see also eps(1) and eps(2) below)
see also description of jdpad in PRMRED

must have JD0 lying between JD1 and JD2, unless JD0 is
negative indicating checkpoint restart, in which case   -JD0
must be between JD1 and JD2, or unless JD0 is zero
indicating no integration or less than or equal to -3000000
indicating no integration and input JD2 to override that on
second record of integration data set in COMPAR link (needed
because such data set might have been produced by check point
restart going beyound point indicated on second record)
(if JD0=-1, consistency check for -JD0 between JD1,JD2 not made

if JD0 .lt.0 & .gt.-3000000, output tape from a previous integration
is read to the time -JD0 (or to the record just before end
of file if JD0=-1 ) and the coordinates on the tape at
that epoch are used to restart that integration, proceeding
from -JD0 to JD2. Integration cannot go in both
directions from epoch -JD0 in restart mode and a
previous integration cannot be restarted at a time in the
first direction of the previous integration if it was in the
two direction mode, since output tape is not written until it
gets to second direction. In the checkpoint restart mode with
JD0 negative, the only data which can be input for the
integration are JD1,JD0,JD2,K(99),EPS(2),EPS(6) because rest
of the integration data is taken from the second record of the
tape of the previous integration.  If K(99) .lt.0, the tape
of the previous integration is rewound, whereas if
K(99) .ge.0, the restarted integration continues writing
the tape of the previous integration.

if JD0 is positive between JD1 and JD2, we have
(a) if JD0=JD1, numerical integration starts at JD0
and proceeds to JD2, writing output tape as integration
proceeds if K(99).ge.0
(b) if JD0 is strictly between JD1 and JD2, integration
starts at JD0 and goes in direction of JD1 (forward
in time if JD0 .lt. JD1, backward in time if JD0 .gt.
JD1). The integration polynomial coefficients determined
at JD0 by the starting proceedure are saved and the
numerical integration proceeds towards JD1, writing onto
disk buffer IBUF if K(99).ge.0.  When JD1 is reached,
IBUF is backspaced and read to write output tape ITAPE in
time direction from JD1 to JD0 if K(99).ge.0.  Then
using the saved starting polynomial coefficients, the
integration goes from JD0 to JD2, writing output
tape ITAPE as the integration proceeds if K(99).ge.0.
When the integration is completed, we have an output tape
from JD1 to JD2, even though integration went from
JD0 to JD1 and then from JD0 to JD2. These
contortions are gone through in order to simplify the logic
needed to read the tape in the COMPAR link.

Note: in a simultaneous integration of lunar orbit,
libration, and core rotation, three separate buffers are
needed for the two-direction mode.  IBUF is used for the
orbit, while KK(13) and KK(12) of the Moon rotation are
used for the libration and core rotation, respectively.

INT  gives interval between tabular points.
if  INT.gt.0, tabular interval is  INT days (only value allowed
if NPLNT=1,...,9,11,...,30 and NCENTR=0)
if  INT.le.0, tabular interval is 2**INT days (for the Moon, NPLNT=10,
it is expected that INT=-1 for half day tabular intervl)
These constant output tabular intervals apply throughout the
integration for integration methods K(88)=1,2,3. for K(88)=0 method,
output interval from initial epoch to first tabular point is
controlled by INT. Thereafter, tabular output is every step of
variable step Nordsieck integration. Output every step immediately
after the initial epoch is supressed because the starting interval for
the Nordsieck method of integration is very small, but it grows to the
size warranted by the stability requirements of the differential
equations integrated and the input accuracy constant EPS(3) as the
integration proceeds away from the initial epoch.

ITAPE =   data set number for output tape
first record is 88 characters describing computer run,second
record is data about body integration, including
IPAR=number of partials+1 (calculated in setup for integ.)
third and subsequent records have integration results.

(1) if NPLNT=1,...,9,11,...,30 with NCENTR=0 there are five
tabular points per record with format
JD,FRACT,IVL,(((CRD(i,j,k),i=1,IVL),j=1,IPAR),k=1,5)
where JD is Julian day number of the data at the first
tabular point,FRACT is double precision fraction of day from
midnight beginning of day for this JD (0 if INT.gt.0)
and where for k=1,5 the double precision CRD array is
CRD(i,1,k) = position,velocity at tabular point k, i=1,6
CRD(i,j,k) = partial derivative of position,velocity with
respect parameter j at tabular point k,
i=1,6 and j=2,IPAR
if  IPAR=1, no partials on tape
IPAR  calculated in integration setup routine from KI(1-n)
the meaning of the partials j=2,ipar  is given by  KI(1-n).
IVL=6 in integration output routine so velocity will be on
tape for checkpoint restart purposes. To save space on tape
or disk output can be copied with IVL=3 before recent times
since optical observations do not need accurate velocity
and the velocity they need is gotten by numerical
differentiation from position tabulation.
Everett eighth difference interpolation is used to
interpolate from this tape in the COMPAR link. There are
three records in storage in the COMPAR link and there are
five tabular points per record so that can interpolate for
any time in the middle record (need more than four tabular
points per record because of the existence of receive,
reflection and send times of an observation)

(2) if NPLNT=10 for the Moon there are eight tabular points
per record with format
JD,FRACT,IVL,(((CRD(i,j,k),i=1,IVL),j=1,IPAR),k=1,8),
nutation,libration,k=1,8)
where the meaning of JD,FRACT,IVL,CRD is the same as in
case(1) with FRACT =0 if INT=-1,-2,-3, where the
nutation of the Earth and the physical libration of the
Moon are copied from the perturbing planet data set for
use in the COMPAR link.

(3) if NPLNT=1,...,9,11,...,30 with NCENTR=-1 or
NCENTR.gt.0, or if NPLNT.gt.30 there are two cases.
if K(88).gt.0, it is the same as case(1) above.
if K(88)=0 (Nordsieck variable output interval) there is
one tabular point per record with format
JD,FRACT,IORD,JORD,NUMPRT,NUMPT1,((SPRB(i,j),i=1,3),
j=1,IORD),(((DSPRB(i,j,k),i=1,3),j=1,JORD),
k=1,NUMPT1),HC1,HC2
where meaning of JD,FRACT is the same as above and where
for the given tabular point at time JD,FRACT
IORD  = order of time derivative - 1 for motion on tape
JORD  = order of time derivative - 1 for partial on tape
NUMPRT= number of partials on tape
NUMPT1=NUMPRT if NUMPRT.gt.0, 1 otherwise
SPRB(i,1)= postion coordinate i, i=1,2,3
SPRB(i,2)= velocity coordinate
SPRB(i,3)= acceleration coordinate
SPRB(i,4)= jerk coordinate
etc. up to j=IORD
DSPRB(i,j,k)= similar meaning for partial of motion with
respect to parameter k
HC1   = time distance from previous tabular point to
this tabular point
HC2   = time distance from this tabular point to next
next tabular point
There are six records in storage in compar link for Hermite
interpolation for a time point in the midst of the middle
four tabular points so that either a two or four point
Hermite interpolation method can be used of order limited
only by the number of derivatives on the tape and by our
being able to derive the formulas. The derivatives after
acceleration come from the Nordsieck intgration coefficients
which are not the real higher derivatives in so far as
they tend to lag behind the current tabular point, the jerk
lagging least of all.

COND(i), i=1,6,= initial conditions for equations of motion
(a) Normally, COND(1-6) are the
initial osculating eliptic orbital elements with non-time
variable gravitational constant and with augmented mass of
central body (so that equations for partial of motion with
respect to given body mass will have non-zero init.cond.)
COND(1)=A    =semi-major axis (astronomical units)
COND(2)=E    =eccentricity
COND(3)=INC  =inclination (deg)
COND(4)=ASC  =right ascension of ascending node (deg)
COND(5)=PER  =argument of perihelion   (deg)
COND(6)=ANOM =initial mean anomaly (deg)
(b) If NPLNT.lt.0 then COND(1-6) are initial Euler angles and
rates (rad and rad/day).
(c) Input controls may specify other types of initial conditions.
The equations of motion are integrated in a coordinate system
referred to the mean equinox and equator of the reference epoch.
(See JCT(13).)

The above conventions for initial conditions are superceded by the
following (default value of icnd is given by above)
ICND =-1 initial conditions are cartesian position and velocity
(units are astronomical units and astronomical units per day)
(referred to mean equinox and equator of reference epoch)
ICND = 0 initial conditions are elliptic orbital elements
ICND = 1 initial conditions are elliptic orbital elements with sums of
angles   ASC, ASC+PER, ASC+PER+ANOM
ICND = 2 initial conditions are
COND(1) = distance from central body (astronomical units)
COND(2) = right ascension of position (epoch coordinates)
(between 0 and 360 deg)
COND(3) = declination of body (epoch coordinates)
(between -90 and 90 deg)
COND(4) = magnitude of velocity minus the velocity in
circular orbit at radius COND(1) (au/day)
COND(5) = azimuth of velocity vector on plane normal to
radius vector (between 0 and 360 deg)
COND(6) = flight path angle measured from radius vector to
velocity vector (between 0 and 180 deg)
ADBARV
a = right ascension = COND(2)
d = declination     = COND(3)
b = beta            = COND(5)
a = angle           = COND(6)
r = radius          = COND(1)
v = velocity        = COND(4)+circular orbit velocity
ADBARV initial conditions can be read in &nmlst2 namelist by the names
ALPHA,DELTA,BETA,ANGLE,R,V with appropriate values for ICND,KCND,LCND
and also JTYPE

The initial conditions must be in one of the above forms in the input
stream, but several variants are allowed, namely, distance and vel-
ocity may be in other units, the reference frame may be other than
1950 mean equatorial, and the velocity 'v' of adbarv notation may be
the total velocity.  The variants, if any, are all resolved to the
standard forms in the input procedure.  For convenience, the elements
may be input in one form (say, Cartesian) and converted to another
(say, elliptic).  The desired final form is always specified via
'ICND' (as above), and the input form, if different, via 'INCND' with
(with the same logic as 'ICND' above).

INCND=-1 input elements are Cartesian
INCND= 0 input elements are elliptic
INCND= 1 input elements are elliptic with angle sums
INCND= 2 input elements are ADBARV

INUNIT=-1 input distances, velocities are au, au/40 d
INUNIT= 0 input distances, velocities are au, au/d
INUNIT= 1 input distances, velocities are km, km/s
INUNIT= 2 input distances, velocities are ft, ft/s

For Moon rotation there is special logic that applies only to the
core's angle rates:
INCND=0 input rates are Euler angle rates
INCND=3 input rates are angular velocity components of the core in
the frame of the mantle (to be converted to "0" if ICND=0)
If ICND=3, then retain these values and integrate only the
core angular velocity, ignoring what would otherwise be the
Euler angles

* * * * * obsolete * * * * * see 'INUNIT' and 'INCND' above
ICND switches internal to BODRED are
ICND =-8 change input position and velocity in km and km/sec to
the ICND=1 initial conditions
ICND =-7 change input position and velocity in km and km/sec to
the ICND=0 initial conditions
ICND =-6 change input position and velocity in au and au/day to
the ICND=1 initial conditions
ICND =-5 change input position and velocity in au and au/day to
the ICND=0 initial conditions
ICND =-4 same as ICND=1 but unit of COND(1) is kilometers instead of
au.  COND(1) is changed to astronomical units (au) and icnd
is set equal to 1
ICND =-3 same as ICND=0 but unit of COND(1) is kilometers instead of
au.  COND(1) is changed to astronomical units (au) and icnd
is set equal to 0
ICND =-2 same as INCD=-1 but units of position and velocity are
kilometers and kilometers per second instead of au and au/day
COND(1-6) are transformed to au and au/day and ICND is set
equal to -1
ICND = 3 same as ICND=2 but units for COND(1) & COND(4) are kilometers
and km/sec instead of au and au/day, respectively.
COND(1) is transformed to au and COND(4) is transformed to
au/day and icnd is set equal to 2
ICND = 4 change input position and velocity in au and au/day to
the ICND=2 initial conditions in au and au/day and set
ICND equal to 2
ICND = 5 change input position and velocity in km and km/sec to
the ICND=2 initial conditions in au and au/day and set
ICND equal to 2
* * * * * obsolete * * * * *

* * * * * obsolete * * * * * see 'INUNIT' and 'INCND' above
JCND switches internal to BODRED for preprocessing of initial
conditions before ICND preprocessing is done
JCND = 0 no element preprocessing in subroutine CHNCNE
JCND = 1 change initial osculating orbital elements (ICND=0 type)
to Cartesian coordinates (semi-major axis unit is
astronomical units)
JCND = 2 same as JCND=1 but semi-major axis unit is kilometers
(change to astronomical units)
* * * * * obsolete * * * * *

* * * * * obsolete * * * * * see 'INUNIT' and 'INCND' above and JTYPE
KCND switches internal to BODRED for preprocessing of initial
conditions before JCND,ICND preprocessing is done
KCND = 0 no coordinate system preprocessing in CHNCNE
KCND = 1 change initial osculating orbital elements (ICND=0 type)
from Euler angles inc,asc,per referred to true
equinox and equator of date to mean equinox and equator
of reference epoch
KCND = 2 change initial cartesian coordinates referred to true equinox
and equator of date to mean reference equinox and equator
KCND = 3 change initial ADBARV coordinates referred to true equinox
and equator of date to initial cartesian coordinates referred
to the mean reference equinox and equator (ICND.gt.2 can
change back to ADBARV)
(units are au & au/day)
KCND = 4 units changed from km & km/sec and then KCND=3 is done
* * * * * obsolete * * * * *

JTYPE overrides above change to mean equinox and equator of epoch.
Instead of true equinox and equator of date, original coordinates are
JTYPE =-1 mean equinox and equator of 1950.0
JTYPE = 0 mean equinox and equator of date
JTYPE = 1 true equinox and equator of date (default assumed if
obsolete KCND is given)
JTYPE = 2 mean equinox and ecliptic of date
JTYPE = 3 mean equinox and ecliptic of 1950.0
JTYPE = 4 mean lunar plane of date (x axis along intersection of mean
lunar and ecliptic planes of date)
JTYPE = 5 mean equinox and true equator of date JDTYPE (0 hr)
JTYPE = 6 mean equinox and equator of date JDTYPE (12 hr)
JTYPE = 7 mean equinox and ecliptic of date JDTYPE (12 hr)

LCND = 0 nothing
LCND = 1 input ADBARV COND(4) is velocity and not difference between
velocity and circular orbit velocity  (units au & au/day)
change to difference velocity
LCND = 2 change from km & km/sec to au & au/day and then do LCND=1
(obsolete, see 'INUNIT' above)

* * * * * obsolete * * * * * see 'INUNIT' and 'INCND' above
Before any of the above are done, the following is done.
IFTKM = 0 nothing
IFTKM = 1 for elliptic orbital element initial condtions change semi-
major axis from feet to kilometers
IFTKM = 2 for ADBARV initial conditions change radius COND(1) from
feet to kilometers and velocity COND(4) from ft/sec to
km/sec
IFTKM = 3 for cartesian initial conditions change position and
velocity COND(1-6) from ft and ft/sec to km and km/sec
* * * * * obsolete * * * * *

MEAN = 0 just the usual &nmlst2 namelist read for this body
MEAN = 1 read the NORAD 2 card mean element set after this &nmlst2.
The mean elements are changed to Cartesian coordinates, and
set JTYPE=0, INUNIT=0, INCND=-1.  ICND should have been set
to the desired output type of orbital elements (in this case
the default is -1).

MCENTR =-2 NCENTR gives the central body, no changes made
MCENTR.ge.-1 NCENTR gives the central body that the initial conditions
are for. Change initial conditions to the MCENTR central
body and set NCENTR=MCENTR. Before this is done, by use of
ICND,JCND,KCND,LCND,JTYPE we must have ended up at the
ICND=-1 type of initial conditions (Cartesian coordinates
referred to the mean reference equinox and equator in au
and au/day). After the change of central body, the initial
conditions can be changed to the ICND=0,1 or 2 type by
MCND switch, which does not apply unless MCENTR.ge.0
the preferred values of MCND are (-1,0,1,2), but the
following (obsolete) values from the (obsolete) ICND
conventions can also be used:
MCND =-6 change ICND=-1 initial conditions to ICND=1 i.c.
MCND =-5 change ICND=-1 initial conditions to ICND=0 i.c.
MCND = 4 change ICND=-1 initial conditions to ICND=2 i.c.

For a planet, asteroid or natural satellite (NPLNT.le.30 and
NPLNT not equal 3 or 10 nor negative)
CON(1)=radius of body
CON(i), i=1,12,  adjustable parameters giving shape and rotation
of planet
CON(2),CON(3),CON(4) = j1,c11,s11 for spherical harmonic shape
model
CON(2)               = flattening for fourier series shape model
one of CON(2),CON(3),CON(4) must be non-zero in order to
include planetary shape in processing radar observations of
planets.  shape coefficients are specified as harmonic input
for body -NPLNT.
CON(6)  = phase angle of rotation at epoch CON1(1) (this has
different meanings in ellipsoidal shape routines in
compar link and in the calculation of the effect of this
body's gravitational potential harmonics on another
body in subroutines sbfn,sbfn1)
for special mars rot. code, this angle is referred to
the inertial frame defined by CON1(7) and CON1(8)
CON(7)  = rotation period in days at epoch CON1(1)
CON(8)  = declination of rotation pole at epoch CON1(1) referred
to the mean equinox and equator of reference epoch
CON(9)  = right ascension of rotation pole at epoch CON1(1)
referred to the mean reference equinox and equator
CON(10) = i0, inclination at epoch CON1(1) referred to
intermediate inertial frame defined by CON1(7) and
CON1(8)
CON(11) = psi0, precession phase angle at epoch CON1(1)
referred to inertial frame defined by CON1(7) and
CON1(8)
CON(12) = mu, precession constant
CON(13) = time derivative of CON(7) in days per days past CON1(1)
CON(14) = a: CON(14) through CON(17)
CON(15) = b: are coeficients in the seasonal change of
CON(16) = c: spin period of mars. see memo by rdr dated 2/14/78
CON(17) = d: units are in radians
CON(14) = A1 nongrav. force coefficient for comets (see K(83))
CON(15) = A2 nongrav. force coefficient for comets
CON(18) = atmospheric scale height (km) of this planet (in patctl)
CON(19) = atmospheric zenith delay (s) at this planet's surface,
assuming a constant scale height
CON(20) = deltap = principle of equivalence violation parameter
eta*delta(body) where eta = 4*beta-gamma-3,
delta = ratio of gravitational to total energy (neg.)
see CON1(10).
CON(21) = atmospheric scale height (km) of this planet (in plnorb)
[Yes, this is conceptually the same as CON(18), but the
two are measured from completely different phenomena
and may need to be estimated independently.]
CON(22) = ref. density (gm/cm**3) of atmosphere: see CON1(12)
CON(23) = radius of body in kilometers for transit or occultation
observations
CON1(1) = epoch from which rotation is measured
for asteroids or satellites 11,...,30 the epoch is
also used for initial conditions for elliptic orbits
which perturb motion in integration in subroutines
SBFN,SBFN1.  Epoch is Julian ephemeris date.
CON1(2-4)= components of unit vector pointing from earth to given
body at epoch defining zero longitude on body.
CON1(5)  =constant (degrees) to be added to longitude calculated
with CON1(1-4) to get IAU system of longitude (applied
in subroutine deldop, affects shape model of planet)
CON1(6)  =constant (degrees) to be added in subroutine RADAR to
IAU longitude on data card for time delay observation
away from the subradar point to get longitude that would
be calculated with CON1(1-5)
(only one of CON1(5) and CON1(6) should be nonzero)
(we assume the external world is on the IAU east (not
west) longitude system and we either convert PEP to the
external world with CON1(5) or convert the external
world to PEP with CON1(6), presuming we were not able
to get the IAU sytsem with CON1(1-4) alone).
CON1(7) = n, an angle defining intermediate inertial frame.
Relates inertial frame defined by a mean orbit of date
to the reference inertial frame. See Reasenber and King,
"The Rotation of Mars" for more precise definition.
CON1(8) = j, same comments as for CON1(7)
CON1(9)=drgeps,if exponent of density function .gt.drgeps
then print exponent with warning: see CON(21-22)
CON1(9) could be converted to a TCON if the space is
needed for another purpose
CON1(10) = ratio of gravitational to total energy (negative,
for principle of equivalence calculations): see CON(20)
CON1(12)  ref. altitude for CON(22)  (km)

For all probes (NPLNT.gt.30) in the COMPAR link
CON(16)= transmitter frequency (hz)
CON(17)= transmitter frequency rate (hz/sec)
CON(18)= transmitter frequency accelerating (hz/sec**2)
CON1(3)= initial epoch for CON(17) and CON(18) (PEP jd.fract)

For a probe (NPLNT.gt.30) with logical center the sun
CON(i),i=1,11,  adjustable parameters giving low thrust and other
forces which affect motion of probe
gas leak force parameters
CON(1) = alpha 1
CON(2) = alpha 2
CON(3) = f 1
CON(4) = f 2
CON(5) = f 3
radiation pressure force parameters
CON(6) = g 1
CON(7) = g 2
CON(8) = g 3

For a probe (NPLNT.gt.30) with logical center a planet (low thrust
forces computed in plnorb)
CON(1)- CON(4) magnitudes of gas leak terms (dynes)
CON(5)- CON(8) inverse time constants for gas leaks (/day)
CON1(1) rt. asc. of canopus or other ref. star (deg.)
CON1(2) dec. of canopus or ref. star (deg.)
CON1(3) beta:rotation of gas leak coor. system from sun-star plane
CON1(4) epoch for time varying gas leak forces (j.e.d.)

radiation pressure parameters
CON(9)-CON(15) radiation pressure constants
CON1(5) beta2:rotation of rad. pres. system from sun-star plane
CON1(6) epoch for time varying radiation pressure quan. (j.e.d.)

atmospheric drag parameters
CON1(10)= cd1 drag coefficient of probe
CON1(11)=cd2 lift coefficient of probe
in COMPAR link we also have:
CON(23) transponder delay for time delay and counted doppler
calculation ((seconds)
CON(1) for instantaneous Doppler (NCODF=1), assumed radius
of observed body for time delay (km)

For a probe (NPLNT.gt.30) with logical center the Earth (low-thrust
forces computed in ERTORB)
radiation pressure parameters (see also KK(78))
CON(9) = scale factor for direct acceleration
CON(10)= scale factor (x total direct accel.) for acceleration
in y-axis direction (for GPS satellites with isotropic
or nds model)
CON(11)= angle between solar panels and y-axis (for GPS NDS model
only) (radians)

for LES-8 and LES-9 Earth satellites in subroutine LESTHP
CON1(1) = weight of satellite in pounds (if zero, default launch
values used in in computing radiation pressure and
thrust accelerations)
CON(17) = solvable parameter multiplying LES-8/9 thrust
(if zero, assumed one)
CON(18) = solvable parameter multiplying LES-8/9 radiation
pressure (if zero, assumed one)

for LES-8/9 station keeping simulation (setup at end of subroutine
SBSET, used in subroutines LESTHP and STATON)
CON1(2) = satellite orbit plane inclination to ecliptic (deg)
CON1(3) = angle between ascending node of satellite orbit
plane on ecliptic and sun perigee location (deg)
CON1(4) = desired station east longitude (deg)
CON1(5) = damping constant (dimensionless)
CON1(6) = coasting constant (deg/day)
CON1(7) = saturating drift constant (deg/day)
CON1(8) = station keeping thrust level (millipounds)
CON1(9) = bias to be added to east longitude sought by station
keeping system (to correct for discrepency between
longitude derived from sun transits and actual longitude
of center of LES-8/9 figure eight ground track)
CON1(10)= Julian ephemeris date of enabling of sun transit
measurements (half day less than jd,fract)

The elliptic orbit elements COND(1-6) are referred to the plane of
the sky, with the three coordinate axes pointing north, east, and
toward the observer.  For observations of the primary, the orbit
position (suitably scaled by the mass factor) must be subtracted.
The initial epoch of the orbit is given by JD0 of the primary.

CON(3) = orbital period in days, explicitly imposed by scaling the
"mass" used in the elliptic orbit formulae as the cube of
the input semimajor axis.
CON(4) = mass of body divided by the sum of masses of the body and
its primary.  (Not yet implemented.  The mass factor for
the body is assumed to be incorporated in COND(1).)

For any integration which uses SBFN (probe, satellite or even embary
if forced there by ICT(40)=1 or planet if forced their by NCENTR=-1)
CON1(12)= distance below which distance to target body and
if KK(81).gt.0 distance to central body is printed out
(unit is astronomical units)

For the Earth (NPLNT=3)
CON(1) = mean equatorial radius in kilometers
CON(2) = degree of flattening
CON(3) = Love number h for solid Earth tides
CON(4) = Love number l for solid Earth tides
CON(5) = lag time (sec) for solid Earth tides

For the Moon (NPLNT=10)
CON(1) = mean radius pointing towards mean subradar point
CON(3) = Love number h for solid Moon tides
CON(4) = Love number l for solid Moon tides
CON(5) = lag time (sec) for solid Moon tides
CON(6) = ad hoc node precession rate (rad/day about ecliptic)
CON(16-18) constants for effect of tidal friction on moon motion
sin(2*lag angle)=CON(16)+CON(17)*t+CON(18)*t**2
CON(19) = dels = principle of equivalence violation parameter
eta*(delta(Earth)+delta(Moon)) where
eta=4*beta-gamma-3 and delta (negative) is ratio of
gravitational to total energy
CON(20) = deld = eta*(delta(Earth)-delta(Moon))

For the Earth rotation (NPLNT=-3)
CON(1-21) are shared by two mutually exclusive sets of parameters
if JCT(29) = 0:
CON(1) = psidot(1) ad hoc small rotation angular rate
CON(2) = psidot(2)
CON(3) = psidot(3)
CON(4) = psi(1)    ad hoc small rotation angle
CON(5) = psi(2)
CON(6) = psi(3)
CON(11-13) polynomial terms in CT-UT (before 1956 Jan 17.0,
applied if JCT(34).gt.0) or in A1-UT1 (after 1956 Jan
17.0, applied if erotat CON1(1).gt.0).  See JCT(34)
CON(11) = r0  constant term for polynomial variation in A1-UT1
(CT-UT before 1956) (sec).  See A1UT1F or CTUTF
CON(12) = r1  linear term of this polynomial variation (sec/cy)
epoch for A1-UT1 is CON1(1), epoch for CT-UT is
1956 Jan 17.0.
CON(13) = r2  quadratic term of this polynomial variation
CON(14-21) seasonal terms in ut2-ut1 (before 1956 Jan 17.0) or
a1-ut1 (after 1956 Jan 17.0), applied if ICT(34).gt.0
CON(14)= coefficient of cosine annual term
CON(15)= time variation of CON(14)
CON(16)= coefficient of cosine semi-annual term
CON(17)= time variation of CON(16)
CON(18)= coefficient of sine annual term
CON(19)= time variation of CON(18)
CON(20)= coefficient of sine semi-annual term
CON(21)= time variation of CON(20)
if JCT(29) > 0:
CON(1-2)= free core nutation, cos & sin components.  See JCT(28).
period (approx. 460 days) is given by CON1(2)
CON(3-10) coefficients of adjustable terms in nutation series
CON(3)  = coefficient of cos(2F +2 omega) term in deps  (13.7 days
CON(4)  = coefficient of sin(2F +2 omega) term in dpsi
CON(5)  = coefficient of cos(2L -2D +2 omega) term in deps (183 da
CON(6)  = coefficient of sin(2L -2D +2 omega) term in dpsi
CON(7)  = coefficient of cos(Lp) term in deps  (365 days)
CON(8)  = coefficient of sin(Lp) term in dpsi
CON(9)  = coefficient of cos(omega) term in deps  (18.6 yrs)
CON(10) = coefficient of sin(omega) term in dpsi
CON(11-18) coefficients of analogous out-of-phase terms--
sin terms in deps, and cos  terms in dpsi
CON(19) = k/c for fortnightly terms in ut1
CON(20) = k/c for monthly terms in ut1
For any value of JCT(29):
CON(22)= precession constant (if zero nominal value used)
CON(23)= obliquity constant (if zero, nominal value used)
CON(24)= nutation constant (if zero, nominal value used)
CON1(1) = epoch for polynomial variation, CON(11-13).
and for seasonal parameters, CON(14-21). for before 1956
hardwired in ctutf as 1956 jan 17.0
should be jan 17.0 of whatever year is chosen, so that
the pre- and post-1956 seasonal terms are in phase
= 0 corrections not applied(see a1ut1f & wobblf)
CON1(2) = period(days) of free mode nutation, CON(1-2).

moment of inertia ratios
CON(2) = alpha
CON(3) = beta
CON(4) = gamma

CON(5) = meqinc, mean inclination of equator to orbit
CON(6) = k2, lunar love number
CON(7) = k2*t or k2/q, lunar dissipation parameter.  see k(83).
CON(8-10) = coefficients of analytic dissipation terms
CON(11-16) = initial conditions of lunar core rotation
CON(17)= core-mantle rotation coupling constant
(solar-mass-AU-squared-per-day)
CON(18)= core moment of inertia (solar-mass-AU-squared)
CON(19)= core flattening (ignored in computing mantle moments)

For limited asteroids (SMALL=t)
initial conditions are specified as for other bodies
NPLNT = asteroid number
CON(1) = radius (km)
CON1(1)= epoch of initial conditions (must be midnight)
denptr = index into prmter array for density parameter.
logic elsewhere limits denptr to 34-38.

For pulsars (PULSAR=t)
NAME - 1st four characters are pulsar spot name
JD0   = reference epoch for pulse phase model
NTYPE = code for kind of phase model (not used)
CON1(2)= approximate pulse period (sec)
CON(1)= parallax (arcsec)
CON(2)= ra proper motion (arcsec/yr)
CON(3)= dec proper motion (arcsec/yr)
CON(4)= dispersion (sec-hz**2)
CON(5)= pulse phase at reference epoch jd0 (cycles)
CON(6)= pulse period - CON1(2) (sec)
CON(7)= pulse period drift rate (sec/sec)
CON(8)= pulse period acceleration (sec/sec**2)
CON(9)= time derivative of dispersion (sec-Hz**2/d)
CON(10)= amplitude of the planet/companion's signature (sec)
CON(11)= planet orbit eccentricity
CON(12)= planet argument of periapse (deg)
CON(13)= planet mean anomaly at JD0 (deg)
CON(14)= planet orbital period (days)
CON(15)= planet orbital inclination (deg)
CON(16)= planet mass (Msun)
CON(17)= time derivative of planet's signature (sec/d)
CON(18)= time derivative of planet's periapse (deg/d)

L(1) to L(30) control whether initial conditions COND(1-6) and
parameters CON(1-24) belonging to given body are adjusted in the
iterative least squares process.
(a) for i=1,6, l(i)=0 or 1 according to whether the
corresponding initial condition  COND(i) is not adjusted
or is adjusted in the least-squares analysis
(b) for i=7,30,l(i)=k means that parameter CON(k) is
adjusted,k=1,24.if l(i)=0,no such parameter is adjusted.
L(7-30) must satisfy same logical conditions as
LPRM(1-100), namely
(1) if j.gt.i and L(i)=0, then L(j)=0
(2) if j.gt.i and L(i),L(j) not 0, then L(j).gt.L(i)

EPS(1) = start of integration output tape is Julian day number JD1
and fraction of day EPS(1) (EPS(1) assumed zero in
non satellite-probe case)
EPS(2) = end of integration output tape is Julian day number JD2
and fraction of day EPS(2) (EPS(2) assumed zero in
non satellite-probe case)
EPS(3) = accuracy constant to control numerical integration step size
for Nordsieck method (K(88)=0,1 or K(88)=2,3 in starting
proceedure)
EPS(6) = if JD0 negative for checkpoint restart, numerical integration
restarted at Julian day number -JD0 and fraction of day
EPS(6) or as close after this time as is practicable.

K(1) to K(100) control integration of equations of motion and
equations for partials of motion for body with planet number NPLNT.
K(1-30) were formerly used to select partial derivatives for
integration.  This function is now filled by KI(1-NUMKI).  If NUMKI
is not specified, then the program logic sets NUMKI to the index of
the highest-numbered non-zero KI.  Old-format input streams will be
recognized by the lack of KI's coupled with values .le.1 for all of
K(1-7).  Old-format input streams will be translated automatically
(and internally) to the new format.

KI(1)
KI(1) =-1 equations of motion are integrated and partial derivatives
w.r.t. initial conditions are determined from elliptic
or mean orbit formulas (depending on value of  K(100))
KI(1) = 0 equations of motion are integrated
KI(1) = 1 equations of motion and equations for partial derivatives
w.r.t. initial conditions are integrated

KI(2-7)
KI(1+j)=-1   partial derivatives with respect to initial orbital
element j (or initial position or velocity coordinate
j if NPLNT.gt.30) are not determined, j=1,...,6
KI(1+j)= 0,1 partial derivatives with respect to initial orbital
element j (or initial position or velocity coordinate
j if NPLNT.gt.30) are determined if KI(1) is not zero,
j=1,...,6

KI(8-NUMKI)
KI(j), j=8,...,NUMKI, determine parameters for which equations for
partial derivatives are integrated. we must have
(1) if j greater than i, then KI(j) greater than KI(i)
unless (a) KI(j)=0 or unless (b) KI(j) is negative,
in which case KI(i) must also be negative with KI(i)
greater than KI(j)  except that KI(k+1),KI(k+2)
are not considered in this logic if KI(k)=100*n+31
nor are KI(k+1),KI(k+1),KI(k+3),KI(k+4) considered
in this logic if KI(k)=100*n+41 or KI(k)=100*n+51
(k=8,...,NUMKI, n=1,2,...)
(2) if j greater than i and KI(i)=0, then KI(j)=0
with these conventions, we then have for j=8,...,NUMKI
KI(j) = 0 no equations integrated
KI(j) =-m (m positive) equations for partial derivatives with respect
to parameters  CON(m) are integrated (m=1,24)
KI(j) = m equations for partial derivatives w.r.t. PRMTER(m) are
integrated, m=1,...,100 (really m=1,...,50 because
prmter(1-50) are supposed to affect motion and
prmter(51-100) are supposed not to affect motion but
only observations)
KI(j) = 100*n+i    equations for partial derivatives with respect to
initial condition i (i=1,6) or parameter i-6 (i=
7,30) for perturbing target planet or central
planet n (n=1,10) are integrated
KI(j) = 100*n+nn
nn=31      equations for partial derivatives w.r.t. zonal
harmonic coefficients between orders KI(j+1) and
KI(j+2) inclusive for target or central planet n
are integrated
nn=40-44   equations for partial derivatives w.r.t. tesseral
harmonic cosine coefficients between orders
KI(j+1),KI(j+2) and KI(j+3),KI(j+4) inclusive for
target or central planet n are integrated.
Warning: only 41 working in COMPAR (PARTL1/HPARTL)
nn=45-49   equations for partial derivatives w.r.t. resonant
tesseral harmonic cosine coefficients of order
KI(j+1) and of degree KI(j+2) thru KI(j+3).
Not working in COMPAR.
nn=50-54   equations for partial derivatives w.r.t. tesseral
harmonic  sine  coefficients between orders
KI(j+1),KI(j+2) and KI(j+3),KI(j+4) inclusive for
target or central planet n are integrated.
Warning: only 51 working in COMPAR (PARTL1/HPARTL)
nn=55-59   equations for partial derivatives w.r.t. resonant
tesseral harmonic sine coefficients of order
KI(j+1) and of degree KI(j+2) thru KI(j+3).
Not working in COMPAR.

K(30)=-1  effect of any input limited asteroids not included
K(30)= 0  effect of all input limited asteroids (if any) included in
equations of motion but not partial derivatives
K(30)= 1  effect of all input limited asteroids (if any) included in
equations of motion and partial derivatives (not implem.)
K(31-60)
K( 30+j)=-1 effect of mass(j) not included in equations of motion or
in equations for partial derivatives, j=1,...,30
K( 30+j)= 0 effect of mass(j) included in equations of motion but not
in equations for partial derivatives, j=1,...,30
K( 30+j)= 1 effect of mass(j) included in equations of motion and
in equations for partial derivatives, j=1,...,30
if NPLNT between 1 and 30, K(30+NPLNT) is ignored
if ncentr is between 1 and 30 (ncentr assumed 3 for NPLNT=10) we have
K(30+ncentr)=-1 effect of sun not included in equations of motion nor
in equations for partials of satellite
K(30+ncentr)= 0 effect of sun included in equations of motion but
not in equations for partials of satellite
K(30+ncentr)= 1 effect of sun included in equations of motion and
in equations for partials of satellite

K(61-80)
K( 30+j)=-1 effect of PRMTER(j) not included in equations of motion
or in equations for partial derivatives, j=31,...,50
K( 30+j)= 0 effect of PRMTER(j) included in equations of motion but
not in equations for partial derivatives, j=31,...,50
K( 30+j)= 1 effect of PRMTER(j) included in equations of motion and
in equations for partial derivatives, j=31,...,50
Note: whether  K(61) controls effect of PRMTER(31) or controls effect
of  CON(24) depends on value of  K(97).
Also, for sun-centered bodies, K(61)=2 selects a simplified form for
the general relativity correction to the force which only includes
the sun and omits terms proportional to the velocity of the sun
relative to the solar system center of mass.  This simplified
formulation is the only one available for non-sun-centered bodies.
K(61)=2 includes relativity in partials.
Note: K(80) controls whether asteroid ring is included in equations
of motion in solprb (PRMTER(50) is mass of asteroid ring).
Also coded in PLNORB.

For a planet (NPLNT.gt.0 and NPLNT.le.30)
K(82)=-1,0,1   according to whether the interaction of the
planet's second zonal harmonic with the sun
is not included, is included in the motion but
not the partials, or is included in motion and
partials (this last option not yet coded)

For probe (NPLNT.gt.30) with logical central body the sun (SOLPRB)
(also for planet, asteroid, or comet)
K(81) =-1   no low thrust forces
K(81) = 0   low thrust forces included in motion, not in partials
K(81) = 1   low thrust forces included in motion and partials

For a probe with logical central body a planet (ERTROB or PLNORB)
K(81)- K(84) control for various low thrust forces
K(81)= 10*kreflt+ kdirct
kdirct.le.0  no direct radiation pressure
kdirct= 1 direct radiation pressure in motion

reflected radiation not yet coded
kreflt.le.0 no radiation from central body
kreflt=1 reflected radiation from planet in motion
kreflt=2 emitted infared radiation from planet in motion
kreflt=3 reflected and emitted radiations included in motion
K(82).le.0 no atmospheric drag
K(82)=+1   effect of drag on motion and partials
K(82).gt.1 reserved for table lookup of density in atmden
(not coded)
K(83).le.0 no gas leak forces included
K(83)=1  effect of gas leaks on motion included

K(86).lt.0 default number of central body harmonics in sbset included
in integration of equations for partials
default=j2 only (+c22,s22 if moon=central body)
K(86).ge.0, K(86)=ntess*100+nzone   ntess tesserals and nzone zonals
of central body harmonics included in partials integration

K(85) is exactly the same for target bodies as K(86) is for central
bodies. no distinction made in K(85) between target bodies,
the decision as to target body harmonic order included
in partials integration is the same for all of them.
default=no harmonics

For Moon (NPLNT=10)
K(81)=n, where n is the degree of the highest zonal harmonic
included in the Earth's potential.
K(82)=m, where m is the maximum degree of the lunar gravity field
used. both zonals and tesserals are carried through the
mth order.
K(84)=-1,0,1 according to whether tidal friction does not affect
motion or partials, affects motion but not partials, or
affects motion and partials.
K(85)=-1,0,1 according to whether a violation of the principle of
equivalence does not affect motion or partials, affects
motion but not partials, or affects motion and partials
(this switch controls the contribution from CON(19)
and CON(20), plus the contribution from the metric
parameters beta and gamma, PRMTER(43-44) )

For Moon rotation (NPLNT=-10)
K(81)=-1,0,1 according to whether the sun does not affect the
rotation, affects the rotation e.o.m. only, or affects
the rotation eom and the variational equations.
K(82)=-1,0 according to whether or not the Earth-figure torque
is included in the equations of motion.
K(83)=-1,0,2 implying, respectively, no lunar elasticity &
dissipation (rigid-body), elas. & constant-time-lag diss.,
elas. & constant-q dissipation model.
1: same as 0, but calculate lunar mean motion from masses
K(84)=-1,0,1 same as K(81), but for planets Venus through Saturn
K(86)  controls lunar harmonics for partials, same as in probes
(but use K(86) of Moon if higher)

note-- if orbit and rotation are simultaneously integrated
(both JD0s non-zero) then the integration control info.
comes from NPLNT=10 stream. this includes JD0, JD1, JD2,
EPS(3), K(87-92), and also K(31-40).

K(87)   = positive, step size for Adams-Moulton or second sum
numerical integration is  K(87) days
K(87)  .le.zero, interval size for Adams-Moulton or second sum
numerical integration is 2** K(87) days

K(88) = Indicator of the type of integration, if any, to be done.
Also saved in NAMTIM common as INTTYP for each body.
The default value is 3.
K(88) =-9 no integration: orbit is approximated by a quasi-elliptic
model specified in COND(1-6) and TCON(1-10). JD0 had
better be 0!
K(88) =-8 no integration: orbit is approximated by a purely elliptic
model with elements COND(1-6) and epoch CON1(1).  JD0 had
better be 0!
K(88) = 0 Nordsieck variable step size method used for numerical
integration with variable tabular interval output
K(88) = 1 Nordsieck variable step size method used for numerical
integration with constant tabular interval output
K(88) = 2 Adams-Moulton method used for numerical integration
K(88) = 3 Royal Road (second sum) method used for numerical integ.
Both Adams-Moulton and Royal Road are constant step integration
methods with constant tabular interval output, and they both call
Nordsieck method in their starting proceedures (Nordsieck is self
starting).

K(97) = 0 relativity motion factor for force controlled by  K(61)
is PRMTER(31)
K(97) = 1 relativity motion factor for force controlled by  K(61)
is  CON(24) inserted into PRMTER(31) in sub.PLANET for fn

K(98).le.0 no printout of partial derivatives
K(98).gt.0 printout of partial derivatives
K(98) = also controls how often data is printed out, namely
every iabs(K(98)) tabular points if K(98).ne.0
and every tabular point of integration if K(98)=0

K(100)=-1 numerical integration of usual equations of motion
K(100)= 0 numerical integration of equations for difference between
true orbit and initial osculating elliptic orbit
K(100)= 1 numerical integration of equations for difference between
true orbit and mean orbit
(For planets, K(100) assumed zero for fn integration)
(For Moon rotation, K(100)=-1,0,1 signifies Euler angles,
osculating Cassini angles, or mean Cassini angles)

KK(2)  = 0 never any output on KOUT.ne.0 from SBEXP extended print
called from SBOUT during SBFN probe or planet integration
.lt.0 SBEXP called every -KK(2)th apsis
.gt.0 SBEXP called every KK(2)th integration step

KK(3) is a packed bits number
1 cartesian coordinates are output
2 all partials are output
4 apsidal time is printed if near
8 elliptic elements referred to the mean equinox and equator
(of the Earth) of the reference epoch are output
not programmed as yet
16 print elliptic elements  central body of date
32 print elliptic elements  plane of the sky

KK(4) is a packed bits number
1 print  Earth - s/c vector
2 print  Earth - central body vector
4 print   sun  - s/c vecyor
8 print   sun  - central body vector
16 print target - s/c vector

KK(5) is a packed bits number
1 print  latitude & longitude of s/c w.r.t. central body
2 print projected motion of s/c on 'stationary' central body
print projected motion of s/c on   rotating   central body

KK(6) is a packed bits number
1 print vector of s/c w.r.t. Arecibo at highest signal
2 print vector of s/c w.r.t. Earth at close approach
4 print vector of s/c w.r.t. 1st target at close approach
8 ditto 2nd target
...
8192 ditto 12th target

KK(7) is a packed bits number for Kout extra print
1 print accelerations at every tabular point (coords+partials)
2 print relativistic accelerations and (1st only) partials
4 print coordinates and gravity gradients w.r.t. them for
three bodies: integrated, 1st target, center
when integrating embary, 3rd body is moon
8 print partials of coordinates of same three bodies w.r.t.
1st parameter
16 print newtonian indirect partials contributions due to
positions and velocities of same three bodies, and then
relativistic likewise (note newtonian velocity indirect
contribution is always zero)
32 change interval of printout to 100 days
64 reserved for nbody extra print
128 suppress printout of moon libration angles

KK(11) = moon core integration output dataset (for moon rotation)
KK(12) = moon core integration forward/backward buffer dataset
KK(13) = moon libration integration forward/backward buffer dataset
(supersedes ibuf if nonzero)

KK(50) =   number of integration steps per printing of
the dadx# matrices.  if .lt. zero, print for
iteration 1 only.  control is bl pntflg, set
in sbfn at "once per iteration for partials"

extra control for Nordsieck variable stepsize integration
see TCON(1-4) below
KK(51) =  0 stepsize control not nused (default)
-1 make transition to stepsize control
1 lock onto stepsize control

KK(78) = 0 no data set with radiation pressure model on it
KK(78).gt.0 data set number for cards from which mvmsrf interpolates
solar radiation force
KK(78).gt.0 in ERTORB, use NDS radiation pressure model

KK(81) = 0 no printout of distance to central body
KK(81) = 1 printout of distance to central body in kilometers
if this distance is less than CON1(12) au (central
body not sun)

KK(84) = 0  read no center/target body tape during integration
KK(84) = 1  read tape for center/target body position, velocity only
KK(84) = 2  read tape for partials as well
(option 1 above not implemented: treated as option 2)
default is 2 for the Moon, 0 otherwise

KK(85)     hsitb=distance from target body at which the
effect of the harmonics(not j2)on the e. of m. and partials
is disabled.  it is nominally set so that the
maximum possible harmonic-acceleration is less than
the solar-acceleration by hsftr=1.0e15
hsitb=nominal+abs(KK(85))/100.    (a.u.)

KK(88) =   intm in EVAL
0 correct 2 times - old form - default
-1 Royal Road mode - 1 correction
1 like -1 plus pseudo correction with dadx matrix
not tested and evidently not working: KK(88)=1

KK(91) = 0 no special LES-8/9 forces in equations of motion
KK(91) =-1 LES-8/9 radiation pressure, no thrusts, standard attitude
KK(91).gt.0 LES-8/9 radiation pressure with thrust and attitude
history read from data set KK(91), which is written in
subroutine lesin from reading cards following this &nmlst2

KK(93) = 0 no change of numerical integration method in checkpoint
restart
KK(93) = 1 such a change in subroutine sbset with input values of
EPS(3),K(87-92) overriding values on second record of
restart tape

KK(94) = 0 no LES-8/9 station keeping thrusts
KK(94) = 1 LES-8/9 station keeping emulation is turned on
in numerical integration of equations of motion
KK(94) = 2 emulation turned on plus write out thrust events

KK(95) =   starting interval is 2**KK(95) (for KK(95) negative)
for variable step size integration starting proceedure
at thrust initiation and termination  for LES-8/9
station keeping emulation. K(91) is still starting
interval for starting proceedure at initial epoch.

KK(97) = 0 no partials integrated in reintegration of equations of
motion after orbit fit convergence
KK(97) = 1 partials are integrated in reintegration of equations of
motion after orbit fit convergence

KK(100)= 0 partial derivatives are integrated each iteration if
K(1-30) so indicate
KK(100).gt.0 partial derivatives are not integrated if iterat.gt.
KK(100) (iterat is least squares interation conter)
(integration will do what K(1-30) say to do when iterat=1
no matter what the value of KK(100))

the following four TCONs are used for Nordsieck integration
in conjunction with stepsize control option - see KK(51)
TCON(1)= scaling variable, approx eq semi-major axis of probe in au
TCON(2)= period of probe in days
TCON(3)= approximate number of integration steps per period
TCON(4)= exponential index of (r/a)

for non-integrated bodies TCON(1-10) are parameters of
a precessing quasi-elliptic orbit (all deg and day)
TCON(1)= i0, inclination of orbit to laplacian plane
TCON(2)= ldot, mean motion in longitude
TCON(3)= pdot, rate of longitude of periapse
TCON(4)= kdot, rate of argument of node on laplacian plane
TCON(5)= ndot, rate of node of laplacian plane on earth equator
TCON(6)= jdot, rate of inclination of laplacian plane
TCON(7)= ldotdot, 1/2 rate of change of mean motion
TCON(8)= theta, amplitude of term sin(k-phi)
TCON(9)= phi, phase of term sin(k-phi)
TCON(10)=k0, initial argument of node on laplacian plane

NZONE =order of zonal harmonic input for gravitational potential
of given body
Ji    =zonal harmonic coefficient for given body, i=2,3,4,...
LJi   =0,1 according to whether corresponding zonal harmonic is not
or is adjusted in least squares analysis

NTESS =order of tesseral harmonic input for gravitiational potential
of given body
Ci(k) =tesseral cosine harmoic coefficient for given body,
i=2,3,4,... and k=1,...,i
Si(k) =tesseral sine harmononic coefficient for given body,
i=2,3,4,... and k=1,...,i
LCi(k)=0,1 according to whether corresponding tesseral cosine
harmonic is not or is adjusted in least squares analysis
LSi(k)=0,1 according to whether corresponding tesseral sine
harmonic is not or is adjusted in least squares analysis

There are 4 planet shape models programmed in the radar link of PEP.
CON(2,3 or 4) for NPLNT.gt.0 turns on shape for radar observations of
a planet.  The model is given as harmonic input for body -NPLNT in
one of the following manners:
1)shperical harmonic expansion, if nshp=0 (default)
2)fourier series model, if nshp=1
3)altitude grid (local model), if nshp=2
4)external model if nshp=3

In the spherical harmonic expansion over the whole planet surface,
CON(2-4) for NPLNT are J1,C11,S11 with Jn,Cmj,Smj (n.le.NZONE, m.le.
NTESS, j.le.m) for -NPLNT being the remaining harmonic coefficients.
normalized coefficients are used for shape, whereas unnormalized are
used for the gravitational potential. for the spherical harmonic model
planet radius = CON(1) + spherical harmonic series
(units are seconds in solution for harmonic coef,  km for CON(1)  )

If NSHP=1 for -NPLNT the Fourier series model is used with
TLAT(1) and TLAT(2) equal to the lower and upper latitude
values for the planet NPLNT being analyzed.
Let phi=east long on planet,
theta=(latitude-tlat(1))/(tlat(2)-tlat(1))*twopi
planet radius = function of CON(1) mean radius and CON(2) degree of
flattening + fourier series coefficients for -NPLNT
the shape model, as it is currently, is a 2-dimensional fourier series
with 20 longitudinal terms (smallest period = 360/20=18 deg) and
1 latitudinal term.  it is:  h(theta,phi)= summation for m=1,20 of
(am*cos(m*phi)*cos(theta)  +  bm*cos(m*phi)*sin(theta)
+  cm*sin(m*phi)*cos(theta)  +  dm*sin(m*phi)*sin(theta)
+  a'm*cos(m*phi)  +  c'm*sin(m*phi) )
+  a''*cos(theta)  +  b''*sin(theta)
notice that the a' and c' coefficients are just longitude dependent,
and the a'' and b'' are the only two latitude-only terms.
all 122 fourier series and radius & flattening partial derivatives
should be calculated when computing the observed minus theoretical
values of observations.  suggested parameters to be solved for
in addition to radius and perhaps flattening are:
Mercury - only low order longitudinal terms, maybe a'1 thru a'5 and
c'1 thru c'5.
Venus   - how about a'1 thru 20 and c'1 thru 20?
Mars    - probably all 122 fourier shape parameters

Summary: to use the fourier shape section of coding you set NSHP=1,
TLAT(1&2) equal to lower and upper latitude values for the rectangle
(almost) on the planet being analyzed in &NMLST2 for -NPLNT.
Default values presently set are:
Mercury   TLAT(1)=-15., TLAT(2)=15.,
Venus     TLAT(1)=-15., TLAT(2)=15.,
Mars      TLAT(1)=-30., TLAT(2)=30.,
In the corresponding &NMLST2 with positive NPLNT, CON(2)=1.0E-20_10
sets npshp=1 so HARSHP is called, and L(7)=1, as usual, adjusts the
radius. In addition L(8)=2 adjust flattening CON(2), which could
have a real value instead of the dummy 1.E-20_10 used as a switch.
Units of radius are km, flattening is dimensionless, and
units of Fourier seriers coefficients for radius are light-seconds.

The Fourier series coefficients and L-vectors may be input
directly as the following &NMLST2 variables:
coefficients                                       l-vectors
fa(1)  to fa(20)         for a coefficients            lfa
fb(1)  to fb(20)         for b coefficients            lfb
fc(1)  to fc(20)         for c coefficients            lfc
fd(1)  to fd(20)         for d coefficients            lfd
fap(1) to fap(20)        for a' coefficients           lfap
fcp(1) to fcp(20)        for c' coefficients           lfcp
fapp                     for a'' coefficient           lfapp
fbpp                     for b'' coefficient           lfbpp
to adjust coefficients, set appropriate l-vector to 1
***** example *****
mercury fourier series shape
solve for a'1 thru a'5 and c'1 thru c'5
required input is:
&nmlst2
NPLNT= -1,name= 'shpmrcry',
nshp= 1,
lfap= 1,1,1,1,1,     (or lfap= 5*1,)
lfcp= 1,1,1,1,1,     (or lfcp= 5*1,)
&end

if nshp=2 for -NPLNT the shape model is a grid of
altitudes (in kilometers) relative to an adjusted spheroid
for a appropriate region of the planet NPLNT. the region
of applicability for the model is input as the real*4
nmlst2 variables tlat(1),tlat(2),tlon(1),tlon(2) representing
minimum and maximum values for latitude and longitude, respectively
(in degrees). tlatin and tlonin give the increment in
degrees for lat and long (ie the grid spacing). default
values for the lat. range are as given above for the fourier
series model. default for long is 0 to 360 degrees (east)
grid spacing defaults are 10 degrees. (latdim,londim) the dimensions
of the grid and the total number of grid points (ngdpts) are computed
from tlat,tlon,tlatin,tlonin.
to input a priori estimates for grid points (in km),
determine the corresponding index into the (1-dimensional)'grid'
array where the 2-dimensional grid points are stored by
row (ie., in latitude strips). grid(1) corresponds to the
grid point of min. long and max lat (ie., the upper left
point). lgrid is the l-vector for the grid points with
lgrid(i)=0, do not adjust grid(i)
=1,        adjust grid(i)
as in the fourier shape model, the flattening coefficient CON(2)
for the body |NPLNT| should be set to some small non-zero value,
such as 1E-20_10 to enable the shape computations in the COMPAR link.
note that the flattening may also be treated as an adjustable
parameter as well as a switch.