2юнит:
PHP код:
unit Unit2_1512;
interface
uses
Messages;
const
(**) TwoPi = 2 * Pi;
RE = 6371.21; //средний радиус Земли км
ToRad = Pi/180; // Deg2Rad =
ToGrad = 180/Pi; // Rad2Deg
Mu = 0.398603e6; // ГПЗ
(**) We = 7.2921151467e-5; // Угловая скорость вращения Земли (ср солн сут)
(**) JD2000 = 2451545.0; // юлианская дата эпохи 2000 г., январь 1.5. (1 января 12 час)
wm_DrawChart = wm_User + 1;
type
T3DVector = record
X, Y, Z : real;
end;
T3DMatrix = array[1..3, 1..3] of real;
TGeo = record
Pos : T3DVector; // положение ГСК
R, Longitude, Latitude,alfa : real // радиус, долгота, широта
end;
TClassic = record
Rp, e, ArgLat, Incl, AscNode, ArgPerigee: real
end;
TCartesian = record
Pos: T3DVector; // АГЭСК ( положение )
Vel : T3DVector; // АГЭСК ( скорость )
end;
POrbit = ^TOrbit; {New}
TOrbit = record
Epoch : TDateTime; // Дата и время начальных условий
Time : real; // Собственное время КА (c)
Classic : TClassic;
Cartes : TCartesian;
Geo : TGeo // ГСК ( положение, долгота и широта )
end;
// Векторная алгебра
function ScalProd( A, B : T3DVector ) : real;
function Summa( A, B : T3DVector ) : T3DVector;
function Differ( A, B : T3DVector ) : T3DVector;
function Modul( A : T3DVector ): real;
function UnitVect ( A : T3DVector ): T3DVector;
function InitVect( X, Y, Z : real) : T3DVector;
function ZeroVect : T3DVector;
function MxV( M : T3DMatrix; V : T3DVector) : T3DVector;
(**) function MAbs_Geo( Value : real ) : T3DMatrix; // матрица перехода АГЭСК->ГСК
// Вычисление декартовых координат КА по классическим элементам
procedure ClassToXYZ( var AOrbit : TOrbit );
// Вычисление координат КА в земной СК
// Вычисление истинной аномалии КА
function TrueAnomaly( AOrbit : TOrbit ) : real;
// Вычисление параметра эллипса р
function Parameter( AOrbit : TOrbit ) : real;
// Вычисление радиус-вектора КА (r)
function Radius( AOrbit : TOrbit ) : real;
// Вычисление среднего движения (n)
function MeanMotion( AOrbit : TOrbit ) : real;
// Вычисление периода обращения (Т)
function Period( AOrbit : TOrbit ) : real;
// Вычисление радиальной скорости КА Vr
function VRadial( AOrbit : TOrbit ) : real;
// Вычисление большой полуоси орбиты (км)
function SemiMajorAxis( AOrbit : TOrbit ) : real;
// Вычисление эксцентрической аномалии (рад)
function EccentrAnomaly ( AOrbit : TOrbit ): real;
// Вычисление средней аномалии (рад)
function MeanAnomaly ( AOrbit : TOrbit ): real;
// Вычисление момента времени прохождения перигея (с)
function PerigeePassTime( AOrbit : TOrbit ) : real;
function KeplerEquation(const e,M : real):real;
procedure SetTime( ATime : real; var AOrb : TOrbit); overload;
procedure SetTime( ATime : real; IniOrb : TOrbit; var AOrb : TOrbit); overload;
function SetTimeF( ATime : real; IniOrb : TOrbit ) : TOrbit;
(**)function SiderialTime0( Value : TDateTime ) : real; // 0h дата
(**)function SiderialTime ( Value : TDateTime ) : real; // UTC дата+время
//procedure urTrassa(Value:TDateTime ; Orbit:TOrbit; var lat,long,long_ekv:double);
procedure NIP (var Nip:TGeo);
procedure ugolMesta (var alfa:real);
procedure DecKoord(var Nip:TGeo);
implementation
uses
Math,
SysUtils,
(**) DateUtils; {!!!}
var
Pos,razn : T3DVector;
Orb : TClassic;
function ScalProd( A, B : T3DVector ) : real;
begin
Result := A.X * B.X + A.Y * B.Y + A.Z * B.Z;
end;
function Summa( A, B : T3DVector ) : T3DVector;
begin
Result.X := A.X + B.X;
Result.Y := A.Y + B.Y;
Result.Z := A.Z + B.Z;
end;
function Differ( A, B : T3DVector ) : T3DVector;
begin
Result.X := A.X - B.X;
Result.Y := A.Y - B.Y;
Result.Z := A.Z - B.Z;
end;
function Modul( A : T3DVector ): real;
begin
Result := Sqrt( Sqr(A.X) + Sqr(A.Y) + Sqr(A.Z));
end;
function UnitVect ( A : T3DVector ): T3DVector;
var
R : real;
begin
R := Modul(A);
if R<>0 then
begin
R := 1/R;
Result.X := A.X * R;
Result.Y := A.Y * R;
Result.Z := A.Z * R;
end else
begin
end;
end;
function InitVect( X, Y, Z : real) : T3DVector;
begin
Result.X := X;
Result.Y := Y;
Result.Z := Z;
end;
function ZeroVect : T3DVector;
begin
Result := InitVect(0,0,0);
end;
function MxV( M : T3DMatrix; V : T3DVector) : T3DVector;
begin
Result.X := M[1,1] * V.X + M[1,2] * V.Y + M[1,3] * V.Z;
Result.Y := M[2,1] * V.X + M[2,2] * V.Y + M[2,3] * V.Z;
Result.Z := M[3,1] * V.X + M[3,2] * V.Y + M[3,3] * V.Z;
end;
function MAbs_Geo( Value : real ) : T3DMatrix; // матрица перехода АГЭСК->ГСК
begin
FillChar( Result, SizeOf(T3DMatrix), 0);
Result[1,1] := Cos( Value );
Result[1,2] := Sin( Value );
Result[2,1] := - Result[1,2];
Result[2,2] := Result[1,1];
Result[3,3] := 1;
end;
function TrueAnomaly( AOrbit : TOrbit ) : real;
var
r : real;
begin
(**)
r := AOrbit.Classic.ArgLat - AOrbit.Classic.ArgPerigee;
if r < 0 then r := TwoPi + r;
Result := r;
end;
function Parameter ( AOrbit : TOrbit ) : real;
begin
Result := AOrbit.Classic.Rp * ( 1 + AOrbit.Classic.e );
end;
function Radius( AOrbit : TOrbit ) : real;
begin
Result := Parameter( AOrbit ) / ( 1 + AOrbit.Classic.e * Cos( TrueAnomaly(AOrbit) ));
end;
function MeanMotion( AOrbit : TOrbit ) : real;
begin
Result := Sqrt( Mu * Parameter( AOrbit )) / Sqr( Radius( AOrbit ));
end;
function Period( AOrbit : TOrbit ) : real;
begin
Result := TwoPi / MeanMotion( AOrbit );
end;
function VRadial( AOrbit : TOrbit ) : real;
begin
Result := Sqrt( Mu/Parameter(AOrbit) ) * AOrbit.Classic.e * Sin( TrueAnomaly( AOrbit ) );
end;
function SemiMajorAxis( AOrbit : TOrbit ) : real;
var
x : real;
begin
x := Mu / Sqr( MeanMotion( AOrbit ) );
Result := Power ( x, 1/3 );
end;
function EccentrAnomaly ( AOrbit : TOrbit ): real;
var
v, e, r : real;
begin
v := TrueAnomaly( AOrbit );
e := AOrbit.Classic.e;
(**)
r := 2 * ArcTan( Sqrt((1-e)/(1+e)) * Tan( 0.5 * v ));
if r < 0 then r := TwoPi + r ;
Result := r;
end;
function MeanAnomaly( AOrbit : TOrbit ) : real;
var
E, r : real;
begin
E := EccentrAnomaly( AOrbit );
(**)
r := E - AOrbit.Classic.e * Sin( E );
if r < 0 then r := TwoPi + r;
Result := r;
end;
function PerigeePassTime( AOrbit : TOrbit ) : real;
begin
Result := - MeanAnomaly( AOrbit ) / MeanMotion( AOrbit );
end;
procedure ClassToXYZ( var AOrbit : TOrbit );
var
Su, So, Si, Cu, Co, Ci,
Vr_r, rn, r : real;
(**) // Внутренния процедура расчета положения КА в ГСК
procedure _ToGeo ;
var
S, Cl, Sl : real;
t : TDateTime;
M : T3DMatrix;
begin
t := AOrbit.Epoch + AOrbit.Time / 86400 ;
S := SiderialTime( t );
M := MAbs_Geo( S );
AOrbit.Geo.Pos := MxV( M, AOrbit.Cartes.Pos );
AOrbit.Geo.R := Modul( AOrbit.Geo.Pos) ;
r := Sqrt( Sqr( AOrbit.Geo.Pos.X ) + Sqr( AOrbit.Geo.Pos.Y ) );
Sl := AOrbit.Geo.Pos.Y / r ;
Cl := AOrbit.Geo.Pos.X / r ;
AOrbit.Geo.Latitude := ArcSin( AOrbit.Geo.Pos.Z / AOrbit.Geo.R ); // Широта
case Sl > 0 of
true : AOrbit.Geo.Longitude := ArcCos( Cl ); // Долгота
false: AOrbit.Geo.Longitude := - ArcCos( Cl );
end;
end;
begin
r := Radius( AOrbit );
rn := r * MeanMotion( AOrbit );
Vr_r := Vradial( AOrbit) / r;
Su := Sin( AOrbit.Classic.ArgLat);
Cu := Cos( AOrbit.Classic.ArgLat);
So:= Sin( AOrbit.Classic.AscNode);
Co := Cos( AOrbit.Classic.AscNode);
Si := Sin( AOrbit.Classic.Incl);
Ci := Cos( AOrbit.Classic.Incl);
AOrbit.Cartes.Pos.X := r * ( Cu * Co - Su * So * Ci );
AOrbit.Cartes.Pos.Y := r * ( Cu * So + Su * Co * Ci );
AOrbit.Cartes.Pos.Z := r * Su * Si;
AOrbit.Cartes.Vel.X := Vr_r * AOrbit.Cartes.Pos.X
- rn * (Su * Co + Cu * Co * Ci);
AOrbit.Cartes.Vel.Y := Vr_r * AOrbit.Cartes.Pos.Y
- rn * (Su * Co - Cu * Co * Ci);
AOrbit.Cartes.Vel.Z := Vr_r * AOrbit.Cartes.Pos.Z
-rn * Cu * Si;
_ToGeo; // Вычислить координаты в ГСК
end;
function KeplerEquation(const e,M : real):real;
{ Решение уравнения Кеплера в цикле repeat }
const
Eps = 1e-9;
var
Eold,Enew : real;
Stop : boolean;
begin
Eold := M;
repeat
Enew := M + e * Sin(Eold);
Stop := Abs(Eold-Enew)<=Eps;
Eold := Enew;
until Stop;
Result := Enew;
end;
(**) // Служебная функция - вычисляет аргумент широты на заданный момент времени
function GetArgLat( ATime : real; AOrb : TOrbit ) : real;
var
M, E, v, u : real;
n : integer;
begin
// Вычислить среднюю аномалию в момент времени Time
M := MeanMotion( AOrb ) * ( ATime - PerigeePassTime( AOrb ) );
// решаем ур-е Кеплера
E := KeplerEquation( AOrb.Classic.e, M );
// находим истинную аномалию
v := 2 * ArcTan(Sqrt((1-AOrb.Classic.e)/(1+AOrb.Classic.e))*Tan(0.5*E));
//привести v к интервалу 0..2пи
if v < 0 then v := TwoPi + v;
// Вычисляем новый аргумент широты
u := v + AOrb.Classic.ArgPerigee;
// Привести к интервалу 0..2Pi
n := Trunc( u / TwoPi ); // целое число периодов
Result := u - n * TwoPi;
end;
|
25.12.2011, 16:11
Линейный вид
