unit ungiarom ;
interface
uses
  unsagcot ; { TypeDim33 }
var
  EquatorMatRot  : TypeDim33 ; { global }
{ from fixed equator of the Earth to equator of the Planet }
Procedure ToClcMatrOfRotation ( PlanetMain     : byte ;
                                TinJD          : Extended ;
                                Var  MatrRot   : TypeDim33 ) ;
implementation
uses
  unephjpl ; { global const NJupiter for example }
{ from fixed equator of the Earth to equator of the Planet }
Procedure ToClcMatrOfRotation ( PlanetMain     : byte ;
                                TinJD          : Extended ;
                                Var  MatrRot   : TypeDim33 ) ;
Var
  i,j,k              : Byte ;
  dt,dn,ap,dp,
  ax,az,sx,sz,cx,cz  : Extended ;
  xm,zm              : TypeDim33 ;
Begin
  dt:=(TinJD-jd2000)/36525.0; { unsagcot const }
  ap:=0.0;
  dp:=90.0; { default value }
  If  PlanetMain = NJupiter  Then   { unephjpl global }
    Begin
      ap:=GraRad*(268.05-0.009*dt); { unsagcot constants }
      dp:=GraRad*( 64.49+0.003*dt);
    End;
  If  PlanetMain = NSaturn  Then   { unephjpl global }
    Begin
      ap:=GraRad*(40.589-0.036*dt); { unsagcot constants }
      dp:=GraRad*(83.537-0.004*dt);
    End;
  If  PlanetMain = NUranus  Then   { unephjpl global }
    Begin
      ap:=GraRad*257.311; { unsagcot constants }
      dp:=GraRad*(-15.175);
    End;
  If  PlanetMain = NNeptune  Then
    Begin
      dn:=GraRad*(357.85+52.316*dt);
      ap:=GraRad*(299.36+0.70*Sin(dn));
      dp:=GraRad*( 43.46-0.51*Cos(dn));
    End;
  az:=ap-HalfPi; { to descending node }
  ax:=-(HalfPi-dp); { to the rotation axis of the planet }
  { to rotate from Earth equator to planet equator }
  sz:=Sin(az); cz:=Cos(az);  sx:=Sin(ax); cx:=Cos(ax);
  zm[1,1]:=+cz; zm[1,2]:=+sz; zm[1,3]:=0; { Rz(AlfaPole-pi/2)  }
  zm[2,1]:=-sz; zm[2,2]:=+cz; zm[2,3]:=0; { to descending node }
  zm[3,1]:=0;   zm[3,2]:=0;   zm[3,3]:=1;
  xm[1,1]:=1;   xm[1,2]:=0;   xm[1,3]:=0; { Rx(-(pi/2-DeltaPole)) }
  xm[2,1]:=0;   xm[2,2]:=+cx; xm[2,3]:=+sx; { to the rotation }
  xm[3,1]:=0;   xm[3,2]:=-sx; xm[3,3]:=+cx; { axis of the planet }
  For i:=1 To 3 Do
    For j:=1 To 3 Do
      Begin
        MatrRot[i,j]:=0;
        For k:=1 to 3 Do
          MatrRot[i,j]:=MatrRot[i,j]+xm[i,k]*zm[k,j];
      End;
End;
end.