📄 stastrop.pas
字号:
B := B + (A * sqr(T0) / 1.0E+8);
A := 276 * cos(0.595 + 6283.076*T0)
+ 17 * cos(3.14)
+ 4 * cos(0.12 + 12566.15*T0);
B := B + (A * sqr(T0) * T0 / 1.0E+8);
{solar radius vector (astronomical units)}
RV := 100013989
+ 1670700 * cos(3.0984635 + 6283.07585*T0)
+ 13956 * cos(3.05525 + 12566.15170*T0)
+ 3084 * cos(5.1985 + 77713.7715*T0)
+ 1628 * cos(1.1739 + 5753.3849*T0)
+ 1576 * cos(2.8649 + 7860.4194*T0)
+ 925 * cos(5.453 + 11506.770*T0)
+ 542 * cos(4.564 + 3930.210*T0)
+ 472 * cos(3.661 + 5884.927*T0)
+ 346 * cos(0.964 + 5507.553*T0)
+ 329 * cos(5.900 + 5223.694*T0)
+ 307 * cos(0.299 + 5573.143*T0)
+ 243 * cos(4.273 + 11790.629*T0)
+ 212 * cos(5.847 + 1577.344*T0)
+ 186 * cos(5.022 + 10977.079*T0)
+ 175 * cos(3.012 + 18849.228*T0)
+ 110 * cos(5.055 + 5486.778*T0)
+ 98 * cos(0.89 + 6069.78*T0)
+ 86 * cos(5.69 + 15720.84*T0)
+ 86 * cos(1.27 +161000.69*T0)
+ 65 * cos(0.27 + 17260.15*T0)
+ 63 * cos(0.92 + 529.69*T0)
+ 57 * cos(2.01 + 83996.85*T0)
+ 56 * cos(5.24 + 71430.70*T0)
+ 49 * cos(3.25 + 2544.31*T0)
+ 47 * cos(2.58 + 775.52*T0)
+ 45 * cos(5.54 + 9437.76*T0)
+ 43 * cos(6.01 + 6275.96*T0)
+ 39 * cos(5.36 + 4694.00*T0)
+ 38 * cos(2.39 + 8827.39*T0)
+ 37 * cos(0.83 + 19651.05*T0)
+ 37 * cos(4.90 + 12139.55*T0)
+ 36 * cos(1.67 + 12036.46*T0)
+ 35 * cos(1.84 + 2942.46*T0)
+ 33 * cos(0.24 + 7084.90*T0)
+ 32 * cos(0.18 + 5088.63*T0)
+ 32 * cos(1.78 + 398.15*T0)
+ 28 * cos(1.21 + 6286.60*T0)
+ 28 * cos(1.90 + 6279.55*T0)
+ 26 * cos(4.59 + 10447.39*T0);
RV := RV / 1.0E+8;
A := 103019 * cos(1.107490 + 6283.075850*T0)
+ 1721 * cos(1.0644 + 12566.1517*T0)
+ 702 * cos(3.142)
+ 32 * cos(1.02 + 18849.23*T0)
+ 31 * cos(2.84 + 5507.55*T0)
+ 25 * cos(1.32 + 5223.69*T0)
+ 18 * cos(1.42 + 1577.34*T0)
+ 10 * cos(5.91 + 10977.08*T0)
+ 9 * cos(1.42 + 6275.96*T0)
+ 9 * cos(0.27 + 5486.78*T0);
RV := RV + (A * T0 / 1.0E+8);
A := 4359 * cos(5.7846 + 6283.0758*T0)
+ 124 * cos(5.579 + 12566.152*T0)
+ 12 * cos(3.14)
+ 9 * cos(3.63 + 77713.77*T0)
+ 6 * cos(1.87 + 5573.14*T0)
+ 3 * cos(5.47 + 18849.23*T0);
RV := RV + (A * sqr(T0) / 1.0E+8);
L := (L + PI);
L := Frac(L / 2.0 / PI) * 2.0 * Pi;
if L < 0 then
L := L + (2.0*PI);
B := -B;
TX := RV * cos(B) * cos(L);
TY := RV * cos(B) * sin(L);
TZ := RV * sin(B);
Result.X := TX + 4.40360E-7 * TY - 1.90919E-7 * TZ;
Result.Y := -4.79966E-7 * TX + 0.917482137087 * TY - 0.397776982902 * TZ;
Result.Z := 0.397776982902 * TY + 0.917482137087 * TZ;
end;
{--------------------------------------------------------------------------}
function EclipticToRectangular(Longitude, Latitude,
RadiusVector : Double) : TStRectangularCord;
var
var1,
var2,
var3 : Double;
begin
var1 := RadiusVector * cos(Longitude) * cos(Latitude);
var2 := RadiusVector * sin(Longitude) * cos(Latitude);
var3 := RadiusVector * sin(Latitude);
Result.X := var1;
Result.Y := var2 * cos(OB2000) - var3 * sin(OB2000);
Result.Z := var2 * sin(OB2000) + var3 * cos(OB2000);
end;
{--------------------------------------------------------------------------}
function RADec(Planet, Sun : TStRectangularCord;
ComputeElong : Boolean) : TStPlanetsRec;
var
var1,
var2,
var3,
var4,
var5 : Double;
begin
FillChar(Result, SizeOf(TStPlanetsRec), #0);
var1 := Sun.X + Planet.X;
var2 := Sun.Y + Planet.Y;
var3 := Sun.Z + Planet.Z;
var4 := arctan(var2/var1);
var4 := RealAngle(var2, var1, var4) * radcor;
var5 := sqrt(sqr(var1) + sqr(var2) + sqr(var3));
var3 := StInvsin(var3/var5) * radcor;
Result.RA := var4;
Result.DC := var3;
var4 := Result.RA / radcor;
var3 := Result.DC / radcor;
if (ComputeElong) then begin
var1 := sin(SunEQ.DC/radcor) * sin(var3);
var2 := cos(SunEQ.DC/radcor) * cos(var3) * cos(SunEQ.RA/radcor - var4);
Result.Elong := StInvcos(var1+var2) * radcor;
end;
end;
{--------------------------------------------------------------------------}
function MercuryPosition(JD : Double) : TStPlanetsRec;
begin
PlanEC := ComputeMercury(JD);
PlanRC := EclipticToRectangular(PlanEC.L0, PlanEC.B0, PlanEC.R0);
Result := RADec(PlanRC, SunRC, True);
end;
{--------------------------------------------------------------------------}
function VenusPosition(JD : Double) : TStPlanetsRec;
begin
PlanEC := ComputeVenus(JD);
PlanRC := EclipticToRectangular(PlanEC.L0, PlanEC.B0, PlanEC.R0);
Result := RADec(PlanRC, SunRC, True);
end;
{--------------------------------------------------------------------------}
function MarsPosition(JD : Double) : TStPlanetsRec;
begin
PlanEC := ComputeMars(JD);
PlanRC := EclipticToRectangular(PlanEC.L0, PlanEC.B0, PlanEC.R0);
Result := RADec(PlanRC, SunRC, True);
end;
{--------------------------------------------------------------------------}
function JupiterPosition(JD : Double) : TStPlanetsRec;
begin
PlanEC := ComputeJupiter(JD);
PlanRC := EclipticToRectangular(PlanEC.L0, PlanEC.B0, PlanEC.R0);
Result := RADec(PlanRC, SunRC, True);
end;
{--------------------------------------------------------------------------}
function SaturnPosition(JD : Double) : TStPlanetsRec;
begin
PlanEC := ComputeSaturn(JD);
PlanRC := EclipticToRectangular(PlanEC.L0, PlanEC.B0, PlanEC.R0);
Result := RADec(PlanRC, SunRC, True);
end;
{--------------------------------------------------------------------------}
function UranusPosition(JD : Double) : TStPlanetsRec;
begin
PlanEC := ComputeUranus(JD);
PlanRC := EclipticToRectangular(PlanEC.L0, PlanEC.B0, PlanEC.R0);
Result := RADec(PlanRC, SunRC, True);
end;
{--------------------------------------------------------------------------}
function NeptunePosition(JD : Double) : TStPlanetsRec;
begin
PlanEC := ComputeNeptune(JD);
PlanRC := EclipticToRectangular(PlanEC.L0, PlanEC.B0, PlanEC.R0);
Result := RADec(PlanRC, SunRC, True);
end;
{--------------------------------------------------------------------------}
function PlutoPosition(JD : Double) : TStPlanetsRec;
begin
PlanEC := ComputePluto(JD);
PlanRC := EclipticToRectangular(PlanEC.L0, PlanEC.B0, PlanEC.R0);
Result := RADec(PlanRC, SunRC, True);
end;
{--------------------------------------------------------------------------}
procedure PlanetsPos(JD : Double; var PA : TStPlanetsArray);
var
I : Integer;
Sun : TStRectangularCord;
begin
{find Sun's Rectangular Coordinates}
SunRC := SunofDate(JD);
FillChar(SunEQ, SizeOf(TStPlanetsRec), #0);
FillChar(Sun, SizeOf(TStRectangularCord), #0);
{find Sun's RA/Dec}
SunEQ := RADec(SunRC, Sun, False);
PA[1] := PlutoPosition(JD);
{find RA/Dec of each planet}
for I := 1 to 8 do begin
case I of
1 : PA[I] := MercuryPosition(JD);
2 : PA[I] := VenusPosition(JD);
3 : PA[I] := MarsPosition(JD);
4 : PA[I] := JupiterPosition(JD);
5 : PA[I] := SaturnPosition(JD);
6 : PA[I] := UranusPosition(JD);
7 : PA[I] := NeptunePosition(JD);
8 : PA[I] := PlutoPosition(JD);
end;
end;
end;
end.
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -