unit uVsop;

 {$i ah_def.inc }

{  Calculates the planetary heliocentric coordinates according to the
   VSOP87 theory. Calculations according to chapter 32 (31) of Meeus. }

(* $define meeus *)   { Only use the accuracy as in the Meeus book }

(*$ifdef delphi_1 *)
  (*$define meeus *) { Otherwise the code segment will be too small }
(*$endif *)

(*@/// interface *)
interface

(*@/// uses *)
uses
  uAHMath,
  sysutils;
(*@\\\000000020B*)

type
  (*@/// TVSOPEntry=record *)
  TVSOPEntry=record
    A,B,C: extended;
    end;
  (*@\\\*)
  TVSOPCalcFunc = function (nr,index: integer):TVSOPEntry of object;
  (*@/// TVSOP=class(TObject) *)
  TVSOP=class(TObject)
  protected
    FDate: TDateTime;
    function LongitudeFactor(nr,index: integer):TVSOPEntry;  VIRTUAL; abstract;
    function LatitudeFactor(nr,index: integer):TVSOPEntry;  VIRTUAL; abstract;
    function RadiusFactor(nr,index: integer):TVSOPEntry;  VIRTUAL; abstract;
    function CalcLongitude:extended;
    function CalcLatitude:extended;
    function CalcRadius:extended;
    function Calc(factor: TVSOPCalcFunc):extended;
    procedure SetDate(value: TDateTime);
    function Tau:extended;
  public
    procedure DynamicToFK5(var longitude,latitude: extended);
    property Longitude:extended read CalcLongitude;
    property Latitude:extended read CalcLatitude;
    property Radius:extended read CalcRadius;
    property Date:TDateTime write SetDate;
    end;
  (*@\\\0000000E01*)
  TCVSOP=class of TVSOP;
  (*@/// TVSOPEarth=class(TVSOP) *)
  TVSOPEarth=class(TVSOP)
  protected
    function LongitudeFactor(nr,index: integer):TVSOPEntry;  override;
    function LatitudeFactor(nr,index: integer):TVSOPEntry;  override;
    function RadiusFactor(nr,index: integer):TVSOPEntry;  override;
    end;
  (*@\\\0000000607*)
  (*@/// TVSOPJupiter=class(TVSOP) *)
  TVSOPJupiter=class(TVSOP)
  protected
    function LongitudeFactor(nr,index: integer):TVSOPEntry;  override;
    function LatitudeFactor(nr,index: integer):TVSOPEntry;  override;
    function RadiusFactor(nr,index: integer):TVSOPEntry;  override;
    end;
  (*@\\\0000000607*)

procedure earth_coord(date:TdateTime; var l,b,r: extended);
procedure jupiter_coord(date:TdateTime; var l,b,r: extended);
(*@\\\0000000301*)
(*@/// implementation *)
implementation

uses
  uMoon;

(*$ifdef delphi_ge_3 *)
var
(*$else *)
const
(*$endif *)
  datetime_2000_01_01: extended = 0;

(*@/// procedure calc_coord(date: TDateTime; obj_class: TCVSOP; var l,b,r: extended); *)
procedure calc_coord(date: TDateTime; obj_class: TCVSOP; var l,b,r: extended);
var
  obj: TVSOP;
begin
  obj:=NIL;
  try
    obj:=obj_class.Create;
    obj.date:=date;
    r:=obj.radius;
    l:=obj.longitude;
    b:=obj.latitude;
    obj.DynamicToFK5(l,b);
  finally
    obj.free;
    end;
  l:=put_in_360(rad2deg(l));  (* rad -> degree *)
  b:=rad2deg(b);
  end;
(*@\\\0000001111*)

(*@/// procedure earth_coord(date:TdateTime; var l,b,r: extended); *)
procedure earth_coord(date:TdateTime; var l,b,r: extended);
begin
  calc_coord(date,TVSOPEarth,l,b,r);
  end;
(*@\\\0000000116*)
(*@/// procedure jupiter_coord(date:TdateTime; var l,b,r: extended); *)
procedure jupiter_coord(date:TdateTime; var l,b,r: extended);
begin
  calc_coord(date,TVSOPJupiter,l,b,r);
  end;
(*@\\\000000031C*)

(*@/// class TVSOP *)
(*@/// function TVSOP.CalcLongitude:extended; *)
function TVSOP.CalcLongitude:extended;
begin
  result:=calc(Longitudefactor);
  end;
(*@\\\0000000401*)
(*@/// function TVSOP.CalcLatitude:extended; *)
function TVSOP.CalcLatitude:extended;
begin
  result:=calc(Latitudefactor);
  end;
(*@\\\000000031F*)
(*@/// function TVSOP.CalcRadius:extended; *)
function TVSOP.CalcRadius:extended;
begin
  result:=calc(radiusfactor);
  end;
(*@\\\000000031D*)
(*@/// procedure TVSOP.SetDate(value: TDateTime); *)
procedure TVSOP.SetDate(value: TDateTime);
begin
  FDate:=value;
  end;
(*@\\\*)
(*@/// function TVSOP.Tau:extended; *)
function TVSOP.Tau:extended;
begin
  result:=(FDate-datetime_2000_01_01-0.5)/365250.0;
  end;
(*@\\\0000000301*)
(*@/// function TVSOP.Calc(factor: TVSOPCalcFunc):extended; *)
function TVSOP.Calc(factor: TVSOPCalcFunc):extended;
var
  t: extended;
  current: extended;
  r: array[0..5] of extended;
  i,j: integer;
begin
  t:=Tau;
  for j:=0 to 5 do begin
    r[j]:=0;
    i:=0;
    repeat
      WITH Factor(i,j) do
        current:=a*cos(b+c*t);
      r[j]:=r[j]+current;
      inc(i);
    until current=0;
    end;
  result:=(r[0]+t*(r[1]+t*(r[2]+t*(r[3]+t*(r[4]+t*r[5])))))*1e-8;
  end;
(*@\\\0000000E17*)
(*@/// procedure TVSOP.DynamicToFK5(var longitude,latitude: extended); *)
procedure TVSOP.DynamicToFK5(var longitude,latitude: extended);
var
  lprime,t: extended;
  delta_l, delta_b: extended;
begin
  t:=10*tau;
  lprime:=longitude+deg2rad(-1.397-0.00031*t)*t;
  delta_l:=-deg2rad(0.09033/3600)+deg2rad(0.03916/3600)*(cos(lprime)+sin(lprime))*tan(latitude);
  delta_b:=deg2rad(0.03916/3600)*(cos(lprime)-sin(lprime));
  longitude:=longitude+delta_l;
  latitude:=latitude+delta_b;
  end;
(*@\\\*)
(*@\\\0000000226*)
(*@/// class TVSOPEarth *)
(*@/// function TVSOPEarth.RadiusFactor(nr,index: integer):TVSOPEntry; *)
function TVSOPEarth.RadiusFactor(nr,index: integer):TVSOPEntry;
const
  (*@/// vsop87_ear_r0:array[0..525,0..2] of extended = (...); *)
  (*$ifdef meeus *)
  vsop87_ear_r0:array[0.. 39,0..2] of extended = (
  (*$else *)
  vsop87_ear_r0:array[0..525,0..2] of extended = (
  (*$endif *)
  { 4330    1 }  ( 100013988.799, 0.00000000000,      0.00000000000 ),
  { 4330    2 }  (   1670699.626, 3.09846350771,   6283.07584999140 ),
  { 4330    3 }  (     13956.023, 3.05524609620,  12566.15169998280 ),
  { 4330    4 }  (      3083.720, 5.19846674381,  77713.77146812050 ),
  { 4330    5 }  (      1628.461, 1.17387749012,   5753.38488489680 ),
  { 4330    6 }  (      1575.568, 2.84685245825,   7860.41939243920 ),
  { 4330    7 }  (       924.799, 5.45292234084,  11506.76976979360 ),
  { 4330    8 }  (       542.444, 4.56409149777,   3930.20969621960 ),
  { 4330    9 }  (       472.110, 3.66100022149,   5884.92684658320 ),
  { 4330   10 }  (       328.780, 5.89983646482,   5223.69391980220 ),
  { 4330   11 }  (       345.983, 0.96368617687,   5507.55323866740 ),
  { 4330   12 }  (       306.784, 0.29867139512,   5573.14280143310 ),
  { 4330   13 }  (       174.844, 3.01193636534,  18849.22754997420 ),