{~ Changing the equatorial coordinates due to precession of the Earth axis ~}

UNIT Precesse;
(*  Copyright (C) 1999 Jan Hollan and 1991 Petr Pravec, ;
  by "program" this unit is further meant.

    This program is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 2 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program; if not, write to the Free Software
    Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
*)

INTERFACE
Function PrecInit(From_Year,To_Year:real):boolean;
                                          {false if From_Year=To_Year}
Procedure Precess(var alfa,delta: real);

(*PRECESSING THE COORDINATE SYSTEM from one equinox to another
                     by Petr Pravec, 1991.
 The accuracy tested on data of the Bright Star Catalogue (conversion
 from B1900 to J2000) and GC (B1950 to J2000) gives standard deviation
 not exceeding one second (of the odrer of 0.0002 degree).
*)

{to a unit converted,
 backward transformation vector->RA,Decl adapted to resist division by zero
 and another minor changes
  (why to make large changes, when it works well)
                                                   by Jan Hollan, 1994}

implementation
const
  mtxdim = 10;
  {lgvdim = 55;    =mtxdim*(mtxdim+1)div2}
  {sqeps = 2E-6;
  ln10 = 2.302585093;}
  deginrad=57.295779513;
type
  {cardinal   = 0..maxint;
  str13      = string[13];
  }
  vector     = array [1..mtxdim] of real;
  intvector  = array [1..mtxdim] of integer;
  {longvector = array [1..lgvdim] of real;
  ilngvector = array [1..lgvdim] of integer;
  slngvector = array [1..lgvdim] of str13;
  }
  matrix     = array [1..mtxdim,1..mtxdim] of real;

var i,j : integer;
    condA : real;
    pm1,pm2 : matrix;
    IxP : intvector;

procedure PRECMAT (    dt   : real;
                   var pmat : matrix);

var pa1,pa2,pa3 : real;

begin
  dt:=dt/100;
  pa1:=(0.6406161+(0.0000839+0.0000050*dt)*dt)*dt/deginrad;
  pa2:=(0.6406161+(0.0003041+0.0000051*dt)*dt)*dt/deginrad;
  pa3:=(0.5567530-(0.0001185+0.0000116*dt)*dt)*dt/deginrad;
  pmat[1,1]:=cos(pa1)*cos(pa2)*cos(pa3)-sin(pa1)*sin(pa2);
  pmat[1,2]:=-sin(pa1)*cos(pa2)*cos(pa3)-cos(pa1)*sin(pa2);
  pmat[1,3]:=-cos(pa2)*sin(pa3);
  pmat[2,1]:=cos(pa1)*sin(pa2)*cos(pa3)+sin(pa1)*cos(pa2);
  pmat[2,2]:=-sin(pa1)*sin(pa2)*cos(pa3)+cos(pa1)*cos(pa2);
  pmat[2,3]:=-sin(pa2)*sin(pa3);
  pmat[3,1]:=cos(pa1)*sin(pa3);
  pmat[3,2]:=-sin(pa1)*sin(pa3);
  pmat[3,3]:=cos(pa3){;
  writeln(pa1*deginrad);
  writeln(pa2*deginrad);
  writeln(pa3*deginrad);
  for i:=1 to 3 do
    for j:=1 to 3 do writeln(pmat[i,j]);}
end;




PROCEDURE Solve  ( dimA :integer; VAR A       :matrix;
                     VAR b:vector;  VAR IndxPvt :intvector);
VAR
  k,m :integer;
  t     :real;
BEGIN  { Solve }
  IF dimA=1
    THEN      { Matice 1x1 }
      b[1]:=b[1]/A[1,1]
    ELSE      { dimA>1 }
      BEGIN
        FOR k:=1 TO dimA-1 DO
          BEGIN
            m:=IndxPvt[k];
            t:=b[m]; b[m]:=b[k]; b[k]:=t;
            FOR i:=k+1 TO dimA DO
              b[i]:=b[i]+A[i,k]*t;
          END ;
                                        { Backward substitution }
        FOR k:=dimA DOWNTO 1 DO
          BEGIN
            b[k]:=b[k]/A[k,k]; t:=-b[k];
            FOR i:=1 TO k-1 DO b[i]:=b[i]+A[i,k]*t
          END
      END

END   { Solve };


PROCEDURE Decomp ( dimA :integer;  VAR A :matrix;
                     VAR condA :real; VAR IndxPvt :intvector);
VAR
  e,t,normaA,normaY,normaZ :real;
  k,kplus1,m           :integer;
  work                 :vector;

BEGIN { Decomp }
  IndxPvt[dimA]:=1;
  IF dimA=1 THEN                                 { Matrix 1x1 }
    BEGIN IF A[1,1]=0 THEN condA:=1E32 END
  ELSE
  BEGIN    { dimA>1 }
   normaA:=0;                               { Norm of A }
   FOR j:=1 TO dimA DO
    BEGIN
      t:=0;
      FOR i:=1 TO dimA DO t:=t+abs(A[i,j]);
      IF t>normaA THEN normaA:=t
    END;
      { Gauss elimination with partial choose of the main element}
   FOR k:=1 to dimA-1 DO
    BEGIN
      kplus1:=k+1; m:=k;                 { Main element }
      FOR i:=kplus1 TO dimA DO
        IF abs(A[i,k])>abs(A[m,k]) THEN m:=i;
      IndxPvt[k]:=m;
      IF m<>k THEN IndxPvt[dimA]:=-IndxPvt[dimA];
      t:=A[m,k]; A[m,k]:=A[k,k]; A[k,k]:=t;
      IF t<>0 THEN
        BEGIN
          FOR i:=kplus1 TO dimA DO   { multiplicatr }
            A[i,k]:=-A[i,k]/t;
          FOR j:=kplus1 TO dimA DO
            BEGIN           { Vymena a eliminace po sloupcich }
              t:=A[m,j]; A[m,j]:=A[k,j]; A[k,j]:=t;
              FOR i:=kplus1 TO dimA DO A[i,j]:=A[i,j]+A[i,k]*t;
            END
        END   { t<>0 }
    END    { cycle k };
                                       { condA estimate }
  condA:=1;
  FOR k:=1 TO dimA DO IF A[k,k]=0 THEN condA:=1.0E32;
  IF condA=1 THEN
    BEGIN    { regular matrix }
     FOR k:=1 TO dimA DO
      BEGIN
       t:=0;
       IF k>1 THEN
         FOR i:=1 TO k-1 DO t:=t+A[i,k]*work[i]; {work gets value at k=1}
       e:=1; IF t<0 THEN e:=-1;
       work[k]:=-(e+t)/A[k,k]
      END;
     FOR k:=dimA-1 DOWNTO 1 DO
      BEGIN
       t:=0;
       FOR i:=k+1 TO dimA DO t:=t+A[i,k]*work[k];
       work[k]:=t; m:=IndxPvt[k];
       IF m<>k THEN
         BEGIN t:=work[m]; work[m]:=work[k]; work[k]:=t END;
      END;
     normaY:=0;
     FOR i:=1 TO dimA DO normaY:=normaY+abs(work[i]);
               { solving A*z=y system by procedure Solve }
     Solve(dimA,A,work,IndxPvt);
     normaZ:=0;
     FOR i:=1 TO dimA DO normaZ:=normaZ+abs(work[i]);
                             { estimate of "condition number" of condA }
     condA:=normaA*normaZ/normaY;
    END    { Regular A }
  END      { dimA>1 }

END   { Decomp };


Function PrecInit(From_Year,To_Year:real):boolean;
begin
     From_Year:=From_Year - 2000.0 ;
     To_Year:=To_Year - 2000.0 ;
     if From_Year<>0 then PRECMAT(From_Year,pm1)
      else
        for i:=1 to 3 do
          for j:=1 to 3 do
            if i=j then pm1[i,j]:=1 else pm1[i,j]:=0;
     DECOMP(3,pm1,condA,IxP);
     if To_Year<>0 then PRECMAT(To_Year,pm2)
      else
        for i:=1 to 3 do
          for j:=1 to 3 do
            if i=j then pm2[i,j]:=1 else pm2[i,j]:=0;
      if From_Year=To_Year then
       PrecInit:=false
      else
       PrecInit:=true;
end;


procedure Precess  ( var alfa, delta : real);
                     {alfa,delta : real;
                        pm1,pm2 : matrix;
                        IxP : intvector);}

var x1,x2 : vector;

begin
  alfa:=alfa/deginrad; delta:=delta/deginrad;
  x1[1]:=cos(delta)*cos(alfa);
  x1[2]:=cos(delta)*sin(alfa);
  x1[3]:=sin(delta);
  for i:=1 to 3 do x2[i]:=0;
  SOLVE(3,pm1,x1,IxP);
  for i:=1 to 3 do
    for j:=1 to 3 do x2[i]:=x2[i]+pm2[i,j]*x1[j];
  if x2[3]=1 then
    begin
     delta:=90; alfa:=360
    end
   else
    if x2[3]=-1 then
      begin
       delta:=-90; alfa:=360
      end
     else
      begin
       IF X2[3]<0.7 THEN
         delta:=deginrad*arctan(x2[3]/sqrt(SQR(X2[1])+sqr(X2[2])))
       ELSE
         DELTA:=90-deginrad*arctan(sqrt(SQR(X2[1])+sqr(X2[2]))/x2[3])  ;
       IF ABS(X2[2])<ABS(X2[1]) THEN        {AROUND 0 OR 180}
        BEGIN
         ALFA:=DEGINRAD*ARCTAN(X2[2]/X2[1]);
         if x2[1]<0 then alfa:=alfa+180;
         IF (X2[2]<0) AND (X2[1]>=0) THEN ALFA:=ALFA+360
        END
       else                          {AROUND 90 OR 270}
        begin
         alfa:=90-deginrad*arctan(X2[1]/X2[2]);
         if x2[2]<0 then alfa:=180+alfa;
        end
      end;
end;

end.