[mmmm]

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;