diff --git a/New Model/3DQSO-Lb_5rev_right.mat b/New Model/3DQSO-Lb_5rev_right.mat
new file mode 100644
index 0000000000000000000000000000000000000000..dfc001877144cf66e87bffed85cbe7aca578a9ab
Binary files /dev/null and b/New Model/3DQSO-Lb_5rev_right.mat differ
diff --git a/New Model/3DQSO-Lc_5rev_right.mat b/New Model/3DQSO-Lc_5rev_right.mat
new file mode 100644
index 0000000000000000000000000000000000000000..2e794e2334f06358e2ab5ad74a4037ba44c7a91d
Binary files /dev/null and b/New Model/3DQSO-Lc_5rev_right.mat differ
diff --git a/New Model/BSP_generator.m b/New Model/BSP_generator.m
new file mode 100644
index 0000000000000000000000000000000000000000..6a2434745d4c623b394e743534a6743a46f8e272
--- /dev/null
+++ b/New Model/BSP_generator.m
@@ -0,0 +1,144 @@
+clear
+close all
+clc
+
+%% Definition of the points to start the continuation for the QSO
+
+ warning('off', 'all'); format longG;
+ set(0,'DefaultTextInterpreter','latex');
+ set(0,'DefaultAxesFontSize', 16);
+
+% Mac
+ restoredefaultpath
+ addpath('../../paper/MMX_Fcn_CovarianceAnalyses');
+ addpath('../../paper/MMX_Product/MMX_BSP_Files_GravLib/');
+ addpath('../../paper/Reference_Definition/');
+ addpath('../../paper/');
+ addpath(genpath('../../mice/'));
+ addpath(genpath('../../generic_kernels/'));
+ addpath(genpath('../../useful_functions/'));
+ addpath(genpath('../../dynamical_model/'));
+ MMX_InitializeSPICE
+ cspice_furnsh(which('mar097.bsp'));
+ cspice_furnsh('MARPHOFRZ.txt');
+ cspice_furnsh(which('Phobos_826891269_828619269.bsp'));
+
+
+ alpha = 13.03; %km
+ beta = 11.4; %km
+ gamma = 9.14; %km
+ Geom.alpha = alpha;
+ Geom.beta = beta;
+ Geom.gamma = gamma;
+ muPh = 7.1139e-4; %km^3/s^2
+ muM = 4.2828e4; %km^3/s^2
+ pPh = 9377.2; %km
+ omega = sqrt((muM + muPh)/pPh^3); %rad/s
+ T = 1/omega;
+ load("Tpo.mat")
+ [pars, units] = MMX_DefineNewModelParametersAndUnits_Right;
+
+% Initial MMX's state vector
+ startdate = '2026-03-16 00:00:00 (UTC)';
+ data = cspice_str2et(startdate);
+ pars.et0 = data;
+ [Ph, pars] = Phobos_States_NewModel(data,pars);
+
+
+ % QSO-Lb
+ load("SwingQSO-Lc_5rev_right.mat")
+
+% If MMX comes from the optimization
+ X_MMX = X_good(1:6,1);
+
+% If MMX comes from a .bsp file
+ % cspice_furnsh(which('MMX_QSO_031_2x2_826891269_828619269.bsp'));
+ % X_MMX = cspice_spkezr('-34', data, 'MARSIAU', 'none', '499')./units.sfVec;
+
+
+% Initial Phobos's state vector
+ Ph0 = Ph./units.sfVec2;
+ Iz = pars.IPhz_bar + pars.ni*Ph0(1)^2;
+ K = Iz*Ph0(4) + pars.IPhz_bar*Ph0(8);
+ pars.K = K;
+
+% Initial state vector
+ St0 = [X_MMX; Ph0(1:2); Ph0(7:8); pars.IPhx_bar; pars.IPhy_bar; pars.IPhz_bar];
+
+% Integration
+ day = 86400;
+ period_Def = 30*day;
+ tspan = [0, period_Def]*units.tsf;
+ RelTol = 1e-13;
+ AbsTol = 1e-16;
+ opt = odeset('RelTol',RelTol,'AbsTol',AbsTol);
+ [t,X] = ode113(@(t,X) Dynamics_MPHandMMX_Relative2(t,X,pars,units),tspan,St0,opt);
+ Xmmx = [pars.PB', zeros(3,3); zeros(3,3), pars.PB']*X(:,1:6)'.*units.sfVec;
+ t = t/units.tsf;
+
+
+
+%% BSP Generation
+
+ Model.epoch = data;
+ Model.pars = pars;
+ Model.units = units;
+ micePATH = '~/Documents/MATLAB/mice/';
+
+ TargetFolder = './MMX_BSP_Files_Relative/';
+ if (exist(TargetFolder, 'dir') ~= 7)
+ command = sprintf('mkdir %s', TargetFolder);
+ system(command);
+ end
+
+ MMX_GenerateBSPFiles_Relative(t, Xmmx, 'SwingQSO', [27, 12], period_Def, Model, micePATH)
+
+
+%% Test risultati
+
+ addpath(genpath('./MMX_BSP_Files_Relative/'))
+ cspice_furnsh(which('MMX_SwingQSO_027_012_2x2_826891269_829483269_relative_right.bsp'));
+ MMX_BSP = cspice_spkezr('-35', data+t', 'IAU_Phobos', 'none', '-401');
+ err_MMX = Xmmx - MMX_BSP;
+
+ figure(1)
+ subplot(2,3,1)
+ plot(t/3600,err_MMX(1,:)*units.lsf)
+ ylabel('$\Delta X [km]$',Interpreter='latex')
+ grid on
+ subplot(2,3,2)
+ plot(t/3600,err_MMX(2,:)*units.lsf)
+ ylabel('$\Delta Y [km]$',Interpreter='latex')
+ grid on
+ subplot(2,3,3)
+ plot(t/3600,err_MMX(3,:)*units.lsf)
+ ylabel('$\Delta Z [km]$',Interpreter='latex')
+ grid on
+ subplot(2,3,4)
+ plot(t/3600,err_MMX(4,:)*units.vsf)
+ ylabel('$\Delta Vx [km/s]$',Interpreter='latex')
+ grid on
+ grid on
+ subplot(2,3,5)
+ plot(t/3600,err_MMX(5,:)*units.vsf)
+ ylabel('$\Delta Vy [km/s]$',Interpreter='latex')
+ grid on
+ subplot(2,3,6)
+ plot(t/3600,err_MMX(6,:)*units.vsf)
+ ylabel('$\Delta Vz [km/s]$',Interpreter='latex')
+ grid on
+
+
+ figure(2)
+ plot3(Xmmx(1,:), Xmmx(2,:), Xmmx(3,:));
+ hold on;
+ grid on;
+ plot3(MMX_BSP(1,:), MMX_BSP(2,:), MMX_BSP(3,:));
+ planetPlot('Asteroid',[0;0;0],Geom,1);
+ xlabel('x [km]');
+ ylabel('y [km]');
+ zlabel('z [km]');
+ axis equal
+
+
+
\ No newline at end of file
diff --git a/New Model/Constraints.m b/New Model/Constraints.m
new file mode 100644
index 0000000000000000000000000000000000000000..3fc9ded5dc7bc0b65f976918b84d4f0f47767221
--- /dev/null
+++ b/New Model/Constraints.m
@@ -0,0 +1,74 @@
+function [c, ceq] = Constraints(X,pars,units,X_patch)
+
+
+ [Ph, pars] = Phobos_States_NewModel(pars.et0,pars);
+ Ph = Ph./units.sfVec2;
+
+ RelTol = 1e-13;
+ AbsTol = 1e-16;
+ St0 = [Ph(1:2); Ph(7:8)];
+ tspanp = X(end,:);
+ opt = odeset('RelTol',RelTol,'AbsTol',AbsTol);
+ [~,X_Ph0] = ode113(@(t,X) MPH_Relative(t,X,pars),tspanp,St0,opt);
+ X_Ph0 = X_Ph0(:,1:4)';
+
+% Non-linear constraint imposing that the two propagations of the two
+% branches match
+
+% Init dei constraints
+ DeltaX = zeros(6*(size(X,2)-2),1);
+ Perp = zeros((size(X,2))-2,1);
+ times = zeros((size(X,2)-1),1);
+ boxPrimo_ultimo = zeros(6,1);
+ box = zeros(3*size(X,2),1);
+
+ opt = odeset('RelTol',RelTol,'AbsTol',AbsTol,'event',@(t,X) landing_Phobos(t,X,pars,units));
+
+
+ for i = 2:size(X,2)
+
+ St1meno = [X(1:6,i-1); X_Ph0(1:end,i-1); pars.IPhx_bar; pars.IPhy_bar; pars.IPhz_bar];
+ tspan1meno = [0, X(end-1,i-1)];
+ [~,X1meno] = ode113(@(t,X) Dynamics_MPHandMMX_Relative2(t,X,pars,units),tspan1meno,St1meno,opt);
+
+ St1plus = [X(1:6,i); X_Ph0(1:end,i); pars.IPhx_bar; pars.IPhy_bar; pars.IPhz_bar];
+ tspan1plus = [0, X(end-2,i)];
+ [~,X1plus] = ode113(@(t,X) -Dynamics_MPHandMMX_Relative2(t,X,pars,units),tspan1plus,St1plus,opt);
+
+ DeltaX(1+(i-2)*6:(6+(i-2)*6)) = (X1meno(end,1:6)' - X1plus(end,1:6)');
+
+ end
+
+
+% Constraint imposing that position vector and velocity vector are
+% perpendicular at patching points
+
+
+ for i = 1:size(X,2)
+
+ Perp(i) = dot(X(1:3,i), X(4:6,i));
+
+ if i==1
+ boxPrimo_ultimo(1:3) = abs(X_patch(:,i) - X(1:3,i)*units.lsf) - .1;
+ elseif i==size(X,2)
+ boxPrimo_ultimo(4:6) = abs(X_patch(:,i) - X(1:3,i)*units.lsf) - 1;
+ end
+
+ box(1+(i-1)*3:3+(i-1)*3) = abs(X_patch(:,i) - X(1:3,i)*units.lsf) - 1;
+
+ end
+
+
+% Constraint imposing that forward and backward propagation times are
+% meeting somewhere between half orbit revolution time
+
+ for i = 2:size(X,2)
+
+ times(i-1) = (X(end,i) - X(end,i-1) - (X(end-2,i) + X(end-1,i-1)));
+
+ end
+
+ c = [];
+ ceq = [DeltaX; times];
+
+end
\ No newline at end of file
diff --git a/New Model/CostFunction.m b/New Model/CostFunction.m
new file mode 100644
index 0000000000000000000000000000000000000000..eb6d4560fcfc16fd2f6396d4b6f2e4eeb59c88c2
--- /dev/null
+++ b/New Model/CostFunction.m
@@ -0,0 +1,35 @@
+function DeltaV = CostFunction(X,pars,units)
+
+
+ [Ph, pars] = Phobos_States_NewModel(pars.et0, pars);
+ Ph = Ph./units.sfVec2;
+
+ RelTol = 1e-13;
+ AbsTol = 1e-16;
+ St0 = [Ph(1:2); Ph(7:8)];
+ tspanp = X(end,:);
+ opt = odeset('RelTol',RelTol,'AbsTol',AbsTol);
+ [~,X_Ph0] = ode113(@(t,X) MPH_Relative(t,X,pars),tspanp,St0,opt);
+ X_Ph0 = X_Ph0(:,1:4)';
+
+% L'ottimizzazione minimizza il DeltaV richiesto nei punti in cui le
+% propagazioni si incontrano per far incontrare le traiettorie
+
+ DeltaV = 0;
+ opt = odeset('RelTol',1e-13,'AbsTol',1e-16,'event',@(t,X) landing_Phobos(t,X,pars,units));
+
+ for i = 2:size(X,2)
+
+ St1meno = [X(1:6,i-1); X_Ph0(1:end,i-1); pars.IPhx_bar; pars.IPhy_bar; pars.IPhz_bar];
+ tspan1meno = [0, X(end-1,i-1)];
+ [~,X1meno] = ode113(@(t,X) Dynamics_MPHandMMX_Relative2(t,X,pars,units),tspan1meno,St1meno,opt);
+
+ St1plus = [X(1:6,i); X_Ph0(1:end,i); pars.IPhx_bar; pars.IPhy_bar; pars.IPhz_bar];
+ tspan1plus = [0, X(end-2,i)];
+ [~,X1plus] = ode113(@(t,X) -Dynamics_MPHandMMX_Relative2(t,X,pars,units),tspan1plus,St1plus,opt);
+
+ DeltaV = DeltaV + norm(X1meno(end,4:6)' - X1plus(end,4:6)')*units.vsf;
+
+ end
+
+end
\ No newline at end of file
diff --git a/New Model/Dynamics_MPHandMMX_Relative.m b/New Model/Dynamics_MPHandMMX_Relative.m
index 6754cc76c8cdd3b6a9edd9c8792f4facdccffc7c..815204446646d23346d3d1a2218a3732072bec45 100644
--- a/New Model/Dynamics_MPHandMMX_Relative.m
+++ b/New Model/Dynamics_MPHandMMX_Relative.m
@@ -56,12 +56,7 @@ function dx = Dynamics_MPHandMMX_Relative(~,x,pars,~)
IPhy_bar = x(16);
IPhz_bar = x(17);
-% Position of Mars wrt Phobos
- rPh = RPh*[cos(Xsi); sin(Xsi); 0];
-
-% Rotation matrix from Phobos fixed to Mars pointing reference coordinates frame
- PhM = [cos(Xsi), sin(Xsi), 0; -sin(Xsi), cos(Xsi), 0; 0, 0, 1];
- MPh = PhM';
+
%% Phobos's orbit
@@ -75,18 +70,22 @@ function dx = Dynamics_MPHandMMX_Relative(~,x,pars,~)
RPh_ddot = theta_dot^2*RPh - dVdRPh/pars.ni;
theta_ddot = dVdXsi/(pars.ni*RPh^2) - 2*RPh_dot*theta_dot/RPh;
Phi_ddot = - dVdXsi/(pars.ni*RPh^2) + 2*RPh_dot*theta_dot/RPh;
- Xsi_ddot = -(1 + pars.ni*RPh^2/IPhz_bar)*dVdXsi/(pars.ni*RPh^2) +...
+ Xsi_ddot = - (1 + pars.ni*RPh^2/IPhz_bar)*dVdXsi/(pars.ni*RPh^2) +...
+ 2*RPh_dot*theta_dot/RPh;
%% MMX 2BP with Mars + 3rd body Phobos
+
+% Position of Mars wrt Phobos
+ rPh = RPh*[cos(-Xsi); sin(-Xsi); 0];
+
% The gravitational effect of Mars on MMX is 2BP+J2
- r = rPh - rsb;
+ r = rsb - rPh;
R = norm(r);
- gM1 = pars.G*pars.MMars*r/R^3;
- gJ2M = 1.5*pars.G*pars.MMars*pars.J2*pars.RMmean/R^5*...
+ gM1 = - pars.G*pars.MMars*r/R^3;
+ gJ2M = - 1.5*pars.G*pars.MMars*pars.J2*pars.RMmean/R^5*...
[(1-5*(r(3)/R)^2)*r(1); (1-5*(r(3)/R)^2)*r(2);...
(3-5*(r(3)/R)^2)*r(3)];
@@ -95,23 +94,25 @@ function dx = Dynamics_MPHandMMX_Relative(~,x,pars,~)
% The gravitational effect of Phobos on MMX includes the Ellipsoid
% harmonics
- C_22_Ph = .25*(IPhy_bar - IPhx_bar);
- C_20_Ph = -.5*(2*IPhz_bar - IPhx_bar - IPhy_bar);
-
- pars.SH.Clm = [1, 0, 0; 0, 0, 0; C_20_Ph, 0, C_22_Ph]./pars.SH.Norm(1:3,1:3);
- [gPh,~,~,~] = CGM(rsb, pars.SH);
+ % C_22_Ph = .25*(IPhy_bar - IPhx_bar);
+ % C_20_Ph = -.5*(2*IPhz_bar - IPhx_bar - IPhy_bar);
+ %
+ % pars.SH.Clm = [1, 0, 0; 0, 0, 0; C_20_Ph, 0, C_22_Ph]./pars.SH.Norm(1:3,1:3);
+ [gPh,~,~,~] = CGM(rsb, pars.SH);
-% Combined effect on MMX
- W = [0; 0; theta_dot-Xsi_dot];
- W_dot = [0; 0; theta_ddot-Xsi_ddot];
+% MMX motion
+ W = [0; 0; theta_dot + Xsi_dot];
+ W1 = [0; 0; theta_dot - Xsi_dot];
+ W_dot = [0; 0; theta_ddot + Xsi_ddot];
- g_tot = gPh - cross(W_dot,rsb) - 2 * cross(W,vsb) - cross(W,cross(W,rsb)) + gM;
+ r_dot = vsb;
+ v_dot = gPh - cross(W_dot,rsb) - 2 * cross(W,vsb) - cross(W,cross(W,rsb)) + gM;
%% Output for the integrator
- dx = [vsb; g_tot; RPh_dot; RPh_ddot; theta_dot; theta_ddot;...
+ dx = [r_dot; v_dot; RPh_dot; RPh_ddot; theta_dot; theta_ddot;...
Phi_dot; Phi_ddot; Xsi_dot; Xsi_ddot; 0; 0; 0];
diff --git a/New Model/Dynamics_MPHandMMX_Relative2.m b/New Model/Dynamics_MPHandMMX_Relative2.m
index 6fa278a1d3fc42c26357d2c0fd4c63e4603d733d..da475f73d162007099359598bf054548613c6bf7 100644
--- a/New Model/Dynamics_MPHandMMX_Relative2.m
+++ b/New Model/Dynamics_MPHandMMX_Relative2.m
@@ -53,11 +53,8 @@ function dx = Dynamics_MPHandMMX_Relative2(~,x,pars,~)
IPhz_bar = x(13);
% Position of Mars wrt Phobos
- rPh = RPh*[cos(Xsi); sin(Xsi); 0];
+ rPh = RPh*[cos(-Xsi); sin(-Xsi); 0];
-% Rotation matrix from Phobos fixed to Mars pointing reference coordinates frame
- PhM = [cos(Xsi), sin(Xsi), 0; -sin(Xsi), cos(Xsi), 0; 0, 0, 1];
- MPh = PhM';
%% Phobos's orbit
@@ -78,30 +75,39 @@ function dx = Dynamics_MPHandMMX_Relative2(~,x,pars,~)
% The gravitational effect of Mars on MMX is 2BP+J2
+ % r = rsb - rPh;
+ % R = norm(r);
+ % gM1 = - pars.G*pars.MMars*r/R^3;
+ % gJ2M = - 1.5*pars.G*pars.MMars*pars.J2*pars.RMmean/R^5*...
+ % [(1-5*(r(3)/R)^2)*r(1); (1-5*(r(3)/R)^2)*r(2);...
+ % (3-5*(r(3)/R)^2)*r(3)];
+ %
+ % gM = gM1 + gJ2M - pars.G*pars.MMars*rPh/norm(rPh)^3;
+
r = rPh - rsb;
R = norm(r);
gM1 = pars.G*pars.MMars*r/R^3;
gJ2M = 1.5*pars.G*pars.MMars*pars.J2*pars.RMmean/R^5*...
[(1-5*(r(3)/R)^2)*r(1); (1-5*(r(3)/R)^2)*r(2);...
(3-5*(r(3)/R)^2)*r(3)];
-
- gM = gM1 + gJ2M - pars.G*pars.MMars*rPh/norm(rPh)^3;
+
+ gT = pars.G*pars.MMars*rPh/RPh^3;
+
+ gM = gM1 + gJ2M - gT;
% The gravitational effect of Phobos on MMX includes the Ellipsoid
% harmonics
- C_22_Ph = .25*(IPhy_bar - IPhx_bar);
- C_20_Ph = -.5*(2*IPhz_bar - IPhx_bar - IPhy_bar);
-
- pars.SH.Clm = [1, 0, 0; 0, 0, 0; C_20_Ph, 0, C_22_Ph]./pars.SH.Norm(1:3,1:3);
[gPh,~,~,~] = CGM(rsb, pars.SH);
% Combined effect on MMX
theta_dot = (pars.K - IPhz_bar*Xsi_dot)/Iz;
- W = [0; 0; theta_dot-Xsi_dot];
+ theta_ddot = dVdXsi/(pars.ni*RPh^2) - 2*RPh_dot*theta_dot/RPh;
+ W = [0; 0; theta_dot + Xsi_dot];
+ W_dot = [0; 0; theta_ddot + Xsi_ddot];
- g_tot = gPh - 2 * cross(W,vsb) - cross(W,cross(W,rsb)) + gM;
+ g_tot = gPh + gM - cross(W_dot,rsb) - 2 * cross(W,vsb) - cross(W,cross(W,rsb));
%% Output for the integrator
diff --git a/New Model/MMX_DefineNewModelParametersAndUnits_Right.m b/New Model/MMX_DefineNewModelParametersAndUnits_Right.m
new file mode 100644
index 0000000000000000000000000000000000000000..c9ba3bbb25a7d643058008273fe4ac46a797825e
--- /dev/null
+++ b/New Model/MMX_DefineNewModelParametersAndUnits_Right.m
@@ -0,0 +1,108 @@
+function [pars, units] = MMX_DefineNewModelParametersAndUnits_Right
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% [pars, units] = MMX_DefineModelParametersAndUnits(PhobosGravityField)
+%
+% Define model parameters and units for the MMX Toolbox
+%
+% Inputs:
+% PhobosGravityField Order x Degree of Phobos' gravity field (max [8, 8])
+%
+% Outputs:
+% pars.d No. of States
+% pars.p No. of model parameters
+% pars.Body Secondary's body structure:
+% .alp Secondary's largest semi-axis [Normalized Units]
+% .bet Secondary's intermediate semi-axis [Normalized Units]
+% .gam Secondary's smallest semi-axis [Normalized Units]
+%
+% units.lsf Length unit
+% units.tsf Time unit
+% units.vsf Velocity unit
+% units.sfVec Vector of normalizing units [lsf, lsf, lsf, vsf, vsf, vsf]
+%
+% Author: STAG Team
+% Date: Jun 2021
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+
+
+% Parameters
+ G = astroConstants(1);
+ pars.GM2 = cspice_bodvrd('401','GM',1); %[kg]
+ pars.GM1 = cspice_bodvrd('499','GM',1); %[kg]
+ r_Ph = cspice_bodvrd('401','RADII',3); %[km]
+ alpha = r_Ph(1);
+ r_Ph = r_Ph/alpha;
+ a_Ph = r_Ph(1);
+ b_Ph = r_Ph(2); %[km]
+ c_Ph = r_Ph(3); %[km]
+ pars.RPhmean = sqrt(3/(a_Ph+b_Ph+c_Ph));
+ r_M = cspice_bodvrd('499','RADII',3); %[km]
+ r_M = r_M/alpha;
+ a_M = r_M(1); %[km]
+ b_M = r_M(2); %[km]
+ c_M = r_M(3); %[km]
+
+ units.lsf = alpha; % Unit of lenght [km]
+ units.tsf = sqrt((pars.GM1 + pars.GM2)/units.lsf^3); % Unit of time [1/s]
+ units.vsf = units.lsf*units.tsf;
+ units.GMsf = units.lsf^3*units.tsf^2;
+
+ units.sfVec = [units.lsf; units.lsf; units.lsf; units.vsf; units.vsf; units.vsf];
+ units.sfVec2 = [units.lsf; units.vsf; 1; units.tsf; 1; units.tsf; 1; units.tsf];
+ units.M = pars.GM2/astroConstants(1); % Unit of mass [km]
+
+% Bodies Geometry
+ pars.Phobos.alpha = a_Ph*units.lsf;
+ pars.Phobos.beta = b_Ph*units.lsf;
+ pars.Phobos.gamma = c_Ph*units.lsf;
+
+ pars.Mars.alpha = a_M*units.lsf;
+ pars.Mars.beta = b_M*units.lsf;
+ pars.Mars.gamma = c_M*units.lsf;
+
+% Mass fractions
+ pars.m = pars.GM2^2/(pars.GM1 + pars.GM2);
+ pars.mu = (pars.GM1 + pars.GM2);
+ pars.ni = pars.GM1/(pars.GM1 + pars.GM2);
+
+ pars.G = G/(units.lsf^3*units.tsf^2/units.M);
+ pars.MMars = pars.GM1/(G*units.M);
+ pars.MPho = pars.GM2/(G*units.M);
+
+% Inertia parameters
+ pars.IPhx_bar = .2*(b_Ph^2 + c_Ph^2);
+ pars.IPhy_bar = .2*(a_Ph^2 + c_Ph^2);
+ pars.IPhz_bar = .2*(a_Ph^2 + b_Ph^2);
+
+ pars.IMx_bar = .2*pars.GM1*(b_M^2 + c_M^2)/(a_Ph^2*pars.GM1);
+ pars.IMy_bar = pars.IMx_bar;
+ pars.Is = pars.IMx_bar;
+ pars.IMz_bar = .2*pars.GM1*(a_M^2 + b_M^2)/(a_Ph^2*pars.GM1);
+
+% Phobos gravitational harmonics WRT principal axis of inertia
+
+ pars.SH = PhobosSphericalHarmonicCoefficients(2,2);
+ pars.SH.Re = pars.SH.Re/units.lsf;
+ pars.SH.GM = pars.SH.GM/units.GMsf;
+
+% par.MMars gravitational harmonics WRT principal axis of inertia
+ pars.RMmean = mean(r_M);
+ C_20_M = (pars.IMx_bar - pars.IMz_bar)/(pars.MMars*pars.RMmean^2);
+ pars.J2 = - C_20_M;
+
+% Some other useful parameters
+ pars.d = 6; % MMX's States
+ pars.d2 = 8; % Phobos's States
+
+% Ellipsoidal Phobos gravitational harmonics WRT principal axis of inertia
+ C_20_Ph = .25*(pars.IPhy_bar - pars.IPhx_bar);
+ C_22_Ph = -.5*(2*pars.IPhz_bar - pars.IPhx_bar - pars.IPhy_bar);
+
+ pars.SH.Clm = [1, 0, 0; 0, 0, 0; C_20_Ph, 0, C_22_Ph]./pars.SH.Norm(1:3,1:3);
+ pars.SH.Slm = zeros(3,3);
+
+
+end
\ No newline at end of file
diff --git a/New Model/MMX_GenerateBSPFiles_Relative.m b/New Model/MMX_GenerateBSPFiles_Relative.m
new file mode 100644
index 0000000000000000000000000000000000000000..1240fd7178cc8b83b36f027a2ddddec7ea19328a
--- /dev/null
+++ b/New Model/MMX_GenerateBSPFiles_Relative.m
@@ -0,0 +1,31 @@
+function MMX_GenerateBSPFiles_Relative(t, Xmmx, OrbitType, Amplitudes, PropTime, Model, micePATH)
+
+
+
+%% Load SPICE
+ MMX_InitializeSPICE;
+ et0 = Model.epoch;
+ et1 = et0 + PropTime;
+
+
+
+%% Create MMX Ephemeris file
+ if(strcmp(OrbitType, 'QSO'))
+ amp = sprintf('%03.0f', Amplitudes(1));
+ else
+ amp = sprintf('%03.0f_%03.0f', Amplitudes(1), Amplitudes(2));
+ end
+ dates = sprintf('%09.0f_%09.0f', et0, et1);
+ SHmodel = sprintf('%dx%d',size(Model.pars.SH.Clm,1)-1,size(Model.pars.SH.Clm,2)-1);
+ bspfilename = ['MMX_',OrbitType,'_',amp,'_',SHmodel,'_',dates,'_relative_right.bsp'];
+
+ if(exist(bspfilename,'file')~=2)
+ ID = -35;
+ MMX_MakeBSPFile_GravLib(bspfilename, ID, t, Xmmx, '-401', 'IAU_PHOBOS', Model, [micePATH, '/exe/']);
+ movefile(bspfilename, './MMX_BSP_Files_Relative');
+ else
+ fprintf('MMX ephemeris file already exists!\n');
+ end
+
+
+end
diff --git a/New Model/MPH_Relative.m b/New Model/MPH_Relative.m
new file mode 100644
index 0000000000000000000000000000000000000000..8d2346992d7e78a4ee198881069bfaac5600f865
--- /dev/null
+++ b/New Model/MPH_Relative.m
@@ -0,0 +1,64 @@
+function dx = MPH_Relative(~,x,pars,~)
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% dx = MPH_Relative(t,x,parss,units)
+%
+% Calculate the acceleration vector in the N-Body Problem with Mars and
+% Phobos only and the derivative of the coefficients to be estimated
+%
+% Input:
+% .t Elapsed time since epoch (-)
+% .x State Vector
+% .parss Problem parsameters
+% .d Dimension of the Problem
+% .et0 Initial epoch (-)
+% .refcenter Origin of Frame (e.g., '399')
+% .refframe Coordinate Frame (e.g., 'ECLIPJ2000' or 'J2000')
+% .Ntb No. of third bodies
+% .BodiesID Bodies' SPICE ID (e.g., '10','399', ...)
+% .GM Bodies' Gravitational parsameters (-)
+% .Nsh No. of Spherical Harmonic Bodies
+% .SH_BodiesID Spherical Harmonic Bodies' SPICE ID (e.g., '10','399', ...)
+% .SH Spherical Harmonic Structure
+% .lsf Length scale factor (km)
+% .tsf Time scale factor (s)
+% .STM true or false depending if you calculate STM or not
+% .nCoeff Number if coefficients in the state vector
+%
+% Output:
+% .dx Vector containing accelerations
+%
+% Author: E.Ciccarelli
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+
+%% Phobos's orbit
+
+ RPh = x(1);
+ RPh_dot = x(2);
+ Xsi = x(3);
+ Xsi_dot = x(4);
+
+%% Phobos's orbit
+
+% Potential's first order parstials
+ dVdRPh = pars.ni/RPh^2*(1 + 3/(2*RPh^2)*((pars.IMz_bar - pars.Is) -...
+ .5*pars.IPhx_bar - .5*pars.IPhy_bar + pars.IPhz_bar + ...
+ 1.5*(pars.IPhy_bar - pars.IPhx_bar)*cos(2*Xsi)));
+ dVdXsi = 1.5*pars.ni/RPh^3 * (pars.IPhy_bar - pars.IPhx_bar)*sin(2*Xsi);
+ Iz = pars.IPhz_bar + pars.ni*RPh^2;
+
+% Phobos equations of motions
+ RPh_ddot = (pars.K - pars.IPhz_bar*Xsi_dot)^2*RPh/Iz^2 - dVdRPh/pars.ni;
+ Xsi_ddot = -(1 + pars.ni*RPh^2/pars.IPhz_bar)*dVdXsi/(pars.ni*RPh^2) +...
+ + 2*(pars.K - pars.IPhz_bar*Xsi_dot)*RPh_dot/(RPh*Iz);
+
+
+
+%% Output for the integrator
+
+ dx = [RPh_dot; RPh_ddot; Xsi_dot; Xsi_ddot];
+
+
+end
\ No newline at end of file
diff --git a/New Model/Plot_InitialCondition.m b/New Model/Plot_InitialCondition.m
new file mode 100644
index 0000000000000000000000000000000000000000..719d0be2adedf34bbc2a2cd44871ba623d9e4b33
--- /dev/null
+++ b/New Model/Plot_InitialCondition.m
@@ -0,0 +1,87 @@
+function stop = Plot_InitialCondition(x,~,~,pars,units)
+
+ stop = false;
+
+ [Ph, pars] = RealPhobos_States_NewModel(pars.et0,pars);
+ Ph = Ph./units.sfVec2;
+ Geom = pars.Geom;
+
+ RelTol = 1e-13;
+ AbsTol = 1e-16;
+ St0 = [Ph(1:2); Ph(7:8)];
+ tspanp = x(end,:);
+ opt = odeset('RelTol',RelTol,'AbsTol',AbsTol);
+ [~,X_Ph0] = ode113(@(t,X) MPH_Relative(t,X,pars),tspanp,St0,opt);
+ X_Ph0 = X_Ph0(:,1:4)';
+
+
+ for j = 1:size(x,2)
+
+ X_0opt = [x(1:6,j); X_Ph0(1:4,j); pars.IPhx_bar; pars.IPhy_bar; pars.IPhz_bar];
+ tspanp = [0,x(8,j)];
+ tspanm = [0,x(7,j)];
+ opt = odeset('RelTol',RelTol,'AbsTol',AbsTol,'event',@(t,X) landing_Phobos(t,X,pars,units));
+ [~,Xp] = ode113(@(t,X) Dynamics_MPHandMMX_Relative2(t,X,pars,units),tspanp,X_0opt,opt);
+ [~,Xm] = ode113(@(t,X) -Dynamics_MPHandMMX_Relative2(t,X,pars,units),tspanm,X_0opt,opt);
+ XPhp = Xp(:,7:10);
+ XPhm = Xm(:,7:10);
+
+
+ X_PhobosMARSIAU = zeros(3,size(Xp,1)+size(Xm,1));
+ X_QSO = zeros(3,size(Xp,1)+size(Xm,1));
+
+ % Braccio in avanti dal patch point
+
+ for i = 1:size(Xp,1)
+
+ X_QSO(:,i+size(Xm,1)) = Xp(i,1:3);
+
+ end
+
+ % Braccio indietro dal patch point
+
+ for i = 1:size(Xm,1)
+
+ X_QSO(:,size(Xm,1)-i+1) = Xm(i,1:3);
+
+ end
+
+
+ X_QSO_patch = zeros(3,size(x,2));
+
+ for i = 1:size(x,2)
+
+ X_QSO_patch(:,i) = x(1:3,i);
+
+ end
+
+% figure(11)
+% plot3(Xp(:,1)*units.lsf,Xp(:,2)*units.lsf,Xp(:,3)*units.lsf);
+% hold on;
+% grid on;
+% axis equal
+% plot3(Xm(:,1)*units.lsf,Xm(:,2)*units.lsf,Xm(:,3)*units.lsf);
+% plot3(x(1,j)*units.lsf,x(2,j)*units.lsf,x(3,j)*units.lsf,'r*');
+% % plot3(X_PhobosMARSIAU(1,:),X_PhobosMARSIAU(2,:),X_PhobosMARSIAU(3,:))
+% planetPlot('Mars',[0;0;0],1);
+% xlabel('x [km]');
+% ylabel('y [km]');
+% zlabel('z [km]');
+% legend('MMX','Patch points','Phobos');
+ figure(10)
+ plot3(X_QSO(1,:)*units.lsf,X_QSO(2,:)*units.lsf,X_QSO(3,:)*units.lsf);
+ hold on;
+ grid on;
+ plot3(X_QSO_patch(1,j)*units.lsf,X_QSO_patch(2,j)*units.lsf,X_QSO_patch(3,j)*units.lsf,'*');
+ planetPlot('Asteroid',[0;0;0],Geom,1);
+ xlabel('x [km]');
+ ylabel('y [km]');
+ zlabel('z [km]');
+ axis equal
+
+
+ end
+
+ hold off
+
+end
\ No newline at end of file
diff --git a/New Model/QSO-Lb_5rev.mat b/New Model/QSO-Lb_5rev.mat
new file mode 100644
index 0000000000000000000000000000000000000000..3478d94397e6d1fca576b38de240f643421a0abb
Binary files /dev/null and b/New Model/QSO-Lb_5rev.mat differ
diff --git a/New Model/QSO-Lb_5rev_right.mat b/New Model/QSO-Lb_5rev_right.mat
new file mode 100644
index 0000000000000000000000000000000000000000..1048dddaa249a796d86886ecc131dd335d58f035
Binary files /dev/null and b/New Model/QSO-Lb_5rev_right.mat differ
diff --git a/New Model/QSO-Lc_5rev_right.mat b/New Model/QSO-Lc_5rev_right.mat
new file mode 100644
index 0000000000000000000000000000000000000000..7fb29cf4a6b0d7dfccc30fd7876794fa6408c2f1
Binary files /dev/null and b/New Model/QSO-Lc_5rev_right.mat differ
diff --git a/New Model/QSOs_Definition.m b/New Model/QSOs_Definition.m
new file mode 100644
index 0000000000000000000000000000000000000000..11f0c24282b1f94bb5795b43d1041d392f77cc68
--- /dev/null
+++ b/New Model/QSOs_Definition.m
@@ -0,0 +1,322 @@
+clear
+close all
+clc
+
+%% Definition of the points to start the continuation for the QSO
+
+ warning('off', 'all'); format longG;
+ set(0,'DefaultTextInterpreter','latex');
+ set(0,'DefaultAxesFontSize', 16);
+
+% Mac
+ addpath('../../paper/');
+ addpath('../../paper/MMX_Fcn_Miscellaneous/');
+ addpath('../../paper/MMX_Product/');
+ addpath('../../paper/MMX_Fcn_CovarianceAnalyses/');
+ addpath('../../paper/Reference_Definition/');
+ addpath(genpath('../../mice/'));
+ addpath(genpath('../../generic_kernels/'));
+ addpath(genpath('../../useful_functions/'));
+ addpath(genpath('../../dynamical_model/'));
+ MMX_InitializeSPICE
+ cspice_furnsh(which('MARPHOFRZ.txt'));
+ cspice_furnsh(which('mar097.bsp'));
+
+%% Patch points
+
+ alpha = 13.03; %km
+ beta = 11.4; %km
+ gamma = 9.14; %km
+ Geom.alpha = alpha;
+ Geom.beta = beta;
+ Geom.gamma = gamma;
+ muPh = 7.1139e-4; %km^3/s^2
+ muM = 4.2828e4; %km^3/s^2
+ pPh = 9377.2; %km
+ omega = sqrt((muM + muPh)/pPh^3); %rad/s
+ eps = (muPh/muM)^(1/3);
+
+ L = eps*(sqrt(muM)/omega)^(2/3);
+ T = 1/omega;
+
+% Adimensionalized
+
+
+ Phobos = struct('alpha',alpha/L,'beta',beta/L,'gamma',gamma/L,...
+ 'mu',muPh/L^3*T^2,'omega',omega*T);
+
+%% Da runnare una sola volta
+
+% r0 = [100;0;0]/L; % km
+% V0 = [0;-2*r0(1);0]; % km/s
+%
+% par.ds = 1e-3; % user-defined step-length parameter
+% par.dsMax = 1e-1; % maximum value for step-length parameter...see line 94
+% par.Iter = 15; % Max. no. of Newton's Method iterations
+% par.Nmax = 170; % No. of desired orbits
+% par.Opt = 5; % Optimal no. of Newton's Iteration - Useful for adaptive step-length...see line 94
+% par.Tol = 1e-10; % Newton's method tolerance
+%
+% X0 = [r0;V0];
+% T0 = 2*pi;
+% dX0_tilde_ds = [-1;0;0;0;2;0]/sqrt(5);
+% dT_tilde_ds = 0;
+% dlam_tilde_ds = 0;
+%
+% [Xpo,Tpo,Lampo,p,q,Lam,bif] = Predictor_Corrector(X0,T0,@HillProb_Sphere_Nd,@landing_Nd,dX0_tilde_ds,dT_tilde_ds,dlam_tilde_ds,Phobos,par);
+% save('Xpo','Xpo')
+% save('Tpo','Tpo')
+
+%% QSOs Patch points
+
+% % Portforlio di orbite
+% load("Xpo.mat")
+% load("Tpo.mat")
+% n_rev = 5;
+%
+% % QSO-H
+% % idx = 1;
+%
+% % QSO-M
+% % idx = 53;
+%
+% % QSO-La
+% % idx = 79;
+%
+% % % QSO-Lb
+% idx = 86;
+%
+% % QSO-Lc
+% % idx = 90;
+%
+% % Definition patch points
+% Phi0 = eye(6);
+% X0 = Xpo(:,idx);
+% theta0 = zeros(6,1);
+% X0 = reshape([X0 Phi0 theta0],[48,1]);
+% tspan = [0 n_rev*Tpo(idx)];
+% ODE.T = 1;
+% ODE.lambda = 0;
+% RelTol = 1e-13;
+% AbsTol = 1e-16;
+% opt = odeset('RelTol',RelTol,'AbsTol',AbsTol,'event', @(t,X) Patch_points(t,X));
+% [t,X,te,Xe] = ode113(@(t,X) HillProb_Sphere_Nd(t,X,Phobos,ODE),tspan,X0,opt);
+% X = X(:,1:6).*[1,1,1,T,T,T]*L;
+%
+% % Xe = (Xe(:,1:6).*[1,1,1,1/T,1/T,1/T]*L)';
+% % te = te(:,1)*T;
+%
+% Xe = (Xe(1:2:end,1:6).*[1,1,1,1/T,1/T,1/T]*L)';
+% te = te(1:2:end,1)*T;
+%
+% % Patch points
+% X_patch = Xe(1:3,:);
+% V_patch = Xe(4:6,:);
+% t_patch = te;
+%
+% % Plot
+% figure()
+% hold on
+% grid on
+% axis equal
+% planetPlot('asteroid',[0,0,0],Geom,1);
+% xlabel('X')
+% ylabel('Y')
+% plot3(X(:,1),X(:,2),X(:,3));
+% plot3(Xe(1,:),Xe(2,:),Xe(3,:),'r*')
+
+
+
+
+%% 3D Patch points
+
+% % Qualche costante utile
+% alpha = 13.03; %km
+% beta = 11.4; %km
+% gamma = 9.14; %km
+% Geom.alpha = alpha;
+% Geom.beta = beta;
+% Geom.gamma = gamma;
+%
+% % Qualche parametro e units utile
+% [pars, units] = MMX_DefineModelParametersAndUnits([2,2]);
+%
+% inp.Model.dynamics = @CRTBPsh;
+% inp.Model.pars = pars;
+% inp.Model.units = units;
+% inp.Model.epoch = '2026-03-16, 00:00:00 (UTC)';
+% inp.OrbitType = '3D-QSO';
+% inp.Event = @Null_Radial_velocity;
+%
+% % QSO-H
+% % inp.Amplitudes = [198, 0];
+% % QSO-M
+% % inp.Amplitudes = [94, 0];
+% % QSO-La
+% % inp.Amplitudes = [49, 0];
+% % QSO-Lb
+% % inp.Amplitudes = [31, 0];
+% % QSO-Lb
+% inp.Amplitudes = [27, 0];
+%
+%
+% % Propagazione e caricamento dei dati della orbita di interesse
+% filename = sprintf('3DQSOs_0%d_%dx%d.mat',round(inp.Amplitudes(1)),...
+% size(inp.Model.pars.SH.Clm,1)-1, size(inp.Model.pars.SH.Clm,2)-1);
+% load(filename, 'Xqp', 'Wqp', 'GMOS', 'Model');
+%
+%
+% % Integrazione e identificazione dei patch points
+% idx = 10;
+% Ic = Xqp(1:6,idx);
+% Model.pars.T = Wqp(1,idx)*units.tsf;
+% N_rev = 3;
+% Time = [0, N_rev*Model.pars.T/units.tsf];
+% AbsTol = 1e-16;
+% RelTol = 1e-14;
+% opt = odeset('RelTol', RelTol, 'AbsTol', AbsTol, 'event',...
+% @Null_Radial_velocity);
+% [t,x,te,xe] = ode113(@eom, Time, Ic, opt, Model);
+% t_patch = te(2:2:end)*units.tsf;
+% X_patch = xe(2:2:end,1:3)'.*units.lsf;
+% V_patch = xe(2:2:end,4:6)'.*units.vsf;
+
+
+%% SWING Patch points
+
+% Qualche costante utile
+ alpha = 13.03; %km
+ beta = 11.4; %km
+ gamma = 9.14; %km
+ Geom.alpha = alpha;
+ Geom.beta = beta;
+ Geom.gamma = gamma;
+
+% Qualche parametro e units utile
+ [pars, units] = MMX_DefineModelParametersAndUnits([2,2]);
+
+ inp.Model.dynamics = @CRTBPsh;
+ inp.Model.pars = pars;
+ inp.Model.units = units;
+ inp.Model.epoch = '2026-03-16, 00:00:00 (UTC)';
+ inp.OrbitType = 'Swing-QSO';
+ inp.Event = @Null_Radial_velocity;
+
+% QSO-H
+% inp.Amplitudes = [198, 0];
+% QSO-M
+% inp.Amplitudes = [94, 0];
+% QSO-La
+% inp.Amplitudes = [49, 0];
+% QSO-Lb
+ % inp.Amplitudes = [31, 0];
+% QSO-Lc
+ inp.Amplitudes = [27, 0];
+
+
+% Propagazione e caricamento dei dati della orbita di interesse
+ filename = sprintf('SwingQSOs_0%d_%dx%d.mat',round(inp.Amplitudes(1)),...
+ size(inp.Model.pars.SH.Clm,1)-1, size(inp.Model.pars.SH.Clm,2)-1);
+ load(filename, 'Xqp', 'Wqp', 'Model');
+
+% Integrazione e identificazione dei patch points
+ idx = 12;
+ Ic = Xqp(1:6,idx);
+ Model.pars.T = Wqp(1,idx)*units.tsf;
+ N_rev = 2.5;
+ Time = [0, N_rev*Model.pars.T/units.tsf];
+ AbsTol = 1e-16;
+ RelTol = 1e-14;
+ opt = odeset('RelTol', RelTol, 'AbsTol', AbsTol, 'event',...
+ @Null_Radial_velocity);
+ [t,x,te,xe] = ode113(@eom, Time, Ic, opt, Model);
+
+ % t_patch = te(2:2:end)*units.tsf;
+ % X_patch = xe(2:2:end,1:3)'.*units.lsf;
+ % V_patch = xe(2:2:end,4:6)'.*units.vsf;
+
+ t_patch = te*units.tsf;
+ X_patch = xe(:,1:3)'.*units.lsf;
+ V_patch = xe(:,4:6)'.*units.vsf;
+
+% Plot orbita e patch points
+ h = figure();
+ PlotTrajectories(h, x, inp.Model);
+ hold on;
+ plot3(X_patch(1,:),X_patch(2,:),X_patch(3,:),'ro')
+ set(0,'DefaultFigureVisible', 'on')
+
+
+
+%% Phobos states
+
+% Now, this points must be defined with respect to the new model. So we
+% integrate Phobos with the new model, and add the MMX postition with
+% respect to this new Phobos
+
+% kernel with the real position of Phobos
+ cspice_furnsh(which('mar097.bsp'));
+
+% Model parameters
+ [pars, units] = MMX_DefineNewModelParametersAndUnits_Right;
+
+% Initial state
+ data = '2026-03-16 00:00:00 (UTC)';
+ data = cspice_str2et(data);
+ pars.et0 = data;
+ [Ph,pars] = Phobos_States_NewModel(data,pars);
+
+% Initial Phobos's state vector
+ Ph0 = Ph./units.sfVec2;
+ Iz = pars.IPhz_bar + pars.ni*Ph0(1)^2;
+ K = Iz*Ph0(4) + pars.IPhz_bar*Ph0(8);
+ pars.K = K;
+
+
+% Initial Phobos's state vector
+ St0 = [Ph(1:2); Ph(7:8)];
+
+% Integration
+ tspanp = t_patch*units.tsf;
+ opt = odeset('RelTol',RelTol,'AbsTol',AbsTol);
+ [tp,XPhp] = ode113(@(t,X) MPH_Relative(t,X,pars),tspanp,St0,opt);
+ tp = tp/units.tsf;
+
+
+
+%% Setup della optimization per la definizione dell'orbita periodica nel nuovo modello
+
+% Definition of the initial vector
+ X0 = [X_patch/units.lsf; V_patch/units.vsf;...
+ (t_patch(2)-t_patch(1))/2*ones(2,size(t_patch,1))*units.tsf;
+ t_patch'*units.tsf];
+ X_Ph0 = XPhp(:,1:end)';
+ pars.Geom = Geom;
+ pars.X0 = X_patch/units.lsf;
+
+%% Propagazione della vecchia initial condition
+
+ [~] = Plot_InitialCondition(X0,[],[],pars,units);
+
+
+%% Ottimizzazione
+
+% Definition of the constraints
+% There are no linear inequality constraints nor linear equality
+% constraints, but only non-linear constraints that are specified in the
+% function @Constraints
+
+ options = optimoptions('fmincon','Display','iter','Algorithm','interior-point',...
+ 'ConstraintTolerance',1e-9,...
+ 'MaxFunctionEvaluations',1e4,'EnableFeasibilityMode',true,...
+ 'OutputFcn',@(x,optimValues,state)Plot_InitialCondition(x,optimValues,state,pars,units));
+ X_good = fmincon(@(x) CostFunction(x,pars,units),X0,[],[],[],[],[],[], ...
+ @(x) Constraints(x,pars,units,X_patch),options);
+
+ save('SwingQSO-Lc_5rev_right.mat','X_good')
+
+
+%% Test della nuova initial condition found
+
+ [~] = Plot_InitialCondition(X_good,[],[],pars,units);
+
diff --git a/New Model/SwingQSO-Lb_5rev_right.mat b/New Model/SwingQSO-Lb_5rev_right.mat
new file mode 100644
index 0000000000000000000000000000000000000000..b86fcd5c3e3ca0191b20449bfc02e38cb79e616c
Binary files /dev/null and b/New Model/SwingQSO-Lb_5rev_right.mat differ
diff --git a/New Model/SwingQSO-Lc_5rev_right.mat b/New Model/SwingQSO-Lc_5rev_right.mat
new file mode 100644
index 0000000000000000000000000000000000000000..77a0a68a0a9f15fece63f7234c7bec228951b1ff
Binary files /dev/null and b/New Model/SwingQSO-Lc_5rev_right.mat differ
diff --git a/New Model/Test_relativeDyn.m b/New Model/Test_relativeDyn.m
index f8633ea5238bed8918426163f9f8a61a989efbeb..47008ecbad8975580682915d143332acc9a30c49 100644
--- a/New Model/Test_relativeDyn.m
+++ b/New Model/Test_relativeDyn.m
@@ -10,6 +10,7 @@ clc
% Mac
restoredefaultpath
+ addpath('./MMX_BSP_Files_Relative/');
addpath('../../useful_functions/');
addpath('../../useful_functions/planet_textures/');
addpath('../../dynamical_model/');
@@ -19,17 +20,20 @@ clc
addpath(genpath('../../generic_kernels/'));
addpath(genpath('../../paper/MMX_Fcn_CovarianceAnalyses/'));
addpath(genpath('../../paper/MMX_Product/MMX_BSP_Files_GravLib/'));
+ cspice_kclear
MMX_InitializeSPICE
cspice_furnsh('MARPHOFRZ.txt');
- cspice_furnsh('PHMARFRZ.txt');
cspice_furnsh(which('mar097.bsp'));
- cspice_furnsh(which('MMX_QSO_198_2x2_826891269_828619269.bsp'));
+ % cspice_furnsh(which('MMX_QSO_031_2x2_826891269_828619269.bsp'));
% cspice_furnsh(which('MMX_3DQSO_031_009_2x2_826891269_828619269.bsp'));
% cspice_furnsh(which('MMX_SwingQSO_094_006_2x2_826891269_828619269.bsp'));
cspice_furnsh(which('Phobos_826891269_828619269.bsp'));
-
+ % cspice_furnsh(which('MMX_QSO_027_2x2_826891269_829483269_relative_right.bsp'));
+ cspice_furnsh(which('MMX_3DQSO_027_014_2x2_826891269_829483269_relative_right.bsp'));
+ % cspice_furnsh(which('MMX_SwingQSO_027_011_2x2_826891269_829483269_relative_right.bsp'));
+
% Model parsameters
- [pars, units] = MMX_DefineNewModelParametersAndUnits;
+ [pars, units] = MMX_DefineNewModelParametersAndUnits_Right;
alpha = 13.03; %km
beta = 11.4; %km
gamma = 9.14; %km
@@ -55,15 +59,41 @@ clc
pars.K = K;
% Initial MMX's State vector
- MMX0 = cspice_spkezr('-34', data, 'IAU_PHOBOS', 'none', '-401');
- MMX0 = MMX0./units.sfVec;
+ % MMX0 = cspice_spkezr('-34', data, 'MARSIAU', 'none', '499')./units.sfVec;
+ % MMX0(1:3) = pars.MARSIAU2perifocal*MMX0(1:3);
+ % MMX0(4:6) = pars.MARSIAU2perifocal*MMX0(4:6);
+ % theta = Ph0(3);
+ % Xsi = Ph0(7);
+ %
+ % % From perifocal 2 PA
+ % Rtheta = [cos(theta), sin(theta), 0; -sin(theta), cos(theta), 0; 0, 0, 1];
+ % Rxsi = [cos(Xsi), sin(Xsi), 0; -sin(Xsi), cos(Xsi), 0; 0, 0, 1];
+ % Rtot = Rxsi*Rtheta;
+ % Phobos0 = cspice_spkezr('-401', data, 'MARSIAU', 'none', '499')./units.sfVec;
+ % r_Phobos = pars.MARSIAU2perifocal*Phobos0(1:3);
+ % v_Phobos = pars.MARSIAU2perifocal*Phobos0(4:6);
+ %
+ % mirror = [-1, 0, 0; 0, -1, 0; 0, 0, 1];
+ % MMX0(4:6) = mirror*(Rtot*(MMX0(4:6) - v_Phobos - ...
+ % cross([0; 0; Ph0(4) + Ph0(8)], MMX0(1:3) - r_Phobos)));
+ % MMX0(1:3) = mirror*Rtot*(MMX0(1:3) - r_Phobos);
+
+ MMX0 = cspice_spkezr('-35', data, 'IAU_PHOBOS', 'none', '-401')./units.sfVec;
+ MMX0(1:3) = pars.PB*MMX0(1:3);
+ MMX0(4:6) = pars.PB*MMX0(4:6);
+
+ MMX1 = cspice_spkezr('-35', data, 'IAU_PHOBOS', 'none', '-401')./units.sfVec;
+ MMX1(1:3) = pars.PB*MMX1(1:3);
+ MMX1(4:6) = pars.PB*MMX1(4:6);
+
% Analysis initial state vector
St0 = [MMX0; Ph0; pars.IPhx_bar; pars.IPhy_bar; pars.IPhz_bar];
- St02 = [MMX0; Ph0(1:2); Ph0(7:8); pars.IPhx_bar; pars.IPhy_bar; pars.IPhz_bar];
+ % St01 = [MMX1; Ph0; pars.IPhx_bar; pars.IPhy_bar; pars.IPhz_bar];
+ St02 = [MMX1; Ph0(1:2); Ph0(7:8); pars.IPhx_bar; pars.IPhy_bar; pars.IPhz_bar];
% Integration
- tspan = (0:freq:2*day)*units.tsf;
+ tspan = (0:freq:1*day)*units.tsf;
RelTol = 1e-13;
AbsTol = 1e-16;
opt = odeset('RelTol',RelTol,'AbsTol',AbsTol,'event', @(t,X) landing_Phobos_relative(t,X,pars,units));
@@ -88,107 +118,57 @@ clc
ylabel('Along Track [km]','FontSize',26,'Interpreter','latex')
xlabel('Radial [km]','FontSize',26,'Interpreter','latex')
zlabel('Z [km]','FontSize',26,'Interpreter','latex')
-
-
-% figure(2)
-% subplot(2,3,1)
-% plot(t/3600, X(:,1)-X1(:,1))
-% grid on
-% subplot(2,3,2)
-% plot(t/3600, X(:,2)-X1(:,2))
-% grid on
-% subplot(2,3,3)
-% plot(t/3600, X(:,3)-X1(:,3))
-% grid on
-% subplot(2,3,4)
-% plot(t/3600, X(:,4)-X1(:,4))
-% grid on
-% subplot(2,3,5)
-% plot(t/3600, X(:,5)-X1(:,5))
-% grid on
-% subplot(2,3,6)
-% plot(t/3600, X(:,6)-X1(:,6))
-% grid on
-%
-% figure(3)
-% subplot(2,3,1)
-% plot(t/3600, X(:,1))
-% hold on
-% grid on
-% plot(t1/3600, X1(:,1))
-% subplot(2,3,2)
-% plot(t/3600, X(:,2))
-% hold on
-% grid on
-% plot(t1/3600, X1(:,2))
-% subplot(2,3,3)
-% plot(t/3600, X(:,3))
-% hold on
-% grid on
-% plot(t1/3600, X1(:,3))
-% subplot(2,3,4)
-% plot(t/3600, X(:,4))
-% hold on
-% grid on
-% plot(t1/3600, X1(:,4))
-% subplot(2,3,5)
-% plot(t/3600, X(:,5))
-% hold on
-% grid on
-% plot(t1/3600, X1(:,5))
-% subplot(2,3,6)
-% plot(t/3600, X(:,6))
-% hold on
-% grid on
-% plot(t1/3600, X1(:,6))
-%
-%
-% figure(4)
-% subplot(2,4,1)
-% plot(t/3600, X(:,7)-X1(:,7))
-% grid on
-% subplot(2,4,2)
-% plot(t/3600, X(:,8)-X1(:,8))
-% grid on
-% subplot(2,4,3)
-% plot(t/3600, X(:,9)-X1(:,9))
-% grid on
-% subplot(2,4,4)
-% plot(t/3600, X(:,10)-X1(:,10))
-% grid on
-% subplot(2,4,5)
-% plot(t/3600, X(:,11)-X1(:,11))
-% grid on
-% subplot(2,4,6)
-% plot(t/3600, X(:,12)-X1(:,12))
-% grid on
-% subplot(2,4,7)
-% plot(t/3600, X(:,13)-X1(:,13))
-% grid on
-% subplot(2,4,8)
-% plot(t/3600, X(:,14)-X1(:,14))
-% grid on
-
+ legend('Dynamics 1','Dynamics 2')
+ X_real = cspice_spkezr('-35', data+t', 'IAU_PHOBOS', 'none', '-401');
+ X_real = (pars.PB*X_real(1:3,:))';
+ figure(2)
+ plot3(X(:,1)*units.lsf,X(:,2)*units.lsf,X(:,3)*units.lsf,'LineWidth',1.3)
+ hold on;
+ grid on;
+ plot3(X_real(:,1),X_real(:,2),X_real(:,3),'LineWidth',1.3)
+ planetPlot('asteroid',[0;0;0],pars.Phobos,1);
+ view(2)
+ axis equal
+ ylabel('Along Track [km]','FontSize',26,'Interpreter','latex')
+ xlabel('Radial [km]','FontSize',26,'Interpreter','latex')
+ zlabel('Z [km]','FontSize',26,'Interpreter','latex')
+ legend('Propagated','BSP')
+
-%% Difference wrt SPICE
+%% Phobos
-
- % MMX_t = cspice_spkezr('-34', data+t', 'MARSIAU', 'none', '499')./units.sfVec;
- %
- %
- %
% figure()
- % plot(t/3600, X(:,1) - MMX_t(1,:)')
+ % subplot(2,2,1)
+ % plot(t, X(:,7)*units.lsf,"LineWidth",1.3)
+ % hold on
% grid on
+ % plot(t, X1(:,7)*units.lsf,"LineWidth",1.3)
+ % xlabel('Epoch [s]','FontSize',26,'Interpreter','latex')
+ % ylabel('$R_{Ph}$ [km]','FontSize',26,'Interpreter','latex')
+ % legend('Mars centered', 'Phobos centered','FontSize',26,'Interpreter','latex')
+ % subplot(2,2,2)
+ % plot(t, X(:,13),"LineWidth",1.3)
% hold on
- % plot(t/3600, X(:,2) - MMX_t(2,:)')
- % plot(t/3600, X(:,3) - MMX_t(3,:)')
- %
- % figure()
- % plot(t/3600, X(:,4) - MMX_t(4,:)')
% grid on
+ % plot(t, X1(:,9),"LineWidth",1.3)
+ % xlabel('Epoch [s]','FontSize',26,'Interpreter','latex')
+ % ylabel('$\Phi_{Ph}$ [rad]','FontSize',26,'Interpreter','latex')
+ % legend('Mars centered', 'Phobos centered','FontSize',26,'Interpreter','latex')
+ % subplot(2,2,3)
+ % plot(t, X(:,8)*units.lsf,"LineWidth",1.3)
% hold on
- % plot(t/3600, X(:,5) - MMX_t(5,:)')
- % plot(t/3600, X(:,6) - MMX_t(6,:)')
- %
\ No newline at end of file
+ % grid on
+ % plot(t, X1(:,8)*units.lsf,"LineWidth",1.3)
+ % xlabel('Epoch [s]','FontSize',26,'Interpreter','latex')
+ % ylabel('$\dot{R}_{Ph}$ [km/s]','FontSize',26,'Interpreter','latex')
+ % legend('Mars centered', 'Phobos centered','FontSize',26,'Interpreter','latex')
+ % subplot(2,2,4)
+ % plot(t, X(:,14),"LineWidth",1.3)
+ % hold on
+ % grid on
+ % plot(t, X1(:,10),"LineWidth",1.3)
+ % xlabel('Epoch [s]','FontSize',26,'Interpreter','latex')
+ % ylabel('$\dot{\Phi}_{Ph}$ [rad/s]','FontSize',26,'Interpreter','latex')
+ % legend('Mars centered', 'Phobos centered','FontSize',26,'Interpreter','latex')
+
\ No newline at end of file
diff --git a/New Model/landing_Phobos_relative.m b/New Model/landing_Phobos_relative.m
index 8d786d1b5db82ffb677b39924cd39bf1269ac33c..67dd654db9f6c62ec8000f37a009ecf44fdca2cc 100644
--- a/New Model/landing_Phobos_relative.m
+++ b/New Model/landing_Phobos_relative.m
@@ -15,7 +15,7 @@ function [value,isterminal,direction] = landing_Phobos_relative(~,X,par,units)
beta = par.Mars.beta/units.lsf; % Mars' intermediate semi-axis (km)
gamma = par.Mars.gamma/units.lsf; % Mars' smallest semi-axis (km)
- rPh = X(7)*[cos(X(13)); sin(X(13)); 0];
+ rPh = X(7)*[cos(-X(13)); sin(-X(13)); 0];
value2 = (x(1) - rPh(1))^2/alpha^2 + (x(2) - rPh(2))^2/beta^2 +...
(x(3) - rPh(3))^2/gamma^2 - 1;
diff --git a/Observables_model_with_Features.m b/Observables_model_with_Features.m
index ccd87e68521f5b1d72bf2b3aff6adef3e88c0787..0a74e4ab7ac4e91e633229e9e867f382ffa68d72 100644
--- a/Observables_model_with_Features.m
+++ b/Observables_model_with_Features.m
@@ -189,7 +189,7 @@ function [G, Rm] = Observables_model_with_Features(et, X, St_ID_Range, St_ID_Rat
R2Rot = [cos(X(9)), sin(X(9)), 0; -sin(X(9)), cos(X(9)), 0; 0, 0, 1];
v_Phobos = pars.perifocal2MARSIAU*(R2Rot'*([X(8); 0; 0]) + ...
- cross([0; 0; X(10)], [X(7)*cos(X(9)); X(7)*sin(X(9)); 0]));
+ cross([0; 0; X(10)], X(7)*[cos(X(9)); sin(X(9)); 0]));
Phobos = [r_Phobos.*units.lsf; v_Phobos.*units.vsf];
[~, Y_pix] = visible_features_model(X_MMX.*units.sfVec, Phobos,...
@@ -220,7 +220,7 @@ function [G, Rm] = Observables_model_with_Features(et, X, St_ID_Range, St_ID_Rat
R2Rot = [cos(X(9)), sin(X(9)), 0; -sin(X(9)), cos(X(9)), 0; 0, 0, 1];
v_Phobos = pars.perifocal2MARSIAU*(R2Rot'*([X(8); 0; 0]) + ...
- cross([0; 0; X(10)], [X(7)*cos(X(9)); X(7)*sin(X(9)); 0]));
+ cross([0; 0; X(10)], X(7)*[cos(X(9)); sin(X(9)); 0]));
Phobos = [r_Phobos.*units.lsf; v_Phobos.*units.vsf];
Limb = Mars_LimbRange_model(X_MMX.*units.sfVec, Phobos, Sun, Xsi, pars, units);
@@ -243,7 +243,7 @@ function [G, Rm] = Observables_model_with_Features(et, X, St_ID_Range, St_ID_Rat
R2Rot = [cos(X(9)), sin(X(9)), 0; -sin(X(9)), cos(X(9)), 0; 0, 0, 1];
v_Phobos = pars.perifocal2MARSIAU*(R2Rot'*([X(8); 0; 0]) + ...
- cross([0; 0; X(10)], [X(7)*cos(X(9)); X(7)*sin(X(9)); 0]));
+ cross([0; 0; X(10)], X(7)*[cos(X(9)); sin(X(9)); 0]));
Phobos = [r_Phobos.*units.lsf; v_Phobos.*units.vsf];