pnitl - j.f.chandler - 1984 august
based on subroutine sate from duxbury and callahan 1982
compute position of satellite using precessing elliptic
elements in standard frame

setup is performed by calling pnitl:
arguments
icnd1 type of initial conditions (should be 1, for longitudes)
cond1 initial conditions (ignoring precession and inclination to the
laplacian plane)
cond2 inclination to laplacian plane, angle rates, etc.
setp  array to be filled with setup quantities

evaluation is performed by calling plipt:
dt    time from initial epoch in days
x     returned array of position or velocity
ic    indicator: 0 -> position, -1 -> velocity

use of array 'setp' --- storage for saved quantities from
setting up.  it is effectively equivalenced to the following:
b(3,2),tsv,a,e,anom0,motn,perrat,krat,nrat,jrat,mot2,fdga,fdgp,
1      7   8 9  10    11    12    13   14   15   16   17   18

note: 'j' and 'n' are the node and inclination of the laplacian plane.
'i' is the inclination of the orbit to the laplacian plane.
'a' and 'e' are the usual semimajor axis and eccentricity.
'k' is the argument of node, not its longitude, while 'per' is
the longitude of periapse and 'long' is the mean longitude.
'fdga' and 'fdgp' represent a solar perturbation term for
deimos (wilken's theory)