{ the forces due to planets, the Moon, the Sun
   the right hand side of differential equations }

 UNIT unequras ;

 INTERFACE

 Uses
   unsagcot ; { const types }

 Procedure PlanetsForces ( TEph  : Extended ;
                           CurP  : tvect3 ;
                           Var  ForceC  : tvect3 ) ;

 Procedure ClcRightHandCoor ( ut   : Extended ;
                              rc   : TVect3 ;
                              Var
                                force  : TVect3 ) ;

 IMPLEMENTATION

 uses
   unsagvar , { global variable PlanetMain }
   unephjpl , { to use planet numbers }
   ungephem , { procedure forgiapos }
   unfortim , { procedures TransUTCtoTDT TransTDTtoTDB }
   ungiarom , { procedure ToClcMatrOfRotation }
   ungiamod ; { global variables for gravitational model of the planet }

 { CurP the planetocentric
   position of the satellite
   in the fixed equatorial system J2000 }

 Procedure PlanetsForces ( TEph  : Extended ; { moment in julian day }
                           CurP  : tvect3 ;   { equatorial position AU }
                           Var  ForceC  : tvect3 ) ; { only the Sun }
 Var
   i  : byte ; { very simple count }
   m  : byte ; { number of main planet }
   fm    : extended ; { instead of FactorMuZ }
   dr,sr : extended ; { temporary for distance }
   posp  : tvect3 ;   { main planet heliocentric position in AU }
   posa  : tvect3 ;   { unsagcot for types const for the Sun }
 Begin
   { posp position of the Main Planet in AU refer to the Sun }
   m:=PlanetMain; { global PlanetMain from unsagvar }
   fm:=FactorMuZ[m]; { for the Sun as main forces in AU^3/day^2 }
   forgiapos(m,teph,posp); { ungephem }
   { only for pur main giant planet  }
   for i:=1 to 3 do     { refer to the main planet }
     posa[i]:=-posp[i]; { for the Sun position }
   dr:=1/Sqrt(+Sqr(posa[1]-CurP[1])   { for current planet }
              +Sqr(posa[2]-CurP[2])   { and the satellite }
              +Sqr(posa[3]-CurP[3])); { 1/delta in km }
   dr:=dr*dr*dr; { 1/delta**3 }
   sr:=1/Sqrt(Sqr(posa[1])+Sqr(posa[2])+Sqr(posa[3])); { 1/r }
   sr:=sr*sr*sr; { 1/r**3 }
   For i:=1 To 3 Do  { simple formulae forces }
     ForceC[i]:=+ForceC[i]
                +fm*(+(posa[i]-CurP[i])*dr-posa[i]*sr);
 End;

 { attraction from the main planet
   planetocentric unsagvar PlanetMain global
   from the main planet with zonal harmonics }

 procedure ForPlanetMain ( tb  : extended ;  { time dynamic barycentric }
                           rc  : tvect3 ;    { planetocentric J2000 }
                           var  force  : tvect3 ) ; { our forces }
 Var
   i,j,NHarm,NPlanet        : Byte ;
   TinDay,
   RRat,RRatN,
   PertJ,PertZ,fmr,
   p,q,r,rr,rrr,zr,ra,rb    : Extended ;
   PolLeg0,
   PolLeg1,PolLeg2,
   DPolLeg1,DPolLeg2        : Extended ; { for the Legandre polinoms }
   RotMatr                  : TypeDim33; { to rotate to planet equator }
   PosFix,PosSat,FSatEq     : TypeDim3 ; { type unsagcot }
 Begin

   TinDay:=tb; { time dynamic barycentric in julian day }

   { matrix from fixed equator J2000 to the equator of the planet }
   ToClcMatrOfRotation(PlanetMain,TinDay,RotMatr); { ungiarom procedure }
                                       { J2000 equator system but }
   For i:=1 To 3 Do  PosFix[i]:=rc[i]; { satellite refer to planet }

   { from fixed equator J2000 to the equator of the planet }
   For i:=1 To 3 Do
     Begin
       PosSat[i]:=0.0;  { to rotate }
       For j:=1 To 3 Do { for each coordinates }
         PosSat[i]:=PosSat[i]+RotMatr[i,j]*PosFix[j];
     End;

   { range of the satellite from the centre of the planet main }
   rr:=1/(Sqr(PosSat[1])+Sqr(PosSat[2])+Sqr(PosSat[3])); { 1/r^2 }
   r:=Sqrt(rr); rrr:=r*rr; { 1/r , 1/r**2 , 1/r**3 }

   { to take into account the influence of zonal harmonics }
   zr:=PosSat[3]*r; { z/r for Legendre polinoms }
   PolLeg0:=1;  PolLeg1:=zr; { initial values }
   DPolLeg1:=1; { for recurrent formula }
   RRat:=PlanetRadius*r;  { r0/r global ungiamod }
   RRatN:=RRat;
   PertJ:=0; PertZ:=0;
   For NHarm:=2 To MaxZonHarm Do
     Begin
       RRatN:=RRat*RRatN; { (r0/r)**NHarm }
       ra:=ZonHarmValue[NHarm]*RRatN;  { global ungiamod }
       p:=NHarm; q:=1/p;
       PolLeg2:=(2*p-1)*q*zr*PolLeg1-(p-1)*q*PolLeg0;
       DPolLeg2:=p*PolLeg1+zr*DPolLeg1;
       rb:=(p+1)*PolLeg2+zr*DPolLeg2;
       PertJ:=PertJ+ra*rb;
       PertZ:=PertZ+ra*DPolLeg2;
       PolLeg0:=PolLeg1; { new values for recurrent process }
       PolLeg1:=PolLeg2;
       DPolLeg1:=DPolLeg2;
     End;                          
                                   { -fm_p/r^3*(x,y,z) the main force }
   fmr:=PlanetFMain*rrr*(1-PertJ); { the main force with part of zonal }
                                   { global values ungiamod }
   { the vector of force reffered to the equator of the Planet }
   For i:=1 To 3 Do  FSatEq[i]:=-fmr*PosSat[i]; { main and part zonal }
   FSatEq[3]:=FSatEq[3]-PlanetFMain*rr*PertZ; { to add to z cooredinate }

   { the vector of force reffered to the fixed equator of the Earth }
   For i:=1 To 3 Do
     Begin
       force[i]:=0.0;    { initial nullo forces }
       For j:=1 To 3 Do  { to rotate fro each coordinates }
         force[i]:=force[i]+RotMatr[j,i]*FSatEq[j]; { to J2000 equator }
     End;

 end ; { the forces due to zonal harmonics for our satellite }

 { to calculate the Forces due to
   the gravitational field of the main planet and all other planets }

 Procedure ClcRightHandCoor ( ut   : Extended ;   { time moment julian day }
                              rc   : TVect3 ;     { planetocentric J2000 }
                              Var                 { from all planets }
                                force  : TVect3 ) ; { our forces }
 Var
   td,dt,tb  : Extended ; { terrestrial time TT barycentric time TDB }
 Begin

   TransUTCtoTDT(ut,td,dt); { unit unfortim terrestrial time TT }
   TransTDTtoTDB(td,tb);    { barycentric time TDB Teph for ephemeris }
   { attraction from the main planet unsagvar PlanetMain global }
   ForPlanetMain(tb,rc,force); { from the main planet with zonal harmonics }
                               { the forces }
   PlanetsForces(tb,rc,force); { due to the Sun only }

 End;

 END.