R. Goldstein   July 1977
compar link implementation of planetary rotations
Computes the spot coordinates referred to equatorial inertial
frame from the body fixed coordinates.
Entry SPOTPR computes the partials of these inertial
spot coordinates wrt the free parameter of the rotation model
(CON(6)-CON(17)).
The array ROT is set up  to
allow SPOTCD to calculate partials of the inertial spot
coords wrt, spot lat, long, and radius.

JDT       Julian day number of time when coord. desired
FRT       fraction of Julian day of same
YSPCD     input body fixed spot coords.
XSPCD     output equatorial inertial coords. of spot
NI        input index (1 or 2) of XSPCD and YSPCD
ROT       output rotation matrix
ICT66     input I*4 version of ICT(66) documented in PRMRED,
mainly controls debug printout
IF        0: interpolation used to calculate INC, PSI,A,B
1: full calculation of above for every observation

see if the time has changed so that calculation
for inc, psi, a and b must be redone
tlim is set so that they are calculated once per
observation. when they are calculated, it is by
interpolation between values calculated at the
beginning and end of a "pass"

set up rot to send to the rest of pep. rot is
needed to calculate partial wrt spot coords. the rot that
is needed is the transpose of the rot in the above
cited references due to fortran row and col. convention.

set up fudge factors
the minus sign in these factors are due to the following:
obs = k*x(s/c-earth)
x(s/c-earth) =  terms  -x(s/c-mars)
the construction of the full partials are done in cpartc

at this point, the vector lrotdp has been set up which
indicates which partials to calculate in the order that
they must be calculated. it does not necessarily correspond to
lpl since partials wrt some parameters depend on calculations
performed when calculating partials wrt other parameters. psav
contains saved results that are used in later partials.
order of partials: phi0, p, i0, psi0, mu, delta0, alpha0

psav:results of operating on yspcd with the following operators
1: prot3(phi)
2: rot1(i)*prot3(phi)
3: rot3(psi)*rot1(i)*prot3(phi)
4: prot1(i)*rot3(phi)
5: rot3(psi)*prot1(i)*rot3(phi)
6: prot3(psi)*rot1(i)*rot3(phi)

apply rfix and calculate final partial derivatives. fill
derem...note...velocity part of derem not calculated as of
july, 1977. this implies that only phase delay method of
doppler calculation can be used  (ict(21)=-1).