r.king and r.cappallo   july 1977   real*10 function morfn
evaluation of right side of differential equations for moon orbit
and rotation.  this routine based on monfn, written by m.ash and
m.slade oct 68, and mrtfn, written by king and cappallo may 76.

ecor  = earth relative to sun (re is distance) = rvec(.,3,10)
mcor  = moon relative to sun (rm is distance) = rvec(.,11,10)
mecor = moon relative to earth (rem is distance) = rvec(.,11,3)
mcor etc. are equatorial coordinates
smcor etc. are selenodetic coordinates

solar radiation pressure model
include effect of absorption and reradiation
assume negligible thermal inertia
solar luminosity is given in equivalent solar mass per year
possible refinements: thermal inertia for earth
some day, the solar wind may be nonnegligible

the following calculations done six times, once for the torque due
to the earth (is=1), once for the torque due to the sun (is=2),
and once for each of four planets
body-specific quantities used:
spcor = selenodetic coordinates of body w.r.t. moon
mpcor = equatorial coordinates of body w.r.t. moon
h     = force on body due to lunar figure
gamtrq= g times mass of body
masfct= scaling for non-earth bodies, due to embary offset

first get partials of angular accelerations
ptrmq already has the partials of the harmonic terms of torque
and is equivalenced with the first half of the array.
Clear the second half of the array and add the partials of the
kinematic terms. For historical reasons,
the angle partials are denoted 2,4,5, and the angular velocity
partials 1,3,6

get mantle-core torque due to friction and flattening, but note:
flattening torque cancels out against kinematic cross term if integrating
core Euler angles
frame is core if integrating core Euler angles, or mantle if core spin