БЭСМ before…

Простые модули

Следующая программная единица содержит исходные тексты простых вычислительных процедур.

UNIT UnForFun ; // for the simple procedure for calculating

INTERFACE

Uses

UnConTyp ;

Function DATAN2 ( SinAngle , CosAngle : Extended ) : Extended ;

Procedure ArcTanTwoArg ( s , c : Extended ; Var a : Extended ) ;

Function ATanDegree ( s , c : Extended ) : Extended ;

procedure PrDegrAngle(Angle : Extended;

Var AngleSign : ShortString;

Var IntDegree,IntMin : integer;

Var AngleSec : Extended);

procedure PrHourAngle(Angle : Extended;

Var AngleSign : ShortString;

Var IntHour,IntMin : integer;

Var AngleSec : Extended);

Function Log10 ( x : Extended ) : Extended ;

Function VectorModul ( a : TVect3 ) : Extended ;

Procedure MatrMultVector ( p : TMatr33 ;

a : TVect3 ;

Var b : TVect3 ) ;

Procedure UMatrMultVect ( p : TMatr33 ;

a : TVect3 ;

Var b : TVect3 ) ;

Procedure MatrMultMatr ( p , q : TMatr33 ;

Var r : TMatr33 ) ;

Procedure SimpleRotMatr ( m : Byte ; a : Extended ;

Var r : TMatr33 ) ;

IMPLEMENTATION

{ DATAN2 is a function

to obtain Arctan for the Angle in interval

from -(1/2)*Pi to +(3/2)*Pi in rad }

Function DATAN2 ( SinAngle , CosAngle : Extended ) : Extended ;

Var

a,c,s,p : Extended;

Begin

c:=CosAngle; s:=SinAngle;

if c <> 0 then Begin

p:=s/c; a:=ArcTan(p);

if c < 0 then a:=a+Pi; End

else if s > 0 then a:=PiTwo else a:=-PiTwo;

DATAN2:=a;

End;

Procedure ArcTanTwoArg ( s , c : Extended ; Var a : Extended ) ;

Begin

a:=DATAN2(s,c);

End;

{ ATanDegree is a function to obtain

Arctan for the Angle in interval from 0 to 360 in degree }

Function ATanDegree ( s , c : Extended ) : Extended ;

Var

a,p : Extended;

Begin

If c <> 0

Then

Begin

p:=s/c;

a:=RadGra*ArcTan(p);

If c < 0 Then a:=a+180;

End

Else

If s > 0 Then a:=90 Else a:=270;

If a < 0 Then a:=360+a;

ATanDegree:=a;

End;

{ PrDegrAngle is quite simple procedure

to change a value for Angle in radians

to the civil form : Sign Degree Minute Second . }

procedure PrDegrAngle(Angle : Extended;

Var AngleSign : ShortString;

Var IntDegree,IntMin : integer;

Var AngleSec : Extended);

Var

r : Extended;

Begin

if Angle < 0 then AngleSign:='-' else AngleSign:='+';

r:=RadGra*Abs(Angle); IntDegree:=Trunc(r);

r:=60*Frac(r); IntMin:=Trunc(r); AngleSec:=60*Frac(r);

if AngleSec > 59.999 then begin AngleSec:=0; IntMin:=IntMin+1; end;

if IntMin > 59 then

begin IntMin:=IntMin-60; IntDegree:=IntDegree+1; end;

End;

{ PrHourAngle is quite simple procedure

to change a value for Angle in radians

to the civil form : Sign Hour Minute Second . }

procedure PrHourAngle(Angle : Extended;

Var AngleSign : ShortString;

Var IntHour,IntMin : integer;

Var AngleSec : Extended);

Var

r : Extended;

Begin

if Angle < 0 then AngleSign:='-' else AngleSign:='+';

r:=RadGra*Abs(Angle)/15; IntHour:=Trunc(r);

r:=60*Frac(r); IntMin:=Trunc(r); AngleSec:=60*Frac(r);

if AngleSec > 59.999 then begin AngleSec:=0; IntMin:=IntMin+1; end;

if IntMin > 59 then

begin IntMin:=IntMin-60; IntHour:=IntHour+1; end;

End;

Function Log10 ( x : Extended ) : Extended ;

Begin

Log10:=Ln(x)/Ln(10); // from the interesting book

End;

Function VectorModul ( a : TVect3 ) : Extended ;

Begin

VectorModul:=Sqrt(Sqr(a[1])+Sqr(a[2])+Sqr(a[3]));

End;

Procedure MatrMultVector ( p : TMatr33 ;

a : TVect3 ;

Var b : TVect3 ) ;

Var

i,j : Byte ;

Begin

For i:=1 To 3 Do

Begin

b[i]:=0;

For j:=1 To 3 Do

b[i]:=b[i]+p[i,j]*a[j];

End;

End;

Procedure UMatrMultVect ( p : TMatr33 ; // to transform matrix

a : TVect3 ;

Var b : TVect3 ) ;

Var

i,j : Byte ;

Begin

For i:=1 To 3 Do

Begin

b[i]:=0;

For j:=1 To 3 Do

b[i]:=b[i]+p[j,i]*a[j]; // to transform matrix

End;

End;

Procedure MatrMultMatr ( p , q : TMatr33 ;

Var r : TMatr33 ) ;

Var

i,j,k : Byte ;

Begin

For i:=1 To 3 Do

For j:=1 To 3 Do

Begin

r[i,j]:=0;

For k:=1 To 3 Do

r[i,j]:=r[i,j]+p[i,k]*q[k,j];

End;

End;

{ angle a in degree, m = 1 , 2 , 3 for X Y Z axis }

Procedure SimpleRotMatr ( m : Byte ; a : Extended ;

Var r : TMatr33 ) ;

Var

i,j : Byte ;

b,s,c : Extended ;

Begin

for i:=1 to 3 do

for j:=1 to 3 do

if i = j then r[i,j]:=1.0 else r[i,j]:=0.0; // unit matrix

b:=GraRad*a;

s:=Sin(b);

c:=Cos(b);

If m = 1

Then

Begin

r[2,2]:=+c;

r[2,3]:=+s;

r[3,2]:=-s;

r[3,3]:=+c;

Exit ;

End;

If m = 2

Then

Begin

r[1,1]:=+c;

r[1,3]:=-s;

r[3,1]:=+s;

r[3,3]:=+c;

Exit ;

End;

If m = 3

Then

Begin

r[1,1]:=+c;

r[1,2]:=+s;

r[2,1]:=-s;

r[2,2]:=+c;

Exit ;

End;

END. // unit UnForFun

Для корректного преобразования строк в числовые данные могут оказаться полезными следующие процедуры и функции.

UNIT UnForStr ; // simple operations with the strings

INTERFACE

Function WordToFixStr ( w : Word ; n : Byte ) : ShortString ;

Function IntegerToFixStr ( w : Integer ; n : Byte ) : ShortString ;

Function FloatToFixStr ( e : Extended ; n , k : Byte ) : ShortString ;

Function BlanksToStr ( s : ShortString ; n : Byte ) : ShortString ;

Procedure ToEraseSomeBlanks ( Var sn : ShortString ) ;

Procedure ToEraseAllBlanks ( Var sn : ShortString ) ;

Procedure ToReplaceBlanks ( Var st : ShortString ) ;

Function WithOutDecimalPoint ( st : ShortString ) : ShortString ;

IMPLEMENTATION

Function WordToFixStr ( w : Word ; n : Byte ) : ShortString ;

Var

i : Byte ;

s : ShortString ;

Begin

Str(w:n,s);

For i:=1 To n Do If s[i] = ' ' Then s[i]:='0';

WordToFixStr:=s;

End;

Function IntegerToFixStr ( w : Integer ; n : Byte ) : ShortString ;

Var

i : Byte ;

s : ShortString ;

Begin

Str(w:n,s);

For i:=1 To n Do If s[i] = ' ' Then s[i]:='0';

IntegerToFixStr:=s;

End;

Function FloatToFixStr ( e : Extended ; n , k : Byte ) : ShortString ;

Var

i,p : Byte ;

s : ShortString ;

Begin

Str(e:n:k,s);

p:=Pos('-',s);

If p > 1

Then

Begin

s[p]:='0';

s[1]:='-';

End;

For i:=1 To n Do

If s[i] = ' '

Then s[i]:='0';

FloatToFixStr:=s;

End;

Function BlanksToStr ( s : ShortString ; n : Byte ) : ShortString ;

Var

a : ShortString ;

l : Byte ;

Begin

a:=s; // our string

l:=Length(a); // length of our sting

If l > n

Then

Delete(a,n+1,l-n) // length is greater than necessary

Else

While l < n Do // length of our string is less than necessary

Begin // to add blanks if it is needed

a:=' '+a; // to add one blank forward

If Length(a) < n Then a:=a+' '; // to add one blank backward

l:=Length(a);

End;

BlanksToStr:=a;

End;

Procedure ToEraseSomeBlanks ( Var sn : ShortString ) ;

Var

i,l : Byte ;

Begin

l:=Length(sn);

For i:=l DownTo 1 Do

If sn[i] = ' '

Then

Delete(sn,i,1)

Else

Exit;

End;

Procedure ToEraseAllBlanks ( Var sn : ShortString ) ;

Var

n : Byte ;

Begin

n:=Pos(' ',sn);

While n > 0 Do

Begin

Delete(sn,n,1);

n:=Pos(' ',sn);

End;

End;

Procedure ToReplaceBlanks ( Var st : ShortString ) ;

Var

i,l,n : Byte ;

Begin

l:=Length(st);

For i:=1 To l Do

If st[i] = ' ' Then st[i]:='0';

n:=Pos('+',st); // may be plus in string

if n > 1

then

begin

st[1]:='+'; // the first character plus

st[n]:='0'; // replace by nullo

end;

n:=Pos('-',st); // may be minus in string

if n > 1

then

begin

st[1]:='-'; // the first character minus

st[n]:='0'; // replace by nullo

end;

End;

Function WithOutDecimalPoint ( st : ShortString ) : ShortString ;

Var

n : Byte ;

s : ShortString ;

Begin

s:=st;

n:=Pos('.',s);

While n > 0 Do

Begin

Delete(s,n,1);

n:=Pos('.',s);

End;

WithOutDecimalPoint:=s;

End;

END. // unit UnForStr

Шкалы времени

Следующие процедуры и функции предназначены для вычисления значений момента времени в системах различных шкал: всемирное координированное время, земное время, эфемеридное время, гринвичское среднее звёздное время.

UNIT UnForTim ;

INTERFACE

Procedure FromUTCtoTT ( JulUTC : Extended;

Var JulTT , DeltaTA : Extended ) ;

Function ToGetTEph ( JulTT : Extended ) : Extended ;

Function ToGetUT1 ( JulUTC , DeltaUT1 : Extended ) : Extended ;

Function ToGetGrMSTime( TinUTC : Extended ) : Extended ;

IMPLEMENTATION

Uses

UnConTyp , // for DifEpoch , JD2000 , JulianC , StrDayWeek

UnForDat ; // for procedure TransDateToJD

{ FromUTCtoTT is a procedure for transition UTC to TT .

Пpевpащение всемиpного кооpдиниpованного вpемени в земное время TT.

Входной паpаметp:

JDutc - юлианская дата, вычисленная по вpемени UTC, не pаньше

1 янваpя 1962 г.

Выходной паpаметp:

JDE - эфемеpидная юлианская дата.

Литеpатуpа: Connaissance des temps 1988, стp. XXVII,

А.А.Коpсунь, Шкалы вpемени, Итоги науки и техники. Астpономия,

1976, т. 12, часть 2, с. 21.

Автоp пpогpаммы пpофессоp К.В.Куимов . }

Procedure FromUTCtoTT ( JulUTC : Extended;

Var JulTT , DeltaTA : Extended ) ;

Var

i : Integer ;

r : Double ;

Const

MaxDTable = 34 ;

Const

TablOne : array[1..7,1..2] of Double =

((35473.00, 31.34), { 1956 1 1 }

(35747.00, 31.56), { 1956 10 1 }

(36112.00, 32.00), { 1957 10 1 }

(36477.00, 32.52), { 1958 10 1 }

(36842.00, 33.00), { 1959 10 1 }

(37208.00, 33.45), { 1960 10 1 }

(37665.00, 33.99)); { 1962 1 1 }

Const

Tabl : array[1..MaxDTable,1..4] of Double =

((37665.0, 34.0298580, 0.0011232, 37665.0), { 1962 01 01 }

(38334.0, 34.1298580, 0.0011232, 37665.0), { 1963 11 01 }

(38395.0, 35.4241300, 0.0012960, 38761.0), { 1964 01 01 }

(38486.0, 35.5241300, 0.0012960, 38761.0), { 1964 04 01 }

(38639.0, 35.6241300, 0.0012960, 38761.0), { 1964 09 01 }

(38761.0, 35.7241300, 0.0012960, 38761.0), { 1965 01 01 }

(38820.0, 35.8241300, 0.0012960, 38761.0), { 1965 03 01 }

(38942.0, 35.9241300, 0.0012960, 38761.0), { 1965 07 01 }

(39004.0, 36.0241300, 0.0012960, 38761.0), { 1965 09 01 }

(39126.0, 36.4971700, 0.0025920, 39126.0), { 1966 01 01 }

(39887.0, 36.5971700, 0.0025920, 39126.0), { 1968 02 01 }

(41317.0, 42.184, 0.0, 0.0), { 1972 01 01 }

(41499.0, 43.184, 0.0, 0.0), { 1972 07 01 }

(41683.0, 44.184, 0.0, 0.0), { 1973 01 01 }

(42048.0, 45.184, 0.0, 0.0), { 1974 01 01 }

(42413.0, 46.184, 0.0, 0.0), { 1975 01 01 }

(42778.0, 47.184, 0.0, 0.0), { 1976 01 01 }

(43144.0, 48.184, 0.0, 0.0), { 1977 01 01 }

(43509.0, 49.184, 0.0, 0.0), { 1978 01 01 }

(43874.0, 50.184, 0.0, 0.0), { 1979 01 01 }

(44239.0, 51.184, 0.0, 0.0), { 1980 01 01 }

(44786.0, 52.184, 0.0, 0.0), { 1981 07 01 }

(45151.0, 53.184, 0.0, 0.0), { 1982 07 01 }

(45516.0, 54.184, 0.0, 0.0), { 1983 07 01 }

(46247.0, 55.184, 0.0, 0.0), { 1985 07 01 }

(47161.0, 56.184, 0.0, 0.0), { 1988 01 01 }

(47892.0, 57.184, 0.0, 0.0), { 1990 01 01 }

(48257.0, 58.184, 0.0, 0.0), { 1991 01 01 }

(48804.0, 59.184, 0.0, 0.0), { 1992 07 01 }

(49169.0, 60.184, 0.0, 0.0), { 1993 07 01 }

(49534.0, 61.184, 0.0, 0.0), { 1994 07 01 }

(50083.0, 62.184, 0.0, 0.0), { 1996 01 01 }

(50630.0, 63.184, 0.0, 0.0), { 1997 07 01 }

(51178.0, 64.184, 0.0, 0.0)); { 1999 01 01 }

Begin

JulTT:=JulUTC;

DeltaTA:=0.0;

If JulTT < 2435473.5 Then Exit ;

If ( JulUTC < ( Tabl[1,1] + DifEpoch ) ) Then

Begin

If ( JulUTC < ( TablOne[1,1] + DifEpoch ) ) Then

Begin

DeltaTA:=TablOne[1,2]*(JulUTC-2415020.5)/(TablOne[1,1]-15020.0);

DeltaTA:=DeltaTA/86400;

End

Else

Begin

For i:=1 To 6 Do

If ( ( TablOne[i,1] <= ( JulUTC - DifEpoch ) )

AND ( ( JulUTC - DifEpoch ) < TablOne[i+1,1] ) )

Then

Begin

r:=(JulUTC-DifEpoch

-TablOne[i,1])/(TablOne[i+1,1]-TablOne[i,1]);

DeltaTA:=TablOne[i,2]+r*(TablOne[i+1,2]-TablOne[i,2]);

DeltaTA:=DeltaTA/86400;

End;

JulTT:=JulUTC+DeltaTA;

Exit;

End;

for i:=1 to MaxDTable-1 do

If ( ( Tabl[i,1] <= ( JulUTC - DifEpoch ) )

AND ( ( JulUTC - DifEpoch ) < Tabl[i+1,1] ) )

then

DeltaTA:=(Tabl[i,2]+Tabl[i,3]

*(JulUTC-DifEpoch-Tabl[i,4]))/86400.e0;

If ( JulUTC - DifEpoch ) > Tabl[MaxDTable,1]

then

DeltaTA:=Tabl[MaxDTable,2]/86400;

JulTT:=JulUTC+DeltaTA;

End;

{ ToGetTEph is a function

to obtain the moment in TE with the use a value for TT }

Function ToGetTEph ( JulTT : Extended ) : Extended ;

Var

d,g : Extended;

Begin

d:=(JulTT-JD2000)/JulianC;

g:=GraRad*(357.528+35999.050*d);

ToGetTEph:=JulTT+(0.001658*Sin(g+0.0167*Sin(g)))/86400.e0;

End;

{ to get Universal Time if UTC and DeltaUT1 are known }

Function ToGetUT1 ( JulUTC , DeltaUT1 : Extended ) : Extended ;

Begin

ToGetUT1:=JulUTC+DeltaUT1;

End;

{ ClcGrMSTime is a function

to obtain Greenwich Mean Sidereal Time in radians .

IAU 1976 for the Greenwich Mean Sidereal Time .

TinUTC in UTC and close to UT1 }

Function ToGetGrMSTime( TinUTC : Extended ) : Extended ;

Var

t,aj,ajj,ajs,s,s0,r,freq : Extended ;

SidTime : Extended ;

IntPart : Integer;

Begin

t:=TinUTC;

aj:=Trunc(t)+0.5;

s:=t-aj;

ajj:=aj-2451545.0;

ajs:=ajj/36525.0;

s0:=1.753368559233266e0+(628.3319706888409e0

+(6.770713944903336e-06-4.508767234318685e-10*ajs)*ajs)*ajs;

freq:=1.002737909350795e0

+(5.900575455674703e-11-5.893984333409384e-15*ajs)*ajs;

s0:=s0+freq*s*2*Pi;

r:=s0/(2*Pi);

IntPart:=Trunc(r);

SidTime:=s0-2*Pi*IntPart;

if SidTime < 0 then SidTime:=SidTime+2*Pi;

ToGetGrMSTime:=SidTime;

End;

END. // unit UnForTim

Параметры нутации

Алгоритм вычисления параметров нутации представлен ниже.

UNIT UnForNut ; { constant 1980 IAU Theory of Nutation ( J. M. Wahr ) }

INTERFACE

Uses

UnConTyp ; // for TFundArg , GraRad , JD2000 , JulianC const

Procedure ClcFundArg ( TCur : Extended;

Var FundArg : TFundArg);

Procedure ClcNut ( TCur : Extended;

Var DeltaPsi,DeltaEps,EpsMean : Extended ) ;

IMPLEMENTATION

{ ClcFundArg is a procedure

to calculate Fundamental Arguments in radians

on moment TCur in julian day.

FundArg[1..5] is massive of the values of arguments.

global const - GraRad - for change from degrees to radians. }

Procedure ClcFundArg ( TCur : Extended ;

Var FundArg : TFundArg ) ;

Var

r,dt,dt2,dt3 : Extended;

i,l : Integer;

Begin

dt:=(TCur-JD2000)/JulianC; // const from UnConTyp

dt2:=sqr(dt); dt3:=dt*dt2;

FundArg[1]:=218.31643250e0+481267.8812772222e0*dt

-0.00161167e0*dt2+0.00000528e0*dt3;

FundArg[2]:=134.96298139e0+477198.8673980556e0*dt

+0.00869722e0*dt2+0.00001778e0*dt3;

FundArg[3]:=357.52772333e0+35999.05034e0*dt

-0.00016028e0*dt2-0.00000333e0*dt3;

FundArg[4]:= 93.27191028e0+483202.0175380555e0*dt

-0.00368250e0*dt2+0.00000306e0*dt3;

FundArg[5]:=297.85036306e0+445267.11148e0*dt

-0.00191417e0*dt2+0.00000528e0*dt3;

for i:=1 to 5 do

Begin

r:=FundArg[i]; l:=Trunc(r/360.e0); FundArg[i]:=GraRad*(r-360.e0*l);

End;

{ ClcNut is a procedure

to calculate the parameters of nutation in longitude and obliquity

and value of the mean obliquity in radians on the moment TCur.

the arguments and coefficients

of nutation theory are contained in the NutArg and CoefNut .

global parameter - SecRad - for change from arcsecs to radians . }

Procedure ClcNut ( TCur : Extended ;

Var DeltaPsi,DeltaEps,EpsMean : Extended ) ;

Var

i,k : Integer ;

r,dt,dt2,dt3 : Extended ;

FundArg : TFundArg ;

Const

{ 1980 IAU Theory of Nutation ( J. M. Wahr ) }

NutArg : Array [1..106,1..5] of ShortInt = { arguments }

( ( 1 , 0 , 0 , 0 , 0 ) ,

( 2 , 0 , 0 , 0 , 0 ) ,

( 1 , -2 , 0 , 2 , 0 ) ,

( 0 , 2 , 0 , -2 , 0 ) ,

( 2 , -2 , 0 , 2 , 0 ) ,

( 0 , 1 , -1 , 0 , -1 ) ,

( 1 , 0 , -2 , 2 , -2 ) ,

( 1 , 2 , 0 , -2 , 0 ) ,

( 2 , 0 , 0 , 2 , -2 ) ,

( 0 , 0 , 1 , 0 , 0 ) ,

( 2 , 0 , 1 , 2 , -2 ) ,

( 2 , 0 , -1 , 2 , -2 ) ,

( 1 , 0 , 0 , 2 , -2 ) ,

( 0 , 2 , 0 , 0 , -2 ) ,

( 0 , 0 , 0 , 2 , -2 ) ,

( 0 , 0 , 2 , 0 , 0 ) ,

( 1 , 0 , 1 , 0 , 0 ) ,

( 2 , 0 , 2 , 2 , -2 ) ,

( 1 , 0 , -1 , 0 , 0 ) ,

( 1 , -2 , 0 , 0 , 2 ) ,

( 1 , 0 , -1 , 2 , -2 ) ,

( 1 , 2 , 0 , 0 , -2 ) ,

( 1 , 0 , 1 , 2 , -2 ) ,

( 0 , 1 , 0 , 0 , -1 ) ,

( 0 , 2 , 1 , 0 , -2 ) ,

( 1 , 0 , 0 , -2 , 2 ) ,

( 0 , 0 , 1 , -2 , 2 ) ,

( 2 , 0 , 1 , 0 , 0 ) ,

( 1 , -1 , 0 , 0 , 1 ) ,

( 0 , 0 , 1 , 2 , -2 ) ,

( 2 , 0 , 0 , 2 , 0 ) ,

( 0 , 1 , 0 , 0 , 0 ) ,

( 1 , 0 , 0 , 2 , 0 ) ,

( 2 , 1 , 0 , 2 , 0 ) ,

( 0 , 1 , 0 , 0 , -2 ) ,

( 2 , -1 , 0 , 2 , 0 ) ,

( 0 , 0 , 0 , 0 , 2 ) ,

( 1 , 1 , 0 , 0 , 0 ) ,

( 1 , -1 , 0 , 0 , 0 ) ,

( 2 , -1 , 0 , 2 , 2 ) ,

( 1 , 1 , 0 , 2 , 0 ) ,

( 2 , 0 , 0 , 2 , 2 ) ,

( 0 , 2 , 0 , 0 , 0 ) ,

( 2 , 1 , 0 , 2 , -2 ) ,

( 2 , 2 , 0 , 2 , 0 ) ,

( 0 , 0 , 0 , 2 , 0 ) ,

( 1 , -1 , 0 , 2 , 0 ) ,

( 1 , -1 , 0 , 0 , 2 ) ,

( 1 , 1 , 0 , 0 , -2 ) ,

( 1 , -1 , 0 , 2 , 2 ) ,

( 0 , 1 , 1 , 0 , -2 ) ,

( 2 , 0 , 1 , 2 , 0 ) ,

( 2 , 0 , -1 , 2 , 0 ) ,

( 2 , 1 , 0 , 2 , 2 ) ,

( 0 , 1 , 0 , 0 , 2 ) ,

( 2 , 2 , 0 , 2 , -2 ) ,

( 1 , 0 , 0 , 0 , 2 ) ,

( 1 , 0 , 0 , 2 , 2 ) ,

( 1 , 1 , 0 , 2 , -2 ) ,

( 1 , 0 , 0 , 0 , -2 ) ,

( 0 , 1 , -1 , 0 , 0 ) ,

( 1 , 2 , 0 , 2 , 0 ) ,

( 0 , 0 , 1 , 0 , -2 ) ,

( 0 , 1 , 0 , -2 , 0 ) ,

( 0 , 0 , 0 , 0 , 1 ) ,

( 0 , 1 , 1 , 0 , 0 ) ,

( 0 , 1 , 0 , 2 , 0 ) ,

( 2 , 1 , -1 , 2 , 0 ) ,

( 2 , -1 , -1 , 2 , 2 ) ,

( 1 , -2 , 0 , 0 , 0 ) ,

( 2 , 3 , 0 , 2 , 0 ) ,

( 2 , 0 , -1 , 2 , 2 ) ,

( 2 , 1 , 1 , 2 , 0 ) ,

( 1 , -1 , 0 , 2 , -2 ) ,

( 1 , 2 , 0 , 0 , 0 ) ,

( 2 , 1 , 0 , 0 , 0 ) ,

( 0 , 3 , 0 , 0 , 0 ) ,

( 2 , 0 , 0 , 2 , 1 ) ,

( 2 , -1 , 0 , 0 , 0 ) ,

( 0 , 1 , 0 , 0 , -4 ) ,

( 2 , -2 , 0 , 2 , 2 ) ,

( 2 , -1 , 0 , 2 , 4 ) ,

( 0 , 2 , 0 , 0 , -4 ) ,

( 2 , 1 , 1 , 2 , -2 ) ,

( 1 , 1 , 0 , 2 , 2 ) ,

( 2 , -2 , 0 , 2 , 4 ) ,

( 2 , -1 , 0 , 4 , 0 ) ,

( 0 , 1 , -1 , 0 , -2 ) ,

( 1 , 2 , 0 , 2 , -2 ) ,

( 2 , 2 , 0 , 2 , 2 ) ,

( 1 , 1 , 0 , 0 , 2 ) ,

( 2 , 0 , 0 , 4 , -2 ) ,

( 2 , 3 , 0 , 2 , -2 ) ,

( 0 , 1 , 0 , 2 , -2 ) ,

( 1 , 0 , 1 , 2 , 0 ) ,

( 1 , -1 , -1 , 0 , 2 ) ,

( 1 , 0 , 0 , -2 , 0 ) ,

( 2 , 0 , 0 , 2 , -1 ) ,

( 0 , 0 , 1 , 0 , 2 ) ,

( 0 , 1 , 0 , -2 , -2 ) ,

( 1 , 0 , -1 , 2 , 0 ) ,

( 1 , 1 , 1 , 0 , -2 ) ,

( 0 , 1 , 0 , -2 , 2 ) ,

( 0 , 2 , 0 , 0 , 2 ) ,

( 2 , 0 , 0 , 2 , 4 ) ,

( 0 , 0 , 1 , 0 , 1 ) ) ;

{ 1980 IAU Theory of Nutation ( J. M. Wahr ) }

CoefNut : Array [1..106,1..4] of Single = { amplitudes }

( ( -17.199600 , -0.017420 , 9.202500 , 0.000890 ) ,

( 0.206200 , 0.000020 , -0.089500 , 0.000050 ) ,

( 0.004600 , 0.000000 , -0.002400 , 0.000000 ) ,

( 0.001100 , 0.000000 , 0.000000 , 0.000000 ) ,

( -0.000300 , 0.000000 , 0.000100 , 0.000000 ) ,

( -0.000300 , 0.000000 , 0.000000 , 0.000000 ) ,

( -0.000200 , 0.000000 , 0.000100 , 0.000000 ) ,

( 0.000100 , 0.000000 , 0.000000 , 0.000000 ) ,

( -1.318700 , -0.000160 , 0.573600 , -0.000310 ) ,

( 0.142600 , -0.000340 , 0.005400 , -0.000010 ) ,

( -0.051700 , 0.000120 , 0.022400 , -0.000060 ) ,

( 0.021700 , -0.000050 , -0.009500 , 0.000030 ) ,

( 0.012900 , 0.000010 , -0.007000 , 0.000000 ) ,

( 0.004800 , 0.000000 , 0.000100 , 0.000000 ) ,

( -0.002200 , 0.000000 , 0.000000 , 0.000000 ) ,

( 0.001700 , -0.000010 , 0.000000 , 0.000000 ) ,

( -0.001500 , 0.000000 , 0.000900 , 0.000000 ) ,

( -0.001600 , 0.000010 , 0.000700 , 0.000000 ) ,

( -0.001200 , 0.000000 , 0.000600 , 0.000000 ) ,

( -0.000600 , 0.000000 , 0.000300 , 0.000000 ) ,

( -0.000500 , 0.000000 , 0.000300 , 0.000000 ) ,

( 0.000400 , 0.000000 , -0.000200 , 0.000000 ) ,

( -0.000400 , 0.000000 , 0.000000 , 0.000000 ) ,

( 0.000100 , 0.000000 , 0.000000 , 0.000000 ) ,

( -0.000100 , 0.000000 , 0.000000 , 0.000000 ) ,

( 0.000100 , 0.000000 , 0.000000 , 0.000000 ) ,

( -0.000100 , 0.000000 , 0.000000 , 0.000000 ) ,

( -0.227400 , -0.000020 , 0.097700 , -0.000050 ) ,

( 0.071200 , 0.000010 , -0.000700 , 0.000000 ) ,

( -0.038600 , -0.000040 , 0.020000 , 0.000000 ) ,

( -0.030100 , 0.000000 , 0.012900 , -0.000010 ) ,

( -0.015800 , 0.000000 , -0.000100 , 0.000000 ) ,

( 0.012300 , 0.000000 , -0.005300 , 0.000000 ) ,

( 0.006300 , 0.000000 , -0.000200 , 0.000000 ) ,

( 0.006300 , 0.000010 , -0.003300 , 0.000000 ) ,

( -0.005800 , -0.000010 , 0.003200 , 0.000000 ) ,

( -0.005900 , 0.000000 , 0.002600 , 0.000000 ) ,

( -0.005100 , 0.000000 , 0.002700 , 0.000000 ) ,

( -0.003800 , 0.000000 , 0.001600 , 0.000000 ) ,

( 0.002900 , 0.000000 , -0.000100 , 0.000000 ) ,

( 0.002900 , 0.000000 , -0.001200 , 0.000000 ) ,

( -0.003100 , 0.000000 , 0.001300 , 0.000000 ) ,

( 0.002600 , 0.000000 , -0.000100 , 0.000000 ) ,

( 0.002100 , 0.000000 , -0.001000 , 0.000000 ) ,

( 0.001600 , 0.000000 , -0.000800 , 0.000000 ) ,

( -0.001300 , 0.000000 , 0.000700 , 0.000000 ) ,

( -0.001000 , 0.000000 , 0.000500 , 0.000000 ) ,

( -0.000700 , 0.000000 , 0.000000 , 0.000000 ) ,

( 0.000700 , 0.000000 , -0.000300 , 0.000000 ) ,

( -0.000700 , 0.000000 , 0.000300 , 0.000000 ) ,

( -0.000800 , 0.000000 , 0.000300 , 0.000000 ) ,

( 0.000600 , 0.000000 , 0.000000 , 0.000000 ) ,

( 0.000600 , 0.000000 , -0.000300 , 0.000000 ) ,

( -0.000600 , 0.000000 , 0.000300 , 0.000000 ) ,

( -0.000700 , 0.000000 , 0.000300 , 0.000000 ) ,

( 0.000600 , 0.000000 , -0.000300 , 0.000000 ) ,

( -0.000500 , 0.000000 , 0.000300 , 0.000000 ) ,

( 0.000500 , 0.000000 , 0.000000 , 0.000000 ) ,

( -0.000500 , 0.000000 , 0.000300 , 0.000000 ) ,

( -0.000400 , 0.000000 , 0.000000 , 0.000000 ) ,

( 0.000400 , 0.000000 , 0.000000 , 0.000000 ) ,

( -0.000400 , 0.000000 , 0.000000 , 0.000000 ) ,

( -0.000300 , 0.000000 , 0.000000 , 0.000000 ) ,

( 0.000300 , 0.000000 , 0.000000 , 0.000000 ) ,

( -0.000300 , 0.000000 , 0.000100 , 0.000000 ) ,

( -0.000200 , 0.000000 , 0.000100 , 0.000000 ) ,

( -0.000300 , 0.000000 , 0.000100 , 0.000000 ) ,

( 0.000200 , 0.000000 , -0.000100 , 0.000000 ) ,

( -0.000200 , 0.000000 , 0.000100 , 0.000000 ) ,

( 0.000200 , 0.000000 , -0.000100 , 0.000000 ) ,

( -0.000200 , 0.000000 , 0.000100 , 0.000000 ) ,

( 0.000200 , 0.000000 , 0.000000 , 0.000000 ) ,

( 0.000200 , 0.000000 , -0.000100 , 0.000000 ) ,

( 0.000100 , 0.000000 , -0.000100 , 0.000000 ) ,

( -0.000100 , 0.000000 , 0.000000 , 0.000000 ) ,

( 0.000100 , 0.000000 , -0.000100 , 0.000000 ) ,

( -0.000200 , 0.000000 , 0.000100 , 0.000000 ) ,

( -0.000100 , 0.000000 , 0.000000 , 0.000000 ) ,

( 0.000100 , 0.000000 , -0.000100 , 0.000000 ) ,

( -0.000100 , 0.000000 , 0.000100 , 0.000000 ) ,

( 0.000100 , 0.000000 , 0.000000 , 0.000000 ) ,

( 0.000100 , 0.000000 , -0.000100 , 0.000000 ) ,

( -0.000100 , 0.000000 , 0.000000 , 0.000000 ) ,

( 0.000100 , 0.000000 , 0.000000 , 0.000000 ) ,

( -0.000100 , 0.000000 , 0.000000 , 0.000000 ) ,

( 0.000100 , 0.000000 , 0.000000 , 0.000000 ) ,

( -0.000100 , 0.000000 , 0.000000 , 0.000000 ) ,

( 0.000100 , 0.000000 , 0.000000 , 0.000000 ) ,

( -0.000100 , 0.000000 , 0.000000 , 0.000000 ) ,

( 0.000100 , 0.000000 , 0.000000 , 0.000000 ) ) ;

Begin

dt:=(TCur-JD2000)/JulianC; // from UnConTyp

dt2:=dt*dt;

dt3:=dt*dt2;

EpsMean:=SecRad*(84381.448-46.8150*dt-0.00059*dt2+0.001813*dt3);

ClcFundArg(TCur,FundArg);

DeltaPsi:=0.e0;

DeltaEps:=0.e0;

for k:=1 to 106 do

Begin

r:=NutArg[k,1]*(FundArg[1]-FundArg[4]);

for i:=2 to 5 do

r:=r+NutArg[k,i]*FundArg[i];

DeltaPsi:=DeltaPsi+(CoefNut[k,1]

+dt*CoefNut[k,2])*Sin(r);

DeltaEps:=DeltaEps+(CoefNut[k,3]

+dt*CoefNut[k,4])*Cos(r);

End;

DeltaPsi:=SecRad*DeltaPsi; { in radian }

DeltaEps:=SecRad*DeltaEps;

End;

END. // unit UnForNut

Вычисление рефракции

Для вычисления поправки за рефракцию доктор Константин Владиславович Куимов предлагает следующий алгоритм.

UNIT UnRefrOl ; // procedure to calculate the refraction correction

INTERFACE // has been presented by doctor Konstantin V. Kuimov

Var

wlmkm : Extended ; { length of wave }

temprc : Extended ; { temperature in degree Celsius }

prepas : Extended ; { pressure in pascal 760 mm = 101325 Pa }

humrun : Extended ; { humidity in unit interval [0,1] }

Procedure SimpleDataForRefraction ;

Procedure ToGetDataForRefraction ;

Procedure RefrOlik ( wlmkm : Extended ;

temp , ppas , hum : Extended ;

z : Extended ;

Var ref : Extended ) ;

Function ToGetRefraction ( alt : Extended ) : Extended ;

IMPLEMENTATION

Uses

UnConTyp ; { for GraRad , RadGra and RadSec }

Procedure SimpleDataForRefraction ;

Begin

{ initial value }

wlmkm:=0.5931; { lehgtn of wave in mkm }

temprc:=0.0; { in degree of Celsius }

prepas:=99325.158; { in pascal 745 in mm }

humrun:=0.9; { relative humidity is 90% }

End;

Procedure ToGetDataForRefraction ;

Var

ft : Text ;

tempk,presh,hum : Extended ;

Begin

{ initial value }

wlmkm:=0.5931; { in mkm }

tempk:=273.15; { in Kelvin }

presh:=745.0; { in mm }

hum:=90.0; { in % }

AssignFile(ft,'RefrParm.Dat');

{$I-}

ReSet(ft);

{$I+}

If IOResult = 0

Then

Begin

If NOT EOF(ft) Then ReadLn(ft,tempk);

If NOT EOF(ft) Then ReadLn(ft,presh);

If NOT EOF(ft) Then ReadLn(ft,hum);

CloseFile(ft);

End;

temprc:=tempk-273.15; { to Celsius }

prepas:=101325.0*(presh/760.0); { pressure in pascal }

humrun:=hum/100.0; { to unit interval }

End;

Procedure RefrOlik ( wlmkm : Extended ;

temp , ppas , hum : Extended ;

z : Extended ;

Var ref : Extended ) ;

Const

R = 8314.72e0; { газовая постоянная, дж/(К * кмоль) }

rw = 8314.72e0/18.0152e0; { газовая постоянная/мол.вес воды }

g0 = 9.80665e0; { ускорение силы тяжести, м/с**2 }

Dm = 28.9644e0; { молярная масса сухого воздуха, кг/моль }

r0 = 6371000.e0; { средний радиус Земли, м }

Hpi = 1.570796326794897e0; { pi/2 }

sqpi2 = 0.8862269254527580e0; { sqrt(pi)/2 }

Var

tabs,tabs2,tabs3,t,roln,pw,ps,

ds,dw,s2,s4,s6,alpha0,b0,a,x,y : Extended ;

Begin

Tabs := temp + 273.15e0; { абсолютная температура }

t := temp; { температура по Цельсию }

roln := -5.32917e0 + t * (0.0688825e0 + t * (-2.9815e-4

+ 1.39e-6 * t)); { нат. логарифм плотности водяного пара (см. [2]) }

pw := Hum * rw * Tabs * exp(roln); { давление водяного пара }

Ps := (Ppas - pw)/100; { выразить в миллибарах давление сухого воздуха }

Pw := Pw/100; { выразить в миллибарах давление водяного пара }

tabs2 := Sqr(Tabs);

tabs3 := tabs2 * Tabs;

Ds := Ps/Tabs

* (1 + Ps * (57.90e-8 - 9.3250e-4/Tabs + 0.25844e0/tabs2));

Dw := Pw/Tabs * (1 + Pw * (1 + 3.7e-4 * Pw)

* (-2.37321e-3 + 2.23366e0/Tabs - 710.792e0/tabs2

+ 7.75141e4/tabs3));

s2 := 1/Sqr(Wlmkm);

s4 := s2*s2;

s6 := s2 * s4;

alpha0 :=((2371.34e0 + 683939.7e0/(130-s2) + 4547.3e0/(38.9e0-s2)) * Ds

+(6487.31e0+58.058e0*s2-0.71150e0*s4+0.08851e0*s6) * Dw)/1.e8;

b0 := R * Tabs/(g0 * Dm * r0); { beta0 }

a := b0 - alpha0/2;

x := cos(z)/sqrt(2 * a);

t := 1.0e0 - 7.5e0/(abs(x) + 3.75e0); { NAG13 s15adf.for }

y := (((((((((((((((+3.328130055126039e-10

* t - 5.718639670776992e-10) * t - 4.066088879757269e-9)

* t + 7.532536116142436e-9) * t + 3.026547320064576e-8)

* t - 7.043998994397452e-8) * t - 1.822565715362025e-7)

* t + 6.575825478226343e-7) * t + 7.478317101785790e-7)

* t - 6.182369348098529e-6) * t + 3.584014089915968e-6)

* t + 4.789838226695987e-5) * t - 1.524627476123466e-4)

* t - 2.553523453642242e-5) * t + 1.802962431316418e-3)

* t - 8.220621168415435e-3) * t + 2.414322397093253e-2;

y := (((((y * t - 5.480232669380236e-2) * t + 1.026043120322792e-1)

* t - 1.635718955239687e-1) * t + 2.260080669166197e-1)

* t - 2.734219314954260e-1) * t + 1.455897212750385e-1;

ref := alpha0 * (1 - alpha0/2) * sqrt(2/a) * sin(z) * y * sqpi2;

End;

Function ToGetRefraction ( alt : Extended ) : Extended ;

Var

z,zr,ref : Extended ;

Begin

ref:=0;

z:=90-alt;

If z > 90.0 Then z:=90.0; { if altitude < 0 degree }

If z < 0.00 Then z:=0.00;

zr:=GraRad*z; { to radian }

RefrOlik(wlmkm,temprc,prepas,humrun,zr,ref); { ref in radian }

ToGetRefraction:=RadGra*ref; { in degree }

End;

END.

Преобразование координат

Алгоритмы, представленные ниже, позволяют вычислять матрицы прецессии, нутации, матрицы перехода между системами стандартного экватора и истинного экватора даты. Здесь же вычисляется матрица преобразования от истинного экватора в топоцентрическую систему отсчёта, связанную с пунктом наблюдений. Две процедуры осуществляют переход между геодезическими сферическими и геоцентрическими прямоугольными координатами обсерватории.

UNIT UnForCoo ; // simple procedure to change coordinate systems

INTERFACE

Uses

UnConTyp ; // for type TMatr33 TVect3 const JD2000 JulianC

Procedure ClcPrecMatr ( Epoch , TCur : Extended ;

Var PrecMatr : TMatr33 );

Procedure ClcThreeRotMatr ( Epoch , TCur : Extended ;

PlaceCoor : TPlaceCooRec ;

Var

RotMatr ,

RosMatr ,

TopMatr : TMatr33 );

Procedure ClcSpherCoor( Pos : TVect3 ;

Var Alfa , Delta, Ro : Extended );

Procedure ClcSpherCoorInDegree( Pos : TVect3 ;

Var Alfa , Delta, Ro : Extended );

Procedure ClcTopCoorInDegree ( Pos : TVect3 ;

Var Azt , Alt, Ro : Extended );

Procedure GeoSpherFromCart ( x , y , z : Extended ;

Var f , v , h : Extended ) ;

Procedure GeoCartFromSpher ( f , v , h : Extended ;

Var x , y , z : Extended ) ;

Procedure DescFromSpherCoor ( Az , Alt , Ro : Extended ;

Var PosDescart : TVect3 ) ;

IMPLEMENTATION

Uses

UnForFun , // procedure ArcTanTwoArg for example

UnForTim , // function ToGetGrMSTime

UnForNut , // for nutation procedure

UnRefrOl ; // for ToGetRefraction function

{ ClcPrecMatr is a procedure

to calculate the Matrix of precession

for change from Epoch to TCur or another words

from the fixed equator and ecliptic to the instantaneous ones.

Fundamental EPOCH is J2000.0

but Epoch as parameter

may be any fixed epoch in JD

and TCur is epoch of the Date in JD .

PrecMatr[1..3,1..3] is the resulting massive .

global const - SecRad - for change from arcsec to radians . }

Procedure ClcPrecMatr ( Epoch , TCur : Extended ;

Var PrecMatr : TMatr33 );

Var

Dzita0,Teta,Zet : Extended;

de,de2,dt,dt2,dt3 : Extended;

r,u1,u2,u3,s1,s2,s3,c1,c2,c3 : Extended;

Begin

de:=(Epoch-JD2000)/JulianC; de2:= de*de; // constants from UnConTyp

dt:=(TCur-Epoch)/JulianC; dt2:= dt*dt; dt3:=dt2*dt;

r:=(2306.2181+1.39656*de-0.000139*de2)*dt;

Dzita0:=r+(0.30188-0.000344*de)*dt2+0.017998*dt3;

Teta:=(2004.3109-0.85330*de-0.000217*de2)*dt

-(0.42665+0.000217*de)*dt2-0.041833*dt3;

Zet:=r+(1.09468+0.000066*de)*dt2+0.018203*dt3;

u1:=SecRad*Dzita0; u2:=SecRad*Teta; u3:=SecRad*Zet;

s1:=Sin(u1); s2:=Sin(u2); s3:=Sin(u3);

c1:=Cos(u1); c2:=Cos(u2); c3:=Cos(u3);

PrecMatr[1,1]:=+c3*c2*c1-s3*s1;

PrecMatr[1,2]:=-c3*c2*s1-s3*c1; PrecMatr[1,3]:=-c3*s2;

PrecMatr[2,1]:=+s3*c2*c1+c3*s1;

PrecMatr[2,2]:=-s3*c2*s1+c3*c1; PrecMatr[2,3]:=-s3*s2;

PrecMatr[3,1]:=+s2*c1; PrecMatr[3,2]:=-s2*s1; PrecMatr[3,3]:=+c2;

End;

{ ClcThreeRotMatr is a procedure

to calculate the Matrix to change

from the fixed equator and ecliptic

to the system of the true equiator and instantaneous ecliptic.

Fundamental EPOCH is J2000.0

TCur is epoch of the Date in JD .

PrecMatr and MatrNut is declared in UnConTyp

RotMatr[1..3,1..3] is the resulting massive .

RosMatr[1..3,1..3] is the matrix for the sidereal time

TopMatr[1..3,1..3] is the matrix for change to topocentric system }

Procedure ClcThreeRotMatr ( Epoch , TCur : Extended ;

PlaceCoor : TPlaceCooRec ;

Var

RotMatr ,

RosMatr ,

TopMatr : TMatr33 );

Var

i,j,k : Byte ;

uf,ul,sf,cf,

sl,cl,st,ct,

e,s,c,se,ce,ts,SidTime,

DeltaPsi,DeltaEps,EpsMean : Extended ;

MatrCur,MatrEps,MatrPsi : TMatr33 ;

PosGeodetic : TVect3 ;

Begin

ClcNut(TCur,DeltaPsi,DeltaEps,EpsMean); { in radian from UnForNut}

se:=Sin(EpsMean);

ce:=Cos(EpsMean);

MatrEps[1,1]:=+1;

MatrEps[1,2]:=0;

MatrEps[1,3]:=0;

MatrEps[2,1]:=0;

MatrEps[2,2]:=+ce;

MatrEps[2,3]:=+se;

MatrEps[3,1]:=0;

MatrEps[3,2]:=-se;

MatrEps[3,3]:=+ce;

For i:=1 To 3 Do

For j:=1 To 3 Do

MatrEcl[i,j]:=MatrEps[i,j];

s:=Sin(DeltaPsi);

c:=Cos(DeltaPsi);

MatrPsi[1,1]:=+c;

MatrPsi[1,2]:=-s;

MatrPsi[1,3]:=0;

MatrPsi[2,1]:=+s;

MatrPsi[2,2]:=+c;

MatrPsi[2,3]:=0;

MatrPsi[3,1]:=0;

MatrPsi[3,2]:=0;

MatrPsi[3,3]:=+1;

For i:=1 To 3 Do

For j:=1 To 3 Do

Begin

MatrCur[i,j]:=0;

For k:=1 To 3 Do

MatrCur[i,j]:=MatrCur[i,j]+MatrPsi[i,k]*MatrEps[k,j];

End;

e:=EpsMean+DeltaEps;

s:=Sin(e);

c:=Cos(e);

MatrEps[1,1]:=+1;

MatrEps[1,2]:=0;

MatrEps[1,3]:=0;

MatrEps[2,1]:=0;

MatrEps[2,2]:=+c;

MatrEps[2,3]:=-s;

MatrEps[3,1]:=0;

MatrEps[3,2]:=+s;

MatrEps[3,3]:=+c;

For i:=1 To 3 Do

For j:=1 To 3 Do

Begin

MatrNut[i,j]:=0;

For k:=1 To 3 Do

MatrNut[i,j]:=MatrNut[i,j]+MatrEps[i,k]*MatrCur[k,j];

End;

ClcPrecMatr(Epoch,TCur,PrecMatr);

For i:=1 To 3 Do

For j:=1 To 3 Do

Begin

RotMatr[i,j]:=0;

For k:=1 To 3 Do

RotMatr[i,j]:=RotMatr[i,j]+MatrNut[i,k]*PrecMatr[k,j];

End;

SidTime:=ToGetGrMSTime(TCur); { mean sidereal time unit UnForTim }

ts:=SidTime+DeltaPsi*ce; { true sidereal time }

s:=Sin(ts);

c:=Cos(ts);

RosMatr[1,1]:=+c;

RosMatr[1,2]:=+s;

RosMatr[1,3]:=0;

RosMatr[2,1]:=-s;

RosMatr[2,2]:=+c;

RosMatr[2,3]:=0;

RosMatr[3,1]:=0;

RosMatr[3,2]:=0;

RosMatr[3,3]:=+1;

{ TopMatr[...] to change from true equator to horizont }

With PlaceCoor Do

Begin

uf:=GraRad*f;

ul:=GraRad*l;

PosGeodetic[1]:=1.0e-3*x; { in km }

PosGeodetic[2]:=1.0e-3*y;

PosGeodetic[3]:=1.0e-3*z;

End;

{ station position from greenwich to true equator }

For i:=1 To 3 Do

Begin

StationPos[i]:=0; // global variable UnConTyp

For j:=1 To 3 Do

StationPos[i]:=StationPos[i]+RosMatr[j,i]*PosGeodetic[j];

End;

sf:=Sin(uf);

cf:=Cos(uf);

sl:=Sin(ul);

cl:=Cos(ul);

st:=s*cl+c*sl; { sin ( teta + lambda ) }

ct:=c*cl-s*sl; { cos ( teta + lambda ) }

ArcTanTwoArg(st,ct,LocalSTime); { local sidereal time in radians }

LocalSTime

TopMatr[1,1]:=+sf*ct;

TopMatr[1,2]:=+sf*st;

TopMatr[1,3]:=-cf;

TopMatr[2,1]:=-st;

TopMatr[2,2]:=+ct;

TopMatr[2,3]:=0;

TopMatr[3,1]:=+cf*ct;

TopMatr[3,2]:=+cf*st;

TopMatr[3,3]:=+sf;

End;

{ ClcSpherCoor is a procedure

to obtain the spherical coordinates of the planet from cartesian ones }

Procedure ClcSpherCoor( Pos : TVect3 ;

Var Alfa , Delta, Ro : Extended );

Var

xr : TVect3;

r,rr,q,qq,p,f,s,c,a,b : Extended;

i : Integer;

Begin

for i:=1 to 3 do xr[i]:=Pos[i];

qq:=Sqr(xr[1])+Sqr(xr[2]); rr:=qq+Sqr(xr[3]); r:=Sqrt(rr);

q:=Sqrt(qq); p:=1/q; f:=1/r;

c:=p*xr[1]; s:=p*xr[2];

ArcTanTwoArg(s,c,a); // unit UnForFun

c:=q*f;

s:=xr[3]*f;

ArcTanTwoArg(s,c,b);

if a < 0 then a:=a+Pi22;

Alfa:=a;

Delta:=b;

Ro:=r;

End;

{ ClcSpherCoorInDegree is a procedure

to obtain the spherical coordinates of the planet from cartesian ones }

Procedure ClcSpherCoorInDegree ( Pos : TVect3 ;

Var Alfa , Delta, Ro : Extended );

Var

xr : TVect3;

r,rr,q,qq,p,f,s,c,a,b : Extended;

i : Integer;

Begin

for i:=1 to 3 do xr[i]:=Pos[i];

qq:=Sqr(xr[1])+Sqr(xr[2]); rr:=qq+Sqr(xr[3]); r:=Sqrt(rr);

q:=Sqrt(qq); p:=1/q; f:=1/r;

c:=p*xr[1]; s:=p*xr[2]; ArcTanTwoArg(s,c,a);

c:=q*f; s:=xr[3]*f; ArcTanTwoArg(s,c,b);

if a < 0 then a:=a+Pi22;

Alfa:=RadGra*a;

Delta:=RadGra*b;

Ro:=r;

End;

{ ClcTopCoorInDegree is a procedure

to obtain the topocentric coordinates in spherical form }

Procedure ClcTopCoorInDegree ( Pos : TVect3 ;

Var Azt , Alt, Ro : Extended );

Var

xr : TVect3;

r,rr,q,qq,p,f,s,c,a,b : Extended;

i : Integer;

zr : Extended ;

Begin

for i:=1 to 3 do xr[i]:=Pos[i];

qq:=Sqr(xr[1])+Sqr(xr[2]); rr:=qq+Sqr(xr[3]); r:=Sqrt(rr);

q:=Sqrt(qq); p:=1/q; f:=1/r;

c:=-p*xr[1]; s:=p*xr[2]; ArcTanTwoArg(s,c,a);

c:=q*f; s:=xr[3]*f; ArcTanTwoArg(s,c,b);

if a < 0 then a:=a+Pi22;

Azt:=RadGra*a;

Alt:=RadGra*b;

Ro:=r;

{ to add refraction }

zr:=ToGetRefraction(Alt);

Alt:=Alt+zr;

End;

{ to obtain spherical coordinates ( h in metr )

of the station from the cartesian ones in metr }

Procedure GeoSpherFromCart ( x , y , z : Extended ;

Var f , v , h : Extended ) ;

Var

irf : Byte ;

a,e,r,ff,vv,

pp,hg,gg,sv,cv,

sf,cf,sff,cff,cc : Extended ;

Begin

r:=6378144.11; { in metr }

a:=1.0/298.257;

e:=2*a-a*a;

pp:=Sqrt(x*x+y*y);

If pp < 1.0e4

Then

r:=1.0e-3*r; { in km }

cv:=x/pp;

sv:=y/pp;

hg:=0.0;

For irf:=1 To 4 Do

Begin

sff:=z*(1+hg);

cff:=pp*(1-e+hg);

cc:=Sqrt(sff*sff+cff*cff);

sf:=sff/cc;

cf:=cff/cc;

gg:=r/Sqrt(1-e*sf*sf);

hg:=pp/(gg*cf)-1.0;

End;

ArcTanTwoArg(sf,cf,ff);

f:=RadGra*ff;

ArcTanTwoArg(sv,cv,vv);

v:=RadGra*vv;

If v < 0.0

Then

v:=360+v;

h:=hg*gg;

End;

{ to obtain cartesian coordinates in metr

of the station from the spherical ones }

Procedure GeoCartFromSpher ( f , v , h : Extended ;

Var x , y , z : Extended ) ;

Var

a,e,r,q,

ff,vv,sf,cf,gg : Extended ;

Begin

r:=6378144.11;

a:=1.0/298.257;

e:=2*a-a*a;

q:=1-e;

ff:=GraRad*f;

vv:=GraRad*v;

sf:=Sin(ff);

cf:=Cos(ff);

gg:=r/Sqrt(1-e*sf*sf);

x:=(gg+h)*cf*Cos(vv);

y:=(gg+h)*cf*Sin(vv);

z:=(gg*q+h)*sf;

End;

Procedure DescFromSpherCoor ( Az , Alt , Ro : Extended ;

Var PosDescart : TVect3 ) ;

Var

a,b,cb : Extended ;

Begin

a:=GraRad*Az;

b:=GraRad*Alt;

cb:=Cos(b);

PosDescart[1]:=Ro*Cos(a)*cb;

PosDescart[2]:=Ro*Sin(a)*cb;

PosDescart[3]:=Ro*Sin(b);

End;

END. // unit UnForCoo

Двустрочные элементы

Начальные параметры движения искусственных спутников Земли в виде средних кеплеровских элементов орбиты на заданную дату можно найти в Интернете по адресу http://celestrak.com/NORAD/elements/ . Эти данные ежедневно публикует радиотехническая служба слежения за спутниками США. Формат параметров называют двустрочными орбитальными элементами НОРАД.

Данные для каждого спутника содержат три строки:

AAAAAAAAAAA

1 NNNNNU NNNNNAAA NNNNN.NNNNNNNN +.NNNNNNNN +NNNNN-N +NNNNN-N N NNNNN

2 NNNNN NNN.NNNN NNN.NNNN NNNNNNN NNN.NNNN NNN.NNNN NN.NNNNNNNNNNNNNN

В строке 0 записано имя объекта из одиннадцати символов.

Строка 1

Колонка Описание

01-01 номер линии,

03-07 номер спутника,

10-11 международный указатель (последние две цифры года запуска),

12-14 международный указатель (порядковый номер запуска),

15-17 международный указатель (часть запущенного изделия),

19-20 эпоха элементов (последние две цифры года),

21-32 эпоха элементов (день с дробной частью от начала года),

34-43 первая производная от среднего движения,

45-52 вторая производная от среднего движения,

54-61 эффективный коэффициент отражения,

63-63 тип эфемериды,

65-68 номер элементов,

69-69 контрольная сумма.

Строка 2

колонка Описание

01-01 номер линии,

03-07 номер спутника,

09-16 угол наклонения (градусы),

18-25 прямое восхождение восходящего узла (градусы),

27-33 эксцентриситет (впереди числа ставить десятичную точку),

35-42 аргумент перигея (градусы),

44-51 средняя аномалия (градусы),

53-63 среднее движение (обороты за сутки),

64-68 количество витков на эпоху,

69-69 контрольная сумма.

Пример:

LAGEOS

1 08820U 76039A 99305.13363095 .00000017 00000-0 -14899-1 0 5491

2 08820 109.8476 88.0433 0044671 225.0214 134.6681 6.38664538292505

В трёх строках формата НОРАД даны средние кеплеровские элементы 29250-го витка спутника Лагеос на 305.13363095 день с начала 1999 года. Угол наклонения равен 109.8476 градусов, долгота восходящего узла равна 88.0433 градуса, эксцентриситет орбиты равен 0.0044671, аргумент перигея и средняя аномалия в градусах равны 225.0214 и 134.6681, среднее движение составляет 6.38664538 оборотов за сутки.

Модуль для чтения данных имеет следующий вид.

UNIT UnrfNord ; // to read Norad format of the elements

INTERFACE // http://celestrak.com/NORAD/elements

Function ToGetNoradNumber ( snum : ShortString ) : LongInt ;

function ToGetLaunchNumber ( st : ShortString ) : LongInt ;

Procedure ToReadNoradFormat ( FileName : ShortString ) ;

Function BooFindElemWithNum ( FileName : ShortString ;

NumNorad : LongInt ) : Boolean ;

IMPLEMENTATION

Uses

UnConTyp , // types and constants , global variable SatElemCur

UnForDat , // for procedure TransDateToJD

UnForStr ; // for procedures ToEraseSomeBlanks ToReplaceBlanks

Function ToGetNoradNumber ( snum : ShortString ) : LongInt ;

Var

i,l : Byte ;

err : Integer ;

ch : Char ;

sa : ShortString ;

numl : LongInt ;

nums : LongInt ;

Begin

nums:=100000;

sa:=snum;

l:=Length(sa);

If l < 2

Then

Begin

ToGetNoradNumber:=nums;

Exit ;

end;

For i:=1 To l Do

Begin

ch:=sa[i];

If ch = ' ' Then ch:='1';

If NOT ( ch IN [ '0'..'9' ] )

Then

ch:='0';

sa[i]:=ch;

End;

Val(sa,numl,err);

If err = 0 Then nums:=numl;

ToGetNoradNumber:=nums;

End;

function ToGetPartNumber ( ca : Char ) : LongInt ;

begin

case ca of

'A' : ToGetPartNumber:=1;

'B' : ToGetPartNumber:=2;

'C' : ToGetPartNumber:=3;

'D' : ToGetPartNumber:=4;

'E' : ToGetPartNumber:=5;

'F' : ToGetPartNumber:=6;

'G' : ToGetPartNumber:=7;

'H' : ToGetPartNumber:=8;

'I' : ToGetPartNumber:=9

else ToGetPartNumber:=10;

end; // case

end;

function ToGetLaunchNumber ( st : ShortString ) : LongInt ;

var

sa : ShortString ;

ca : Char ;

ns : LongInt;

begin

sa:=Copy(st,1,5);

ca:=st[6];

ns:=100*ToGetNoradNumber(sa)+ToGetPartNumber(ca);

ToGetLaunchNumber:=ns;

end;

Procedure ToGetDateFromString ( st : ShortString ;

Var el : TElemRec ; // type UnCoTVar

Var b : Boolean ) ;

Var

ErrVal : Integer ;

Day,

Month,Year : Integer ;

PDay,d,v : Extended ;

sa,sy,sd,sn : ShortString ;

Begin

b:=FALSE;

sa:=Copy(st,3,5); // number of a satellite from norad

el.satnord:=ToGetNoradNumber(sa); // NORAD number

sa:=Copy(st,10,6); // number of satellite as date of launch

el.satnuml:=ToGetLaunchNumber(sa);

sy:=Copy(st,19,2);

sd:=Copy(st,21,12);

sn:=Copy(st,33,11);

ToReplaceBlanks(sy); // UnForStr

Val(sy,Year,ErrVal);

If ErrVal <> 0 Then Exit ;

If Year < 51 Then Year:=2000+Year;

If Year < 100 Then Year:=1900+Year;

Day:=31;

Month:=12;

Year:=Year-1;

TransDateToJD(Day,Month,Year,d); // UnForDat

ToReplaceBlanks(sd); // UnForStr

Val(sd,PDay,ErrVal);

If ErrVal <> 0

Then

Exit ;

d:=d+PDay;

el.t:=d; // julian date of elements

ToReplaceBlanks(sn); // UnForStr

Val(sn,v,ErrVal);

If ErrVal <> 0 Then v:=0.0;

el.dn:=v; // dn/dt n mean motion rev/day per day

b:=TRUE;

End;

Procedure ToGetElementsFromString ( st : ShortString ;

Var el : TElemRec ; // UnCoTVar

Var b : Boolean ) ;

Var

ErrVal : Integer ;

sc : ShortString ;

ac : Extended ;

Begin

b:=FALSE;

sc:=Copy(st,9,8); { inclination in degree }

ToReplaceBlanks(sc);

Val(sc,ac,ErrVal);

If ErrVal <> 0 Then Exit ;

el.ai:=ac;

sc:=Copy(st,18,8); { ascending node in degree }

ToReplaceBlanks(sc);

Val(sc,ac,ErrVal);

If ErrVal <> 0 Then Exit ;

el.au:=ac;

sc:=Copy(st,27,7); { pack eccentricity }

ToReplaceBlanks(sc);

sc:='0.'+sc;

Val(sc,ac,ErrVal);

If ErrVal <> 0 Then Exit ;

el.ae:=ac;

sc:=Copy(st,35,8); { argument of perigei in degree }

ToReplaceBlanks(sc);

Val(sc,ac,ErrVal);

If ErrVal <> 0 Then Exit ;

el.ao:=ac;

sc:=Copy(st,44,8); { mean anomaly in degree }

ToReplaceBlanks(sc);

Val(sc,ac,ErrVal);

If ErrVal <> 0 Then Exit ;

el.am:=ac;

sc:=Copy(st,53,11); { mean motion in rev/day }

ToReplaceBlanks(sc);

Val(sc,ac,ErrVal);

If ErrVal <> 0 Then Exit ;

el.an:=ac;

b:=TRUE;

End;

procedure ToReadNoradWithName ( FileName : ShortString ) ;

Var

TmpFile : TextFile ;

stwo : String[2] ;

StrN : String[16] ;

StrTmp : ShortString ;

ba,bb : Boolean ;

begin

AssignFile(TmpFile,FileName); // to assign file with this name

ReSet(TmpFile); // to open file

ReadLn(TmpFile,StrTmp); // to read file string by string

While NOT EOF(TmpFile) Do // if there is no end of file

Begin

StrN:=Copy(StrTmp,1,16); // this may be name of satellite

ReadLn(TmpFile,StrTmp);

stwo:=Copy(StrTmp,1,2);

If ( ( stwo = '1 ' ) // this is the first string

AND ( NOT EOF(TmpFile) ) ) // and the previous string is nullo

Then // and contains the name of satellite

Begin

SatElemCur.satname:=StrN; // this is previous nullo string

ToGetDateFromString(StrTmp,SatElemCur,ba); // the first string

ReadLn(TmpFile,StrTmp); // the second string

ToGetElementsFromString(StrTmp,SatElemCur,bb);

End;

CloseFile(TmpFile);

end;

// to read text file with the name FileName

Procedure ToReadNoradFormat ( FileName : ShortString ) ;

Begin

ToReadNoradWithName(FileName)

End;

Function BooFindElemWithNum ( FileName : ShortString ;

NumNorad : LongInt ) : Boolean ;

Var

TmpFile : TextFile ;

stwo : String[2] ;

StrN : String[16] ;

StrTmp : ShortString ;

ba,bb : Boolean ;

Begin

AssignFile(TmpFile,FileName); // to assign file with this name

ReSet(TmpFile); // to open file

SatElemCur.satnord:=0; // default to compare

ReadLn(TmpFile,StrTmp); // to read file string by string

While ( NOT EOF(TmpFile) ) // if there is no end of file

AND ( NumNorad <> SatElemCur.satnord ) Do // to compare

Begin

StrN:=Copy(StrTmp,1,16); // this may be name of satellite

ReadLn(TmpFile,StrTmp);

stwo:=Copy(StrTmp,1,2);

If ( ( stwo = '1 ' ) // this is the first string

AND ( NOT EOF(TmpFile) ) ) // and the previous string is nullo

Then // and contains the name of satellite

Begin

SatElemCur.satname:=StrN; // this is previous nullo string

ToGetDateFromString(StrTmp,SatElemCur,ba); // the first string

ReadLn(TmpFile,StrTmp); // the second string

ToGetElementsFromString(StrTmp,SatElemCur,bb);

End;

CloseFile(TmpFile);

BooFindElemWithNum:= ( NumNorad = SatElemCur.satnord ) ;

End;

END. // unit UnrfNord

Прогноз орбиты

Здесь представлены исходные тексты процедур для вычисления геоцентрических координат космического объекта в системе истинного экватора даты.

{ for the procedures

of the generalized two fixed centers problem }

UNIT UnitEilP ;

INTERFACE

Uses

UnConTyp ; // TypeDim3 , TElemRec , CoefType , PosVType

Procedure CRONA ( Var C : CoefType ) ;

Procedure ElemC ( ElemR : TElemRec ;

Var C : CoefType ) ;

Procedure ElemX ( TCur : Extended ; // julian date

ElemR : TElemRec ;

Var X : TypeDim3 ) ; // km

Procedure ALEI ( Var C : CoefType ) ;

Procedure XVEilR ( t : Extended ; // julian date

X : PosVType ; // km , km/s

Var ElemR : TElemRec ) ; // UnConTyp

Procedure AnomVM ( AE , AV : Extended ;

Var AM : Extended ) ;

Procedure KeplXV ( AB,AE,AI,AO,AU,AM : Extended ;

Var X : PosVType ) ; // km , km/s

IMPLEMENTATION

Uses

UnForFun ;

Procedure CRONA ( Var C : CoefType ) ;

Var

FM,CJ2,SIG,A,E,S,

O,O2,O3,O4,O5,E2,E4,

RE,S2,RS,S4,

P,P2,PS,P4,

Z2,PM,R1,R2,R3,R4,

F0,F1,F2,B0,B1,B2 : Extended ;

A1,A2,A3,

AB,ASZ,AC,

AD,ADZ,AX,AG,AT,

AG1,AG2,AN,AQ,

AC0,AC1,AC2,AC3,AC4,

AK1,AK2,AQ0,AQ1,AQ2,

AZ,AP0,AP1,AP2,AF0,

AF1,AF2,AF3,AF4,ALAM,

AN0,OH,AMU,AKF2,AKF4,

AKP2,AKP4 : Extended ;

Begin

FM:=C[1];

CJ2:=C[3];

SIG:=C[4];

A:=C[8];

E:=C[9];

S:=C[10];

AC:=SQRT(Abs(1.0e0-S*S));

If C[7] < 0.e0 Then AC:=-AC;

O:=1.0e0;

O2:=2.e0;

O3:=3.e0;

O4:=0.25e0;

O5:=0.5e0;

E2:=E*E;

RE:=O-E2;

E4:=E2*E2;

S2:=S*S;

RS:=O-S2;

S4:=S2*S2;

P:=CJ2/(A*RE);

P2:=P*P;

PS:=P*SIG;

Z2:=SIG*SIG;

AX:=E*(O+P2*RE*(O-O2*S2+P2*(O3

-16.e0*S2+14.e0*S4-O2*E2*RS*RS)));

C[14]:=AX;

AG:=-PS*(O-O2*S2-P2*(O3-12.e0*S2+10.e0*S4+E2*(O-O2*S4)));

AT:=PS*S*(O-P2*(5.e0-6.e0*S2-E2*(O-O2*S2)));

C[15]:=AG;

C[16]:=AT;

AK1:=0.125e0*P2*S2*(RE+Z2-4.e0*P2*RE*RS);

AK2:=0.125e0*P2*E2*(S2-P2*(O-10.e0*S2+11.0e0*S4+E2*S4));

C[26]:=AK1;

C[27]:=AK2;

PM:=SQRT(FM*A*RE);

P2:=O5*P2;

P4:=O5*P2*P2;

AG1:=PM*(O+P2*(RS*(O3+E2)+Z2*(6.e0-7.e0*S2))

-P4*RS*(9.e0+11.0e0*S2+E2*(6.e0+34.e0*S2)+E4*(O+O3*S2)));

AG2:=PM*(O-P2*(O3-4.e0*S2-E2)-P4*(9.e0-72.e0*S2

+64.e0*S4+E2*(O2-40.e0*S2+48.e0*S4)+E4));

AN:=AG1*(O+O2*AK2+9.e0*Sqr(AK2))/(AG2*(O+O2*AK1+9.e0*Sqr(AK1)));

C[17]:=AG1;

C[18]:=AG2;

C[19]:=AN;

PS:=PS*S*AC;

P2:=P*P;

P4:=P2/16.0e0;

C[20]:=O2*PS*(O-P2*(4.e0-5.e0*S2+E2*S2));

R1:=P2*AC;

R2:=O2*P4;

R3:=4.e0*R2;

R4:=O4*R1*P4;

AC0:=-R1*(O+O5*E2+P4*(24.e0-56.e0*S2

-4.e0*E2*(O+16.e0*S2)-E4*(O2+O3*S2)));

C[21]:=AC0;

C[22]:=-O2*R1*E*(O+R2*(4.e0-28.e0*S2-E2*(6.e0+7.e0*S2)));

C[23]:=-O4*R1*E2*(O-R3*(11.0e0+E2*(O+S2)));

C[24]:=16.e0*R4*E*E2*(O2-S2);

C[25]:=R4*E4*(O2+S2);

AQ0:=-O5*R1*(RE+O3*Z2-O2*P4*RE*(30.e0-35.e0*S2+E2*(O2+O3*S2)));

C[28]:=AQ0;

C[29]:=P2*PS*RE;

C[30]:=6.e0*R4*S2*RE*RE;

F2:=O-R3*RE*(4.e0-7.e0*S2)-P2*R2*RE

*(20.e0-136.e0*S2+113.e0*S4-E2*(8.e0+24.e0*S2-47.e0*S4));

F1:=P2*RE*(O2-O3*S2+R3*(O2-22.e0*S2+19.e0*S4

-E2*S2*(10.e0-13.e0*S2)));

F0:=O/(F1+F2);

F0:=F0/SQRT(O-AX*AX);

F0:=F0*RE*RE;

B0:=R3*(-S2+R3*(6.e0-20.e0*S2+15.e0*S4-E2*S4));

B1:=-P2*R3*S2*(O2-S2);

B2:=O3*P2*R2*S4;

AP0:=F0*(B0+B1+B2*(O+O5*AX*AX));

C[32]:=AP0;

C[33]:=AX*(B1+O2*B2)*F0;

C[34]:=O4*E2*B2*F0;

F0:=P2*F0*AG2/AG1;

AF0:=O5*F0*S2*(O+O3*AK1);

C[35]:=AF0;

C[36]:=-F0*S*(O2*AG-1.5e0*S*AT);

C[38]:=-F0*S2*AT/6.e0;

C[37]:=-F0*S2*(O4+AK1);

C[39]:=0.125e0*F0*S2*AK1;

C[40]:=-AP0-AN*AF0;

C[45]:=AN-O;

C[42]:=C[41]/(O-C[40]);

C[43]:=C[45]*C[42];

C[46]:=AC0+AN*AQ0;

C[44]:=C[42]*C[46];

C[47]:=AK1*(O+4.e0*AK1);

C[48]:=-0.75e0*AK1*AK1;

C[49]:=-AN*AK2*(O+4.e0*AK2);

C[50]:=0.75e0*AN*AK2*AK2;

C[23]:=C[23]+AQ0*C[49];

C[30]:=C[30]+AQ0*C[47];

C[34]:=C[34]+AF0*C[49];

C[37]:=C[37]+AF0*C[47];

End ;

Procedure ElemC ( ElemR : TElemRec ; // type TElemRec from UnConTyp

Var C : CoefType ) ;

Var

j : Byte ;

FM,RZ,CJ2,SIG,

S,ASZ,AC,AD,ADZ,F0,F1,F2,FC,

RE,PE,P2,P4,PS,D2,RD,D4,RC,RR,

R1,F3,AB,DZ2,RDZ,DZ4,RCZ,RRZ,

F4,F5,F6,F7,F8,F9,R3,R4,R2,

A1,A2,A3,AN0 : Extended ;

Begin

FM:= 3.98600448000e+5; // km^3/s^2

RZ:= 6.37813700000e+3; // km

CJ2:= 2.09728850857e+2; // km

SIG:=-3.55889244242e-02;

AN0:=TwoPi*ElemR.an/86400.0; // to radian per second

S:=GraRad*ElemR.ai; // constants from unit UnConTyp

ASZ:=SIN(S);

AC:=COS(S);

AD:=ASZ;

ADZ:=-ASZ;

A1:=-Exp(2.0e0/3.0e0*Ln(FM*AN0)); // km^2/s^2

F0:=-FM/A1; // km

AB:=F0; // km

FC:=CJ2/F0;

RE:=1.e0-ElemR.ae*ElemR.ae;

F1:=2.e0-RE;

F2:=RE*RE;

For j:=1 To 4 Do

Begin

PE:=CJ2/(AB*RE);

P2:=PE*PE;

P4:=P2*P2;

PS:=PE*SIG;

D2:=AD*AD;

RD:=1.e0-D2;

D4:=D2*D2;

RC:=1.e0-2.e0*PS*AD-P2*D2*RE;

RR:=1.e0+2.e0*P2*D2*F1+P4*D4*F2;

R1:=RC/RR;

F3:=1.e0-RD*R1;

AB:=F0*(1.e0-P2*RE*RD*R1);

DZ2:=ADZ*ADZ;

RDZ:=1.e0-DZ2;

DZ4:=DZ2*DZ2;

RCZ:=1.e0-2.e0*PS*ADZ-P2*DZ2*RE;

RRZ:=1.e0+2.e0*P2*DZ2*F1+P4*DZ4*F2;

R2:=RCZ/RRZ;

F4:=1.e0-RDZ*R2;

F5:=F3+D2*ADZ*(FC*PE*(2.e0*AD+ADZ)+2.e0*PS);

F6:= F4+DZ2*AD*(FC*PE*(2.e0*ADZ+AD)+2.e0*PS);

F7:=SQRT(ABS((F5/F6)*(F4/F3)));

F9:=R1/R2;

F8:=2.e0*F7*(1.e0+F7)*F9*ASZ;

ADZ:=(-1.e0+F9-F9*Sqr((1.e0+F7)*ASZ)+(1.e0-F9*F7*F7)*DZ2)/F8;

AD:=ASZ*(1.e0+F7)+F7*ADZ;

End;

R3:=FM*AB*RE*RD*(1.e0+2.e0*P2*F1+P4*F2)*R1;

R4:=FM*AB*RE*RDZ*(1.e0+2.e0*P2*F1+P4*F2)*R2;

A3:=SQRT(ABS(R3));

IF AC < 0.e0 Then A3:=-A3;

A2:=-FM*AB*RE*(1.e0+2.e0*P2*F1*RD*R1+P4*F2*RD*R1);

C[8]:=AB;

C[9]:=ElemR.ae;

C[10]:=ASZ;

C[11]:=AC;

C[5]:=A1;

C[6]:=A2;

C[7]:=A3;

C[12]:=AD;

C[13]:=ADZ;

C[1]:=FM;

C[2]:=RZ;

C[3]:=CJ2;

C[4]:=SIG;

C[31]:=-C[8]*C[9]*A1/FM;

C[41]:=AN0;

CRONA(C);

End;

Procedure ElemX ( TCur : Extended ;

ElemR : TElemRec ;

Var X : TypeDim3 ) ; // type from UnConTyp

Var

j : Byte ;

C : CoefType ;

O,W,OL,

A1,A2,A3,

AB,AE,ASZ,AC,

AD,ADZ,AX,AG,AT,

AG1,AG2,AN,AQ,

AC0,AC1,AC2,AC3,AC4,

AK1,AK2,AQ0,AQ1,AQ2,

AZ,AP0,AP1,AP2,AF0,

AF1,AF2,AF3,AF4,ALAM,

AN0,OH,AMU,AKF2,AKF4,

AKP2,AKP4,RXX,RX,PUDL,

SPP,CPP,SPR,RB,

SFF,CFF,SFR,EE,RT,BT,

SW,CW,SU,CU : Extended ;

X1,X2,X3,V1,V2,V3,

ZX,U2,R1,F2,BS,

SE,CE,EA,SP,CP,EP,DE,

SF,CF,EF,S2P,C2P,S4P,

S2F,C2F,S4F,AO0,BW,

AU0,AM0,R2,R3,R4,R5,R6,

FM,RZ,CJ2,SIG,Z2,Z0 : Extended ;

UDL,UDG,UDH,DifT,ddl : Extended ;

Begin

ElemC(ElemR,C);

ddl:=TwoPi*0.5*ElemR.dn*Sqr(TCur-ElemR.t); // +1/2 dn/dt (t-t_0)^2 radian

DifT:=86400.0*(TCur-ElemR.t); { in second }

UDL:=GraRad*ElemR.am+C[42]*DifT+ddl; { in radian }

UDG:=GraRad*ElemR.ao+C[43]*DifT+HalfPi;

UDH:=GraRad*ElemR.au+C[44]*DifT-HalfPi;

O:=1.0e0;

W:=2.00000e0;

OL:=0.0e0;

FM:=C[1];

RZ:=C[2];

CJ2:=C[3];

SIG:=C[4];

A1:=C[5];

A2:=C[6];

A3:=C[7];

AB:=C[8];

AE:=C[9];

ASZ:=C[10];

AC:=C[11];

AD:=C[12];

ADZ:=C[13];

AX:=C[14];

AG:=C[15];

AT:=C[16];

AG1:=C[17];

AG2:=C[18];

AN:=C[19];

AQ:=C[20];

AC0:=C[21];

AC1:=C[22];

AC2:=C[23];

AC3:=C[24];

AC4:=C[25];

AK1:=C[26];

AK2:=C[27];

AQ0:=C[28];

AQ1:=C[29];

AQ2:=C[30];

AZ:=C[31];

AP0:=C[32];

AP1:=C[33];

AP2:=C[34];

AF0:=C[35];

AF1:=C[36];

AF2:=C[37];

AF3:=C[38];

AF4:=C[39];

ALAM:=C[40];

AN0:=C[41];

OH:= C[45];

AMU:=C[46];

AKF2:=C[47];

AKF4:=C[48];

AKP2:=C[49];

AKP4:=C[50];

RXX:=1.e0-AX*AX;

RX:=SQRT(RXX);

Z2:=Sqr(CJ2);

Z0:=CJ2*SIG;

SP:=Sin(UDL);

SE:=SP;

CP:=Cos(UDL);

CE:=CP;

SF:=OL;

CF:=O;

EP:=UDL;

EA:=EP;

PUDL:=EP-UDL;

SPP:=W*SP*CP;

CPP:=CP*CP-SP*SP;

SPR:=W*SPP*CPP;

SFF:=W*SF*CF;

j:=0;

DE:=1.0e0;

While ( ( j < 15 ) AND ( DE > 1.0e-15 ) ) Do

Begin

j:=j+1;

PUDL:=EP-UDL;

SPP:=W*SP*CP;

CPP:=CP*CP-SP*SP;

SPR:=W*SPP*CPP;

SFF:=W*SF*CF;

CFF:=CF*CF-SF*SF;

SFR:=W*SFF*CFF;

EE:=UDL+AZ*SE+ALAM*PUDL-AP1*SP-AP2*SPP

-AF1*SF+AF2*SFF+AF3*(SF*CFF+CF*SFF)-AF4*SFR;

SE:=Sin(EE);

CE:=Cos(EE);

RB:=O/(O-AX*CE);

SP:=RX*RB*SE;

CP:=(CE-AX)*RB;

EP:=EE+ArcTan((SP*CE-CP*SE)/(CP*CE+SP*SE));

EF:=EP+UDG+OH*PUDL+AKP2*SPP+AKP4*SPR-AKF2*SFF+AKF4*SFR;

SF:=Sin(EF);

CF:=Cos(EF);

DE:=Abs(EA-EE);

EA:=EE;

End;

BS:=AB*(O-AE*CE);

RT:=O/(O-AT*CF);

BT:=RT*(-ASZ*CF+AG);

R1:=BS*BS+Z2;

R2:=O-BT*BT;

R3:=SQRT(R1*R2);

R4:=-AC*CF+AQ;

R6:=O/(R4*R4+SF*SF);

R5:=SQRT(R6);

BW:=UDH+1.570796326794896e0+AMU*PUDL+AQ1*SF-AQ2*SFF;

BW:=BW+AC1*SP+AC2*SPP+AC3*(SP*CPP+CP*SPP)+AC4*SPR;

SW:=R4*R5;

CW:=SF*R5;

SU:=SIN(BW);

CU:=COS(BW);

X[1]:=R3*(CW*CU-SW*SU);

X[2]:=R3*(CW*SU+SW*CU);

X[3]:=Z0+BS*BT;

End;

Procedure ALEI ( Var C : CoefType ) ; // type from UnConTyp

Var

j : Byte ;

O,A1,A2,A3,FM,CJ2,SIG,

AB,AE,AD,ADZ,RE,FE,PA,

FA,RF,RT,RD,RDZ,DZ2,

D2,PE,P1,P4,P2,RR,RC,

RRZ,RCZ,R,FD,PS,E2 : Extended ;

Begin

O:=1.0e0;

A1:=C[5];

A2:=C[6];

A3:=C[7];

FM:=C[1];

CJ2:=C[3];

SIG:=C[4];

AB:=-FM/A1;

FA:=AB;

FE:=-A2/FM;

RE:=FE/FA;

PA:=A3*A3;

RD:=-PA/A2;

RDZ:=RD;

FD:=RD;

PA:=PA/FM;

For j:=1 To 5 Do

Begin

PE:=CJ2/(AB*RE);

P2:=PE*PE;

PS:=PE*SIG;

D2:=O-RD;

AD:=SQRT(Abs(D2));

P1:=2.e0*P2*(2.e0-RE);

P2:=P2*RE;

P4:=P2*P2;

RR:=O+P1*D2+P4*D2*D2;

RC:=O-2.e0*PS*AD-P2*D2;

RR:=RC/RR;

RF:=(O+P1+P4)*RR;

RT:=O+(RF-RR)*RD;

AB:=FA*(O-P2*RD*RR);

RE:=FE/(AB*RT);

RD:=FD*RT/RF;

DZ2:=O-RDZ;

ADZ:=-SQRT(Abs(DZ2));

RRZ:=O+P1*DZ2+P4*DZ2*DZ2;

RCZ:=O-2.e0*PS*ADZ-P2*DZ2;

RDZ:=PA*RRZ/(AB*RE*(O+P1+P4)*RCZ);

End;

E2:=O-RE;

C[8]:=AB;

C[9]:=SQRT(Abs(E2));

D2:=O-RD;

C[12]:=SQRT(Abs(D2));

DZ2:=O-RDZ;

C[13]:=-SQRT(Abs(DZ2));

P2:=PE*PE*(3.e0-4.e0*D2+E2);

C[10]:=C[12]+PS*RD*(O-3.e0*PS*C[12]-P2);

C[11]:=SQRT(Abs(O-C[10]*C[10]));

IF A3 < 0.e0 Then C[11]:=-C[11];

RF:=O/FM;

R:=-A1*A1*A1;

C[41]:=SQRT(R)*RF;

C[31]:=-C[8]*C[9]*A1*RF;

End;

Procedure XVEilR ( t : Extended ;

X : PosVType ;

Var ElemR : TElemRec ) ; // type from UnConTyp

Var

C : CoefType ;

Z4,FI,VT,VS : Extended ;

A1,A2,A3,

AB,ASZ,AC,

AD,ADZ,AX,AG,AT,

AG1,AG2,AN,AQ,

AC0,AC1,AC2,AC3,AC4,

AK1,AK2,AQ0,AQ1,AQ2,

AZ,AP0,AP1,AP2,AF0,

AF1,AF2,AF3,AF4,ALAM,

AN0,OH,AMU,AKF2,AKF4,

AKP2,AKP4,PUDL : Extended ;

X1,X2,X3,V1,V2,V3,

ZX,U2,R1,F2,BS,BT,

SE,CE,EA,SP,CP,EP,

SF,CF,EF,SPP,CPP,SPR,

SFF,CFF,SFR,AO0,BW,

AU0,AM0,R2,SR,CR,

FM,RZ,CJ2,SIG,Z2,Z0 : Extended ;

AOR,AUR,AMR : Extended ; // in radian

AM,AO,AU,AI,AE,ANOB : Extended ;

Begin

FM:= 3.98600448000e+5; // km^3/S^2

RZ:= 6.37813700000e+3; // km

CJ2:= 2.09728850857e+2; // km

SIG:=-3.55889244242e-02;

Z2:=Sqr(CJ2);

Z0:=CJ2*SIG;

X1:=X[1];

X2:=X[2];

X3:=X[3]-Z0;

V1:=X[4];

V2:=X[5];

V3:=X[6];

ZX:=X3*X3;

R2:=X1*X1+X2*X2+ZX;

U2:=V1*V1+V2*V2+V3*V3;

R1:=X1*V1+X2*V2+X3*V3;

F2:=0.5e0*(R2-Z2+Sqrt(Sqr(R2-Z2)+4.0e0*Z2*ZX));

BS:=SQRT(F2);

BT:=X3/BS;

Z4:=Z2*Sqr(BT);

FI:=2.0e0*FM/(F2+Z4);

VT:=BS*V3-BT*R1;

VS:=BS*R1+Z2*BT*V3;

A1:=U2-FI*(BS-Z0*BT);

A3:=X1*V2-V1*X2;

A2:=R1*R1-R2*U2+Z2*V3*V3-FI*(BS*Z4+Z0*BT*F2);

C[1]:=FM; C[2]:=RZ; C[3]:=CJ2; C[4]:=SIG;

C[5]:=A1; C[6]:=A2; C[7]:=A3;

ALEI(C);

CRONA(C);

AB:=C[8]; AE:=C[9];

ASZ:=C[10]; AC:=C[11];

AD:=C[12]; ADZ:=C[13];

AX:=C[14]; AG:=C[15];

AT:=C[16]; AG1:=C[17];

AG2:=C[18]; AN:=C[19];

AQ:=C[20]; AC0:=C[21];

AC1:=C[22]; AC2:=C[23];

AC3:=C[24]; AC4:=C[25];

AK1:=C[26]; AK2:=C[27];

AQ0:=C[28]; AQ1:=C[29];

AQ2:=C[30]; AZ:=C[31];

AP0:=C[32]; AP1:=C[33];

AP2:=C[34]; AF0:=C[35];

AF1:=C[36]; AF2:=C[37];

AF3:=C[38]; AF4:=C[39];

ALAM:=C[40]; AN0:=C[41]; OH:= C[45]; AMU:=C[46];

AKF2:=C[47]; AKF4:=C[48]; AKP2:=C[49]; AKP4:=C[50];

CE:=(1.0e0-BS/AB)/AE;

SE:=SQRT(1.0e0-CE*CE);

IF VS < 0.0e0 Then SE:=-SE;

EA:=DATAN2(SE,CE);

CP:=(CE-AX)/(1.0e0-AX*CE);

SP:=Sqrt(1.0e0-CP*CP);

IF VS < 0.e0 Then SP:=-SP;

EP:=DATAN2(SP,CP);

SR:=(BT-AG)/(ASZ-BT*AT);

CR:=Sqrt(1.0e0-SR*SR);

IF VT < 0.e0 Then CR:=-CR;

SF:=CR; CF:=-SR;

EF:=DATAN2(SF,CF);

SPP:=2.e0*SP*CP; CPP:=CP*CP-SP*SP;

SPR:=2.e0*SPP*CPP; SFF:=2.e0*SF*CF;

CFF:=CF*CF-SF*SF; SFR:=2.e0*SFF*CFF;

AM0:=EA-AZ*SE-ALAM*EP+AP1*SP+AP2*SPP

+AF1*SF-AF2*SFF-AF3*(SF*CFF+CF*SFF)+AF4*SFR;

AMR:=AM0/(1-ALAM);

PUDL:=EP-AMR;

AO0:=EF-EP-OH*PUDL-AKP2*SPP-AKP4*SPR+AKF2*SFF-AKF4*SFR;

AOR:=AO0-PiTwo;

R1:=-AC*CF+AQ;

BW:=DATAN2(X2,X1)-DATAN2(R1,SF);

AU0:=BW-PiTwo-AMU*PUDL

-AC1*SP-AC2*SPP-AC3*(SPP*CP+CPP*SP)

-AC4*SPR-AQ1*SF+AQ2*SFF;

AUR:=AU0+PiTwo;

AM:=RadGra*AMR; // mean anomaly in degree

If AM < 0.0 Then AM:=360.0+AM; If AM > 360.0 Then AM:=AM-360.0;

AO:=RadGra*AOR; // argument of perigei in degree

If AO < 0.0 Then AO:=360.0+AO; If AO > 360.0 Then AO:=AO-360.0;

AU:=RadGra*AUR; // the ascending node in degree

If AU < 0.0 Then AU:=360.0+AU; If AU > 360.0 Then AU:=AU-360.0;

AI:=RadGra*DATAN2(ASZ,AC); // function from UnForFun

ANOB:=(AN0/(1-ALAM))*(RADGRA/360.e0)*8.64e4; // mean motion rev. per day

ElemR.t:=t;

ElemR.ao:=AO; // argument of perigei in degree

ElemR.au:=AU; // the ascending node in degree

ElemR.ai:=AI; // inclination in degree

ElemR.ae:=AE;

ElemR.am:=AM; // mean anomaly in degree

ElemR.an:=ANOB; // mean motion in revolution per day

End;

Procedure AnomVM ( AE , AV : Extended ;

Var AM : Extended ) ;

Var

V,CV,SV,E,U,

R,SE,CE,RS,RC : Extended;

Begin

V:=GRARAD*AV;

CV:=COS(V);

SV:=SIN(V);

R:=1.0e0/(1.0e0+AE*CV);

SE:=SQRT(1.0e0-AE*AE)*SV*R;

CE:=(CV+AE)*R;

RS:=SE*CV-CE*SV;

RC:=CE*CV+SE*SV;

E:=V+DATAN2(RS,RC); // function from UnForFun

U:=E-AE*SE;

AM:=RADGRA*U;

If AM < 0 Then AM:=360+AM; // mean anomaly in degree

End ;

Procedure KeplXV ( AB,AE,AI,AO,AU,AM : Extended ;

Var X : PosVType ) ;

Var

j : Byte ;

FM,S,C,O,R,P,Q,

VR,VN,SV,CV,SB,CB,SU,CU,

V,B,U,W,E,SE,CE,RR,RS,RC : Extended ;

Begin

FM:=3.9860047E5; // km^3/s^2

S:=SIN(GRARAD*AI); C:=COS(GRARAD*AI);

O:=GRARAD*AO; U:=GRARAD*AU; W:=GRARAD*AM;

E:=W;

For j:=1 To 10 Do E:=W+AE*SIN(E);

CE:=COS(E); SE:=SIN(E);

RR:=AB*(1.0e0-AE*CE); // km

R:=1.0e0/(1.0e0-AE*CE);

SV:=SQRT(1.0e0-AE*AE)*SE*R; CV:=(CE-AE)*R;

RS:=SV*CE-CV*SE; RC:=CV*CE+SV*SE;

V:=E+DATAN2(RS,RC);

B:=O+V;

SB:=SIN(B); CB:=COS(B);

SU:=SIN(U); CU:=COS(U);

X[1]:=RR*(CB*CU-SB*SU*C); // km

X[2]:=RR*(CB*SU+SB*CU*C); // km

X[3]:=RR*SB*S; // km

P:=AB*(1.0e0-AE*AE); // km

Q:=SQRT(FM/P);

VR:=Q*AE*SV/RR;

VN:=Q*(1.0e0+AE*CV);

X[4]:=VR*X[1]+VN*(-SB*CU-CB*SU*C); // km/s

X[5]:=VR*X[2]+VN*(-SB*SU+CB*CU*C); // km/s

X[6]:=VR*X[3]+VN*CB*S; // km/s

End ;

END. // unit UnitEilP

Решение задачи

Основной модуль программного приложения содержится в файле ‘SatelEph.dpr’ и имеет вид

program SatelEph;

{$APPTYPE CONSOLE}

uses

SysUtils ,

UnConTyp , // for PrtFile

uncomput in 'uncomput.pas';

begin

AssignFile(PrtFile,'sateleph.prt');

ReWrite(PrtFile);

ForSimpleSatelEph ; // the main procedure from unit UnComPut

CloseFile(PrtFile);

end.

Результаты расчётов будут сохранены в наборе данных ‘SatelEph.prt’.

В основном модуле вызывается основная подпрограмма для ввода исходных данных и проведения вычислений. Эта процедура определена в следующей программной единице.

unit uncomput; // main unit for calculating

interface

procedure ForSimpleSatelEph ;

implementation

uses

UnConTyp , // for constants types and global variables PrtFile

UnForCoo , // for procedure GeoSpherFromCart

UnForDat , // for procedure TransDateToJD

UnForSat , // for procedure SatHorizon

UnForStr , // function IntegerToFixStr

UnRefrOl , // for procedure SimpleDataForRefraction

UnrfNord ; // function BooFindElemWithNum

procedure PointCoor ( n : Integer ) ;

var

f : TextFile ; // file with sites positions pointcoo.txt

k : Integer ;

a : Extended ;

b : Extended ;

c : Extended ;

s : ShortString ;

begin

AssignFile(f,'pointcoo.txt');

ReSet(f);

ReadLn(f); // the first text string

ReadLn(f); // the second text string

while not EOF(f) do

begin

ReadLn(f,k,a,b,c,s);

if k = n // to compare numbers

then // this is our point

with PlaceCoor do // global variable PlaceCoor from UnConTyp

begin

num:=k;

name:=s;

x:=a; // x position in meter

y:=b; // y position in meter

z:=c; // z position in meter

GeoSpherFromCart(x,y,z,f,l,h); // UnForCoo

end;

CloseFile(f);

with placecoor do

begin

WriteLn(PrtFile,'point ',name);

WriteLn(PrtFile,x:15:3,y:15:3,z:15:3);

WriteLn(PrtFile,f:15:6,l:15:6,h:15:3);

end;

procedure NoradFile ( s : ShortString ) ; // name of file to stNoradFilw

var

n : Byte ;

a : ShortString ;

begin

a:=s; // in string may be not only name of file

n:=Pos(' ',a); // to delete forward blanks

while n = 1 do

begin

Delete(a,1,1); // to delete forward blanks

n:=Pos(' ',a); // n = 1 for froward blanks

end;

if n > 1 // may be unusefull information after blanks

then // to truncate this part of string

a:=Copy(a,1,n-1);

stNoradFile:=a; // global ShortString form UnConTyp

WriteLn(PrtFile,'name of NORAD file ',stNoradFile);

end;

procedure BooReadStartFile ; // file forsolve.txt

var // file with the data exists

f : TextFile ;

n : Integer ; // number of point

k : LongInt ; // NORAD number of satellite

d : Integer ; // day

m : Integer ; // month

y : Integer ; // year

p : Integer ; // step to print

s : ShortString ; // name of file with elements

begin

AssignFile(f,'forsolve.txt');

ReSet(f);

ReadLn(f); // the first text string

ReadLn(f); // the second text string

ReadLn(f,n); // number of point for observations

PointCoor(n); // to define point and to fill global var PlaceCoor

ReadLn(f,s); // for name of file with elements

NoradFile(s); // name of elements file to stNoradFile

ReadLn(f,k); // number of satellite for ephemeris

BooFindElemWithNum(stNoradFile,k); // UnrfNord to fill global var SatElemCur

ReadLn(f,d,m,y); // for start date

TransDateToJD(d,m,y,StartDate); // proc UnForDat var UnConTyp

ReadLn(f,d,m,y); // for finish date

TransDateToJD(d,m,y,FinishDate); // proc UnForDat var UnConTyp

ReadLn(f,p); // for step of priniting in second

StepEphem:=p/86400.0; // step of printing in day

Close(f); // this is the end of file forsolve.txt

with SatElemCur do

begin

WriteLn(PrtFile,satname:20,satnord:10,satnuml:10);

WriteLn(PrtFile,t:15:6);

WriteLn(PrtFile,ae:12:8,ai:12:6,an:15:10);

end;

procedure SatInform ( t , az , al , ro : Extended ) ; // to PrtFile

var

d,m,y : Integer ; // day month year

h,min : Integer ; // hour minute

ps,pp : Extended ; // second as part of day and pure second

begin

TransJDtoDate(t,d,m,y,ps); // UnForDat with part of day in second

FromSecPartToHourMinSec(ps,h,min,pp); // UnForDat pure time

WriteLn(PrtFile,IntegerToFixStr(y,4):5, // year function from UnForStr

IntegerToFixStr(m,2):3, // month

IntegerToFixStr(d,2):3, // day

IntegerToFixStr(h,2):5, // hour

IntegerToFixStr(min,2):3, // minute

FloatToFixStr(pp,6,3):7, // second

az:12:6,al:12:6,ro:12:3); // azimuth altitude range

end;

procedure ForSimpleSatelEph ;

var

az,al,ro : Extended ; // azimut altitude

begin

SimpleDataForRefraction; // UnRefrOl

BooReadStartFile; // to read file forsolve.txt

CurrentDate:=StartDate; // julian day

while CurrentDate < FinishDate do

begin

SatHorizon(CurrentDate,SatElemCur,PlaceCoor,az,al,ro); // UnForSat

if al > 10.0

then // satellite above horizon

SatInform(CurrentDate,az,al,ro); // to PrtFile

CurrentDate:=CurrentDate+StepEphem; // the next step

end;

end.

Результаты вычислений с исходными данными, приведёнными выше, имеют следующий вид.

point Base

2886365.206 2155941.870 5245817.642

55.699306 36.757528 198.000

name of NORAD file geosat04.txt

STELLA 22824 9306102

2453115.281324

0.00078370 98.256300 14.2725691400

2004 04 22 02 48 00.000 41.255753 13.157505 2166.837

2004 04 22 02 49 00.000 51.165652 17.479755 1910.123

2004 04 22 02 50 00.000 64.233406 21.464420 1715.791

2004 04 22 02 51 00.000 80.472760 24.049352 1607.643

2004 04 22 02 52 00.000 98.214623 24.122307 1603.803

2004 04 22 02 53 00.000 114.626186 21.636321 1705.129

2004 04 22 02 54 00.000 127.908858 17.668551 1894.603

2004 04 22 02 55 00.000 137.985782 13.313253 2148.525

2004 04 22 04 27 00.000 11.931128 12.578454 2206.485

2004 04 22 04 28 00.000 10.136180 19.197010 1823.525

2004 04 22 04 29 00.000 6.970522 28.241252 1461.444

2004 04 22 04 30 00.000 0.361029 41.410155 1142.421

2004 04 22 04 31 00.000 341.054260 59.934979 914.632

2004 04 22 04 32 00.000 271.766306 69.067149 855.585

2004 04 22 04 33 00.000 228.305275 51.384726 996.001

2004 04 22 04 34 00.000 216.809478 34.998861 1271.034

2004 04 22 04 35 00.000 212.102634 23.839760 1612.483

2004 04 22 04 36 00.000 209.619755 15.988229 1985.404

2004 04 22 04 37 00.000 208.117173 10.074260 2373.893

2004 04 22 06 08 00.000 346.091499 12.415791 2217.468

2004 04 22 06 09 00.000 335.089577 15.417243 2027.090

2004 04 22 06 10 00.000 321.958427 17.428114 1913.594

2004 04 22 06 11 00.000 307.661308 17.834948 1891.301

2004 04 22 06 12 00.000 293.811617 16.488153 1963.350

2004 04 22 06 13 00.000 281.753757 13.843938 2119.819

2004 04 22 06 14 00.000 271.963867 10.565738 2343.254

2004 04 22 12 37 00.000 87.337228 10.614430 2338.322

2004 04 22 12 38 00.000 77.462204 13.850848 2118.203

2004 04 22 12 39 00.000 65.347388 16.426110 1965.794

2004 04 22 12 40 00.000 51.507999 17.686520 1898.213

2004 04 22 12 41 00.000 37.302390 17.205910 1924.751

2004 04 22 12 42 00.000 24.308371 15.159059 2041.693

2004 04 22 12 43 00.000 13.439493 12.157551 2234.513

2004 04 22 14 14 00.000 151.385898 10.315530 2354.032

2004 04 22 14 15 00.000 149.776230 16.298698 1966.262

2004 04 22 14 16 00.000 147.108481 24.262295 1594.751

2004 04 22 14 17 00.000 142.019683 35.596240 1256.132

2004 04 22 14 18 00.000 129.456951 52.117472 986.873

2004 04 22 14 19 00.000 83.060599 68.744011 856.593

2004 04 22 14 20 00.000 18.075929 58.427091 926.421

2004 04 22 14 21 00.000 359.697311 40.380593 1160.544

2004 04 22 14 22 00.000 353.158747 27.577523 1482.409

2004 04 22 14 23 00.000 349.972708 18.736844 1845.753