From e55576cce7588ca74ed22fbb93a54c749562ce29 Mon Sep 17 00:00:00 2001
From: Ciccarelli <ec01259@fvffw3gnq05p.home>
Date: Wed, 4 Oct 2023 10:07:11 +0100
Subject: [PATCH] Orbite rigenerate
---
New Model/3DQSO-Lb_5rev_right.mat | Bin 0 -> 912 bytes
New Model/3DQSO-Lc_5rev_right.mat | Bin 0 -> 616 bytes
New Model/BSP_generator.m | 144 ++++++++
New Model/Constraints.m | 74 ++++
New Model/CostFunction.m | 35 ++
New Model/Dynamics_MPHandMMX_Relative.m | 41 +--
New Model/Dynamics_MPHandMMX_Relative2.m | 30 +-
...X_DefineNewModelParametersAndUnits_Right.m | 108 ++++++
New Model/MMX_GenerateBSPFiles_Relative.m | 31 ++
New Model/MPH_Relative.m | 64 ++++
New Model/Plot_InitialCondition.m | 87 +++++
New Model/QSO-Lb_5rev.mat | Bin 0 -> 742 bytes
New Model/QSO-Lb_5rev_right.mat | Bin 0 -> 886 bytes
New Model/QSO-Lc_5rev_right.mat | Bin 0 -> 888 bytes
New Model/QSOs_Definition.m | 322 ++++++++++++++++++
New Model/SwingQSO-Lb_5rev_right.mat | Bin 0 -> 837 bytes
New Model/SwingQSO-Lc_5rev_right.mat | Bin 0 -> 834 bytes
New Model/Test_relativeDyn.m | 188 +++++-----
New Model/landing_Phobos_relative.m | 2 +-
Observables_model_with_Features.m | 6 +-
20 files changed, 992 insertions(+), 140 deletions(-)
create mode 100644 New Model/3DQSO-Lb_5rev_right.mat
create mode 100644 New Model/3DQSO-Lc_5rev_right.mat
create mode 100644 New Model/BSP_generator.m
create mode 100644 New Model/Constraints.m
create mode 100644 New Model/CostFunction.m
create mode 100644 New Model/MMX_DefineNewModelParametersAndUnits_Right.m
create mode 100644 New Model/MMX_GenerateBSPFiles_Relative.m
create mode 100644 New Model/MPH_Relative.m
create mode 100644 New Model/Plot_InitialCondition.m
create mode 100644 New Model/QSO-Lb_5rev.mat
create mode 100644 New Model/QSO-Lb_5rev_right.mat
create mode 100644 New Model/QSO-Lc_5rev_right.mat
create mode 100644 New Model/QSOs_Definition.m
create mode 100644 New Model/SwingQSO-Lb_5rev_right.mat
create mode 100644 New Model/SwingQSO-Lc_5rev_right.mat
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
GIT binary patch
literal 912
zcmeZu4DoSvQZUssQ1EpO(M`+DN!3vZ$Vn_o%P-2c0*X0%nwjV*I2WZRmZYXA<mXu_
zgp{T#_$QYrC>Sdk8d(_{TNxWF7#SEDD-a79V1UunmmkRHU}j*bm}7gXR_98g1j7gU
z>pK@a^UT~nQKwvD^%artf9HB7#XRHCeP5&;duqb<PuoIFo~o)XnfcP-+0*_<`#zl$
zv)3znw)XS7LvyG9c0G}3dwX;G`Oj}YU;nJ3$l;`>&@pF=QhTV>k;QR!_M+c{pJhCM
zD17E2i`+TGNEX3d%UkJ{k00#WznwdH15b{y`tzQ%ZVs7;4!^w+#I@p5<Tb@T3l9mU
z*JeL<nUg;2`Yfl;g8}95f0yz;l&Es`FArE`(NU)P^v4ge`~z{jrR_VTv;*$@o4k?w
z#4PmR?ETtzs~keDSNGjWNWP#aa`o`dvbe?ZR~q?uWhO+vZd`b5&hpl+7rR#^oyvU6
zSk74?o4-VV;~eIFeA4xHt$j?-SmXBH;p8$9O8+kw5}2rSDf#=8!)CR)3opgIx_#=3
zulWAt7ninnJh;_bTC#P+Bdf?3=XLYr&Zj;Sx*~Qq>GRp80n?K%rWU$0`qfOl86v%u
z>3@=zqOF4%SI`Z%p!o;lwlB5dUl1u#zq)5@^DGYOzo!qYG%R<%aQ0I)+q{`(iF$_4
z1`Fm3%1&6#!}R*jvAOxdimycS&u86U#C&)~aghAGX=yLij$IHt-#z7_Y^0l4QGc>u
zBTL);x7uqBMGj5#i;KB%No2;2%fbHV)DMTe{kg=Dk5_@uNxSW*xtyl3SH$l+-f#NV
zt_dr4x^WtR5$Zgs{{7afu!jyd|7S>MDzX$G%6n?~$s(j@hVb-`6vunraXZWOc>P{I
z*tYnD+WT4E*|V=KUa-d5F!_JpMdKq&qK``!#Y!h0vXEJ~B>QSQ|JlQn8K<RBiu8Rj
zb=S(?LtLDGg^qjM#CSVity^^_?-axRHCx11eqU90?dcM$hPDoECF`UqlCzk~eB%ZA
zSa?^)bRFa|unars_xO8rdDFy~@vff|rv6LNU)I3pne}Sj`L*ShNBTwA2JPQ6>6S^R
zSwr2m!}6his|==iEVoeewOI0b|0mx`tsivEH|UzDa|vZ<8Q=TFIYBD2j{oJ7##i6(
e&Z?HW{N>!2Wz7Ho%fF3l+tiuWQ?K2W=K}z@T$&63
literal 0
HcmV?d00001
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
GIT binary patch
literal 616
zcmeZu4DoSvQZUssQ1EpO(M`+DN!3vZ$Vn_o%P-2c0*X0%nwjV*I2WZRmZYXA<mXu_
zgp{T#_$QYrC>Sdk8d;f`S(z9q7#SEDD-a79V1UunmmkP}z{tQ*G3Rk|LPEvC69*2k
zEwPkPFyYwX^m9&tp__)R;QV>>=Q9c3XAp_`{&e!wIp^MgyZxh#Wp93E+Wg-?J6<jF
zY=0{~OUwT5uG3lZl@8Hak5$k9z4cv9F7EpG>ieG)zQ5BId@i$Zo2T7e!HRQ<wW%*n
z-|v4eTcmAL7kRDb?){|;f4=d4dHL_r>tY4>KJVRVAM*IYqwmI*-%snTU$u79wJrRd
z_IBEp=aP3ny}jwcugYt;>r1v;J}<DmuIC*0c&GCP_lsHA&*XnwwYKN3n4bR|?s+9L
z|1K9Dj(c{^|Nj%8yx*(4S!dNKn|}X%YL?5_v-Qd830uD#d;GsTeKYqm8LQ_C6{+(Z
z7kx@+|MzP9uKl(9-uCbL^7V4y-R-~Cr&Me{r+sVt{PRiAY#Y;Seh25D$v$=U@`mqL
z+h^)Ge|&E@S@zLUg}<kJHZQj`vwJV_B=J@Cn@{@eQ{PYienhoqYrN{@ieF#OU;h@~
zSFJfU)^6MWKaORejx#w{^S66`-1J83ccp3W+IvU;Z<%eH{=H<vKj)?F70>QiK3||c
z!`j%c;bP#!Ppq8Vz2m>d*Dw58ezdpzcYgidU!QKzmi|`fE?XM+tNYFTU9WFwJzx3i
k+p)Lli%g61-+cbPk#CK>*u5GVolEn#M#?g7I{l~_0J8=y`~Uy|
literal 0
HcmV?d00001
diff --git a/New Model/BSP_generator.m b/New Model/BSP_generator.m
new file mode 100644
index 0000000..6a24347
--- /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 0000000..3fc9ded
--- /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 0000000..eb6d456
--- /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 6754cc7..8152044 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 6fa278a..da475f7 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 0000000..c9ba3bb
--- /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 0000000..1240fd7
--- /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 0000000..8d23469
--- /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 0000000..719d0be
--- /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
GIT binary patch
literal 742
zcmeZu4DoSvQZUssQ1EpO(M`+DN!3vZ$Vn_o%P-2c0*X0%nwjV*I2WZRmZYXA<mXu_
zxD{n81g91#7+ES9npl|{S{WNE7#SEDD-a79V1UunmmkQEV`5;anDaO}Az?*AQbGdr
zB^e$8DVB(ZHGUdrIi_}Z_@6)T-`HWz@b&H9oTYVd-{=1>7GJ-;HG9sc^Yd*}V{fjR
z@%QtV<;VZ8ug$!ke<;@EZTGu>$IjJHVmdGO=62noSC<aY>E?SUeS^R8#k8|`<{0Hg
ze)C)Uy(lf~&&~YBJ3rq!7G3?_I@Dt6mblY*>brg)IX1mJ?!EOwo9xS4oVLHEndF=G
zUYy-N+t8h9js3U8i<|<nFGb3pC!TX|FEx_hm#hCbd4KNrM(<7h-*1VV<Z3+^Uc0^i
z_>O}IV_(gZonC%xpLSK%@A*=jQs!NmvBlr?-+Rvq|DHJi*K0UyY{OS_-@W}%@9gE8
zQh#pN7k_>0`<Y+%=KFiaj;9|bU;bCyy{B)^y^`$tv#YAQFDJ)V-!#eIxBrjaZ!gpD
z4Sz$oZT}$C!B@Arw0y!RefEv)x-vX9^>5yvxuP+1*YgFn-<Nzo$6GA*MpANm_O(;>
z8>$Z{?0Wwq|L)(9nZGyZhN|6a{i)0RIy_5n7uO>3`t^H%7_?u!_-pl%jhUJC-||oX
zSQJp7dcQ7{xiR>Bo!Gl@<BFP7YL9R3i`~XOeO=+@&GFx#A9$Uf^ZDD$_v!n1F0DUU
zCHiLl_TSZamiTS3s`?qbF8<Ni<8jJ%<r}MS$^Bk`=IHBbGF*4+PhH-j=vBLTu2bUg
zWhSb7|MpDF@t0p8k#KxpJL8)8CwX^x1$1LJ)`U+mn!b2Z^uAZhAMJnt4Cl;>@4CIi
z{?^0)``4d(`t8$s@vb@N_HH}Bv$Xii9sOP3zWeXJU2;C@nw{+Hlnr$ef4tw{%fAqN
M{5c21ArAfh0ClKt?*IS*
literal 0
HcmV?d00001
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
GIT binary patch
literal 886
zcmeZu4DoSvQZUssQ1EpO(M`+DN!3vZ$Vn_o%P-2c0*X0%nwjV*I2WZRmZYXA<mXu_
zgp{T#_$QYrC>Sdk8d@2dTNxWE7#SEDD-a79V1UunmmkP}$Hc%;G3Rk|Lc)WDw1foa
zi#!5SED;N9{1l!sCAqDLajB`QtK)L{!z9ey{CBF$)a#o2-*l_R6z@Ip;Y3|g^rCNX
zj>eW9e^fSaRz>f@dz<5u8GqNm{X6lI+PCUgX1DU}6H8wvo^x$4JvU3+p`_~Vk&t^U
zb-xMmXU7{z{<*dOuS#)|`_H$36*({Ol6Cs|{b~A@zZRi~cieYhsd@g6&+&a7y3tF$
z@?XbI`@Mcz^OxI}qH9?%zQ6mZ&h|{@iw$oZKdx%~zyA&Iy!#)G&Gx0$XMSEg$2Dp1
z?|HB1M24jOxcjbqDNAzF>~sJ4Pae8h8GP-Qm9g<!=9lq@EM6OQ>`OdU{50{LJNvrj
zUGBTL{l8Ve?C|$f7kB;s9l3n{w~OoNXO;iES#^2;`TDok?=+U)JGD0XeAH4a<AaKc
z>)#(ep{(#XKUDX?7ZKlcxsIpy%iQ}`tRfr!CpG`4xzWdjXR^OzuQl4VO;>xjIU-Q#
z|J=rt{|c&I4*ONk{rgURum1UM5nHd{IQxXb<NsCbWe*BI1<xw^mB_dF@!j}F->B(z
zZ}$f)O+6-GTl?QWuPpu=^U^Zkt&&c+mHxl`z0Yd<%=%q@f9G##UGG2teo}Q*<tfn#
zb%x7B`5j7vg0?kZ;S6)%X)x_~fBe5q&#!;cdoq1T?)%MGKCe%WUw!x3e4hL1Q){;F
zHh%nTZ6>GH@6S(q|9}5?Z%6s^XV={miX~R%HN3F?_r07&v)FRq&f4Tjc_-KZb>n*N
zbKWiPetrHe(WUDr1%K1c{JwMd`OnYyy)ECqk@Mz@8_!)crtO(LGc{WE|8q5u|NDNZ
z>N|e1Otk&`&dToi%T4EQYxV!%B>rPt!P&`6vcFw5&a*GloR#!E{#fkxzX#&2mtQ%0
z`}FFPZ>y)|-8h!B>p+<2PSG|0r$+r`70})FLgm8zroQm-cbs;mPp{YMyZL>d{cZdA
z{r4x<>#DD+ysss>tT|N4{=4rvn_VYV{f^(u>#IBR?!nn=*XiO7Uk_h;9KKV!rgU?Q
KEMw=fB1Zr&@5877
literal 0
HcmV?d00001
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
GIT binary patch
literal 888
zcmeZu4DoSvQZUssQ1EpO(M`+DN!3vZ$Vn_o%P-2c0*X0%nwjV*I2WZRmZYXA<mUm!
z^79n@lS>p7j1&wltqd%!jEohG41ff&fB^;=J$?Cs><>%~3>9-8CnqF4NJvXaV7|yB
zAjJ}~u*Of}8B>zmiWrxgs=7KZmp@Fx%!j{k6<^<G@n^&TW4F>OcE4P(_oCkQ%BG^U
z-)i-9>Ymo!Q?`D8ZpM-Ofwz9^zGXeU%6#|xy|#a6{!+RBneVdtUeoy;OMW{Yit{&}
zn^Es>cGCLxv*Z81MBOi{Tg?1>`kLpN>#qfF(w<%1>iF#M>)QN(wXv0IuD9#5R`gBc
zjbD<>cV2$OKlA<5UNxUFzW(|1yeHN77Y07Pv^R34{hQr~XYKiM=i4pmv(ql_p0Y6M
zxT&jp?ecrf)_s5H*Pjl_de6_jHuU~dxsB!1o|)NB_}3!o#HMg7?BKbwN$;!K#nyf}
zy62z$x9OjvR@8Ug+im}K>;BLB$)9Uy_e`7Kb$0)b_XVfciT=O4^=<jSGij&#O16mp
zoVAem(*BhBd&E2JZlBohwf)s@XQ!S`H<MP&e>+ypvsCSy+?M;B?7rrQw$F@x6MF8@
zgjL~X`hEqGex+CLy)SrtVNU(Iylv%j6|9Y4>}N$BVp;S(ZZ2=h`TBO}e3q+zyZ_&i
z5Bd3_>h0y=->RzhITPot+p>OX!Dp-A`tjd7MML(~h%L|EZ+>FUyK|R2_O8_7lu`IQ
zbN+UYri}E-uYSDZJk$4gdGHVIGxay>H=Ou6Z&}%yeXoj+orw+7`YlrJ=hwBpGB5m}
z-~ZEZjNblT{<-J#%u=PCda*m{>@WXsbeb#bQXbZ+b7gH(xk&Aq@9_s;AO9CQ|5NL<
z?|C=vpB~qG>G}K2n_GXRD{4jk-&RNK2Tzq>oLj#u>zL{9nSRx$Ev6j!SDAa|kFZES
zw~Mv>J^yEAArGI-+iv%+y|&jU_pS5Z$KSRt&1Z_!F}PjsnYh*T@%QqwT`_BSZOXg#
z+x63O-yd(@@A{tg`_6lXzk>6_ZZ)0x95`pK!O7|OH$IDZN#Ay(|1Zz_+Gl@u*Cqa6
zD5`U_Hh#*Ro^@V-A7}mC_I`f(nUhM!xz*n^nuX4r3U1R6d(Zm1ed%$@bW<O(vtlJI
J4D%hMx&S{Dz@PvC
literal 0
HcmV?d00001
diff --git a/New Model/QSOs_Definition.m b/New Model/QSOs_Definition.m
new file mode 100644
index 0000000..11f0c24
--- /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
GIT binary patch
literal 837
zcmeZu4DoSvQZUssQ1EpO(M`+DN!3vZ$Vn_o%P-2c0*X0%nwjV*I2WZRmZYXA<mXu_
zgp{T#_$QYrC>Sdknpv3|SQ#297#SEDD-a79V1UunmmkR9%f!G?G3Rk|LW0GCBL@z!
zEwPkPFyYwX^m9(Y8IGyl9scLf`!{x!GJLgtelKd}^Iy03?dkk`<L2(?Z!TL@J=MRr
zbHn?eg^K%AKlyz7x?Cx4ug~S1@Bi9W`MXX#y>HK|w<|SQ_gPB#{z*FUtmMqY!X2mX
zmG3+*Ym|OEJ$&=`_w}DXAOF9#_YM20$oLm2%QW8Ke>HWp+l7$A|BKTn`?W8z-*(yb
zx53ga9!2pp&S!7<*V5drw)dIED}l4sA1dz7XXmSy_#&Zy``^9!C+oHI!rof`j%AmK
z-QB|S?ftCm=ru-DFUD249_(b@Y4uIKZ^@l~rWbF%>-=B8p)UO6Rh^@L)w0)ouFUzd
z^x@XS|B7!n@6s!IBYtsz?YjR)&CPL%w|l>CcYSLgzUR)TXL~p8eLv@L?(Ey^-hAHk
z-)fEf<&%0*`Buyx|Nos%ynbEo-s^-9-;Q<(OUlcB(dSNkQ}{E{zVhcPVb4{!-#*Wl
zyLbA;f*kv$Ps_d+HFlTT-QT13Kk?Gd|DF%skN1nT*3OQ*znlBg{_n52p59sfTjq*w
zZqOT%{Z>2X_|N*7Uj6XjrO@ts%jPFrEI9kgH}US5{bfQnmxK0g*La^_8!USu_D|CN
zonN?L=IP)0wsDTl8n#RGm*{PNS$o2cJ@Mgjnf$tc1so@0*Gi@RsNVl#y2#>XbBupS
z{?C8;uRNe`|6{Ai%yGRuzxV(7zbSjw=Dxp|`Tc|5<!`-bZk~53YeL;mm31rmPxo7i
z{J-~m{pS<%?EE|LUf-v8Z$tfk_1fy2k#G1SgSq_Ef1Y}i9#elTu&B3b?}t@iG<V<8
zzgb#8FYcXYdTzOOfl<VNp;Gtwnd^_s)ZKevP$ZI4GUuUY!P7Ti<G=kemkX=4DoU?^
z`DL#EbK{Pi`o^p8%Lp9YcHi)D-o@?8w|=W@n_qobk)3>8bc5&``DuM#GG%Y-KFzx+
K$S_;NuM_}a%c%1J
literal 0
HcmV?d00001
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
GIT binary patch
literal 834
zcmeZu4DoSvQZUssQ1EpO(M`+DN!3vZ$Vn_o%P-2c0*X0%nwjV*I2WZRmZYXA<mXu_
zgp{T#_$QYrC>SdkT3Q)fTA7$A7#SEDD-a79V1UunmmkR9#l*l+G3Rk|LW0GCBL@z!
zEwPkPFyYwX^m9(Y8IGyl9scLf`!{x!GJHK-bzJR3`XT<`d+pkUw`t7Vod4Foa7q5u
z#kJX?Ura31)P8?ix@Y6G|69J*f1Lj|;=|s9dPR1dzaBfG%%{KYOWG78xxyKW^`Cok
zJ_oXY%X$4Gm;Lv*^LO)a|9_+VnX&HV>@?QG!f&~ECO=;~L*w0jNs;~0?@!cPZxl8>
z9uq42P4x8lhOg(}O#fT);kaUZ{8anD5q$THXPD@!7JQo`|9j`E@2P8k)ZdDG<L5hL
z?WVQ1*-OrDYx>LctxxKfRU*@^Z^HU3ru?$cy`jGGJo_5^cl^O|*&p@V&oC6W9xso#
zm|__BQ%KCb+<ty<-`6jz_J8|x+y8&r9ML`Tzoxu96IwszZL<G|yXU>Fw!c5NzILkh
zwy<}*XCzJIz7(G;_bYY%=e+iRvv$2w4t*~hvtI7JkJyj2{cmortn`|F{l@+7$i2Q1
z3RU&o)4Oc1cNUlBdH;LTwK3!Ub?w+oNss0vZZB7S&OhPb?(5sc_odu;ToAqc>v`X+
z=lL!_c3oERcHW--*Uv_#Io-eicU~M%oByrU@BW(nG25TDe5zX)vGb~xz5JWo8|s?8
z?=Du`ef-W#SB1aQ_a)|Aq;I_482f6+weL&s%j`4LpDYmnCHa5V|63M{_k(|Lzgjd!
z_I1aW?N>~zjh>&{Upt#?`MsJy(Wm#VdT_}1<@%##Uw6&sZ~jtjd&1-2)1n)Rq2=0}
z72fN&&(xf4ZxvSZW`6kX(zoiu-G}eLlv_VbyzR#98Sj42FAF^Q?;Q93$azfSZ?mE&
z=UML6xSekQoF!ua`zt!-v)h#Uj_uz3{Kx&@&lRelO#d9WS^nFxU&bci>vnAW|0B5G
zjk&Yo@A-O}-AbwFcl|%d^7U|7Z`AK{jrW{}mErr=N0+Ka*z3sd^j_8a`hR)}3&Yc^
GPJ93s$e)=2
literal 0
HcmV?d00001
diff --git a/New Model/Test_relativeDyn.m b/New Model/Test_relativeDyn.m
index f8633ea..47008ec 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 8d786d1..67dd654 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 ccd87e6..0a74e4a 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];
--
GitLab