compute general relativistic correction to force on moon.
based on eq. (6.34) of will, theory and experiment in
gravitational physics, cambridge univ. press, 1981.
(actually this equation is in isotropic rather than harmonic
coordinates, but these are the same to 1st post-newtonian order).

ncall=-2 setup once per iteration of a given step for partials
(not used - already done with motion)
ncall=-1 setup once per iteration of a given step for motion
ncall= 0 setup once per step
ncall= k evaluate acceleration for equation k.gt.0

rvec(4-6) not filled in, except for earth-sun-moon combinations
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

note on subscripts in this routine:
a,b,c (integers) match formula (6.34) or (6.78) in will
a is equivalenced to nplnt
1-9 refer to planets (3 = earth)
10  is sun, not moon
11  is integrated body (moon)

However, after computing sumrp (partials w.r.t. masses of these
bodies) in terms of the above subscript convention, the results
for 3 and 11 are converted into partials w.r.t. MASS(3) and
MASS(10) and stored in the 3rd and 10th slots.  Similarly,
the Dadxrp and Dadvrp arrays (partials w.r.t. coordinates) are
converted for 3 and 11 (Earth and Moon vs Sun) into 3 and 10
(Embary vs Sun and Moon vs Earth).