Симулятор НАСА (Еще один путь к Земле)

Наткнулся на интересный проект от НАСА - GMAT (General Mission Analysis Tool, http://gmat.gsfc.nasa.gov).

Вкратце - это очень мощный инструмент для точного моделирования полета практически любой сложности в Солнечной системе. Программа с открытыми исходниками, развернутой документацией, скриптовым языком.

Начал осваивать. Ради интереса решил смоделировать в ней тот полет "Кон-Тики", который в Орбитере повторить было нельзя - полет по низкой круговой орбите Луны с учетом масконов. В GMAT для этого есть модель гравитационного поля Луны LP165P с достаточно высоким разрешением 165x165. Маршрут выбрал примерно тот же, что и у Перепелкина - от Центрального залива по азимуту -65 мимо кратеров Коперника и Аристарха.

Получилось в итоге не совсем так, как описано у Пухова. Гравитационное поле Луны сильно неоднородное, добиться нулей на приборах невозможно даже на короткое время. Поспать бы у Перепелкина не получилось ;) Высота начинает падать с самого начала пути. И это несмотря на то, что самые крупные масконы встретились только в середине (возле кратера Гагарина) и в конце (возле Моря Нектара).

В общем, как оказалось, долго на низкой орбите над Луной по баллистической траектории не полетаешь - очень быстро свалишься. Даже на одном витке высота с 4 км. опустилась ниже 2.5 км.

Выкладываю скрипт.

Исходные условия - корабль на круговой орбите высотой 4 км. над Центральным заливом. Для сравнения там же (отображается на графике зеленой линией) его дубликат, траектория которого считается по точечной гравитации (но с учетом возмущений от Солнца и Земли - как в Орбитере).

Create Spacecraft KonTiki;
GMAT KonTiki.DateFormat = TDBGregorian;
GMAT KonTiki.Epoch = '02 Jul 2086 12:00:00.000';
GMAT KonTiki.CoordinateSystem = LunaInt;
GMAT KonTiki.DisplayStateType = SphericalAZFPA;
GMAT KonTiki.RMAG = 1742.2;
GMAT KonTiki.RA = 177.4157916541891;
GMAT KonTiki.DEC = 4.351387326645394;
GMAT KonTiki.VMAG = 1.677540870580422;
GMAT KonTiki.AZI = -64;
GMAT KonTiki.FPA = 90;

Create Spacecraft KonTikiSimple;
GMAT KonTikiSimple.DateFormat = TDBGregorian;
GMAT KonTikiSimple.Epoch = '02 Jul 2086 12:00:00.000';
GMAT KonTikiSimple.CoordinateSystem = LunaInt;
GMAT KonTikiSimple.DisplayStateType = SphericalAZFPA;
GMAT KonTikiSimple.RMAG = 1742.2;
GMAT KonTikiSimple.RA = 177.4157916541891;
GMAT KonTikiSimple.DEC = 4.351387326645394;
GMAT KonTikiSimple.VMAG = 1.677540870580422;
GMAT KonTikiSimple.AZI = -64;
GMAT KonTikiSimple.FPA = 90;

Create ForceModel MoonProp_ForceModel;
GMAT MoonProp_ForceModel.CentralBody = Luna;
GMAT MoonProp_ForceModel.PrimaryBodies = {Luna};
GMAT MoonProp_ForceModel.PointMasses = {Earth, Jupiter, Sun};
GMAT MoonProp_ForceModel.Drag = None;
GMAT MoonProp_ForceModel.SRP = Off;
GMAT MoonProp_ForceModel.RelativisticCorrection = On;
GMAT MoonProp_ForceModel.ErrorControl = RSSStep;
GMAT MoonProp_ForceModel.GravityField.Luna.Degree = 165;
GMAT MoonProp_ForceModel.GravityField.Luna.Order = 165;
GMAT MoonProp_ForceModel.GravityField.Luna.PotentialFile = 'LP165P.cof';

Create ForceModel PointProp_ForceModel;
GMAT PointProp_ForceModel.CentralBody = Luna;
GMAT PointProp_ForceModel.PointMasses = {Earth, Jupiter, Luna, Sun};
GMAT PointProp_ForceModel.Drag = None;
GMAT PointProp_ForceModel.SRP = Off;
GMAT PointProp_ForceModel.RelativisticCorrection = On;
GMAT PointProp_ForceModel.ErrorControl = RSSStep;


Create Propagator MoonProp;
GMAT MoonProp.FM = MoonProp_ForceModel;
GMAT MoonProp.Type = RungeKutta89;
GMAT MoonProp.InitialStepSize = 1;
GMAT MoonProp.Accuracy = 1e-014;
GMAT MoonProp.MinStep = 0.001;
GMAT MoonProp.MaxStep = 100;
GMAT MoonProp.MaxStepAttempts = 50;
GMAT MoonProp.StopIfAccuracyIsViolated = true;

Create Propagator PointProp;
GMAT PointProp.FM = PointProp_ForceModel;
GMAT PointProp.Type = RungeKutta89;
GMAT PointProp.InitialStepSize = 1;
GMAT PointProp.Accuracy = 1e-014;
GMAT PointProp.MinStep = 0.001;
GMAT PointProp.MaxStep = 100;
GMAT PointProp.MaxStepAttempts = 50;
GMAT PointProp.StopIfAccuracyIsViolated = true;


Create CoordinateSystem LunaInt;
GMAT LunaInt.Origin = Luna;
GMAT LunaInt.Axes = BodyInertial;

Create CoordinateSystem LunaFixed;
GMAT LunaFixed.Origin = Luna;
GMAT LunaFixed.Axes = BodyFixed;


Create XYPlot Altitude;
GMAT Altitude.XVariable = KonTiki.ElapsedSecs;
GMAT Altitude.YVariables = {KonTiki.Luna.Altitude, KonTikiSimple.Luna.Altitude};

Create GroundTrackPlot GroundTrack;
GMAT GroundTrack.Add = {KonTiki};
GMAT GroundTrack.CentralBody = Luna;

Create XYPlot VerticalV;
GMAT VerticalV.XVariable = KonTiki.ElapsedSecs;
GMAT VerticalV.YVariables = {VerticalSpeed};

Create XYPlot VerticalAccel;
GMAT VerticalAccel.XVariable = KonTiki.ElapsedSecs;
GMAT VerticalAccel.YVariables = {VerticalACC};

Create Variable VerticalSpeed VerticalSpeed0 ElapsedSeconds0 VerticalACC Altitude0;
GMAT VerticalSpeed = 0;
GMAT VerticalSpeed0 = 0;
GMAT ElapsedSeconds0 = 0;
GMAT VerticalACC = 0;
GMAT Altitude0 = 0;


BeginMissionSequence;

GMAT VerticalSpeed=cos(DegToRad(KonTiki.LunaFixed.FPA))*KonTiki.LunaFixed.VMAG*1000;
GMAT ElapsedSeconds0 = KonTiki.ElapsedSecs;
GMAT Altitude0 = KonTiki.Luna.Altitude;
While KonTiki.ElapsedSecs < 6585
 Propagate Synchronized MoonProp(KonTiki) PointProp(KonTikiSimple);
 GMAT VerticalSpeed = (KonTiki.Luna.Altitude-Altitude0)/(KonTiki.ElapsedSecs-ElapsedSeconds0)*1000;
 GMAT VerticalACC = (VerticalSpeed-VerticalSpeed0)/(KonTiki.ElapsedSecs-ElapsedSeconds0);
 GMAT ElapsedSeconds0 = KonTiki.ElapsedSecs;
 GMAT VerticalSpeed0 = VerticalSpeed;
 GMAT Altitude0 = KonTiki.Luna.Altitude;
EndWhile;