diff --git a/.DS_Store b/.DS_Store
index a085ba9f2eefe9ea7a47ac443a2bd50311b601f3..8c602a9effaeb0528f8725e0829ba377515ff5b1 100644
Binary files a/.DS_Store and b/.DS_Store differ
diff --git a/Analisys_With_Features.m b/Analisys_With_Features.m
index ca462a3660e0eadfeb97f9c3c1d0bdff2f45ce10..66cba209c9a7921685437a3e4b1894123fead91e 100644
--- a/Analisys_With_Features.m
+++ b/Analisys_With_Features.m
@@ -92,7 +92,6 @@ clc
% @Observables_model_with_Features,...
% par.R,YObs,par,units,file_features);
-
Results(Est, YObs, par, units)
% save('./Results/QSOLb_9days_AllObs_alpha8e1.mat', 'Est')
diff --git a/CovarianceAnalysisParameters.m b/CovarianceAnalysisParameters.m
index ba9bc4b3e2cb493684e77aec93622cd38faf8d50..a55c3880e0920a9b474557dee6095bb16e5234d3 100644
--- a/CovarianceAnalysisParameters.m
+++ b/CovarianceAnalysisParameters.m
@@ -75,8 +75,8 @@ pars.R(8,8) = pars.ObsNoise.pixel^2;
% Seconds between observation
pars.interval_Range = 600e10; % s
pars.interval_Range_rate = 600e10; % s
-pars.interval_lidar = 600; % s
-pars.interval_camera = 600; % s
+pars.interval_lidar = 300; % s
+pars.interval_camera = 300; % s
pars.flag_Limb = 1;
pars.flag_features = 1;
diff --git a/Mars_LimbRange.m b/Mars_LimbRange.m
index 5c20d67f8ebee058face44565ceb652de80aa9fb..ee039da050bf2dfb1016c923d7777a2f8950a94e 100644
--- a/Mars_LimbRange.m
+++ b/Mars_LimbRange.m
@@ -54,6 +54,7 @@ function Limb_Range = Mars_LimbRange(MMX, Phobos, Sun, pars, units)
Sun = R2rot*Sun(1:3);
I_Sun = Sun/norm(Sun);
+ angle_MMXSun = acos(dot(I_sM,I_Sun));
[~, ~, PM] = ProjMatrix(-r_sb_Ph, pars, units);
MP = PM*[-Mars_pole; 1];
@@ -63,11 +64,12 @@ function Limb_Range = Mars_LimbRange(MMX, Phobos, Sun, pars, units)
if ((any(MC<0))||(any(MP<0))||((MC(1)>pars.cam.nPixX))||...
(MP(1)>pars.cam.nPixX)||(MC(2)>pars.cam.nPixY)||...
- (MP(2)>pars.cam.nPixY))||(dot(r_sb_Ph, Mars_centre)>0)
+ (MP(2)>pars.cam.nPixY))||(dot(r_sb_Ph, Mars_centre)>0)||...
+ ((dot(I_sM,I_Sun)<0)&&((angle_MMXSun-pi/2)<pars.FOV))
Limb_Range = [];
else
Limb_radius = MP - MC;
- Limb_Range = (sqrt(Limb_radius(1)^2 + Limb_radius(2)^2) + random('Normal',0, pars.ObsNoise.pixel));
+ Limb_Range = round(sqrt(Limb_radius(1)^2 + Limb_radius(2)^2) + random('Normal',0, pars.ObsNoise.pixel));
% Plotfeatures_Pic(MC, MP, pars);
end
diff --git a/Mars_LimbRange_model.m b/Mars_LimbRange_model.m
index 67048418923f8deda884033bff80c859f2cc509b..5ce7d307fe1357684cbd387c7da692fb109315ac 100644
--- a/Mars_LimbRange_model.m
+++ b/Mars_LimbRange_model.m
@@ -61,16 +61,10 @@ function Limb_Range = Mars_LimbRange_model(MMX, Phobos, Sun, pars, units)
MC = PM*[-Mars_centre; 1];
MC = MC./repmat(MC(3),3,1);
- if ((any(MC<0))||(any(MP<0))||((MC(1)>pars.cam.nPixX))||...
- (MP(1)>pars.cam.nPixX)||(MC(2)>pars.cam.nPixY)||...
- (MP(2)>pars.cam.nPixY))||(dot(r_sb_Ph, Mars_centre)>0)
- Limb_Range = [];
- else
- Limb_radius = MP - MC;
- Limb_Range = sqrt(Limb_radius(1)^2 + Limb_radius(2)^2);
- % Plotfeatures_Pic(MC, MP, pars);
- end
+ Limb_radius = MP - MC;
+ Limb_Range = sqrt(Limb_radius(1)^2 + Limb_radius(2)^2);
+ % Plotfeatures_Pic(MC, MP, pars);
diff --git a/Observations_Generation.m b/Observations_Generation.m
index 20b7dbe7d02c3d2655d2e09b8e03040c5a104d5b..6311ab3c285172f1505a2fcf16815fe807850f6e 100644
--- a/Observations_Generation.m
+++ b/Observations_Generation.m
@@ -30,7 +30,7 @@ clc
% Time of the analysis
data = '2026-03-16 00:00:00 (UTC)';
data = cspice_str2et(data);
- Ndays = 3.5;
+ Ndays = 2;
t = Ndays*86400;
date_end = data+t;
diff --git a/Phobos_cartesian_covariance.m b/Phobos_cartesian_covariance.m
new file mode 100644
index 0000000000000000000000000000000000000000..828f72a1ea9f31488809286bc2ea82578dd0aaab
--- /dev/null
+++ b/Phobos_cartesian_covariance.m
@@ -0,0 +1,58 @@
+function [P_PhXYZ, P_PhXYZ_MMX, P_rel, P_PhVXYZ, P_Vrel] = Phobos_cartesian_covariance(X_t, P_t, pars, ~)
+
+
+% Partials delle cartesian coordinates wrt radial distance and true
+% anomaly
+
+% Position
+ dXdR = pars.perifocal2MARSIAU*[cos(X_t(9,:)); sin(X_t(9,:)); zeros(size(X_t(9,:)))];
+ dXdtheta = pars.perifocal2MARSIAU*[-sin(X_t(9,:)); cos(X_t(9,:)); zeros(size(X_t(9,:)))].*X_t(7,:);
+
+% Velocity
+ dVdR = zeros(3,size(X_t(9,:),2));
+ dVdR_dot = zeros(3,size(X_t(9,:),2));
+ dVdtheta = zeros(3,size(X_t(9,:),2));
+ dVdtheta_dot = zeros(3,size(X_t(9,:),2));
+
+ for i=1:size(X_t(9,:),2)
+
+ dVdR(:,i) = pars.perifocal2MARSIAU*cross([0;0;X_t(10,i)], [cos(X_t(9,i));sin(X_t(9,i));0]);
+ dVdR_dot(:,i) = pars.perifocal2MARSIAU*[cos(X_t(9,i)), sin(X_t(9,i)), 0; -sin(X_t(9,i)), cos(X_t(9,i)), 0; 0, 0, 1]'*([1; 0; 0]);
+ dVdtheta(:,i) = pars.perifocal2MARSIAU*([-sin(X_t(9,i)), cos(X_t(9,i)), 0; -cos(X_t(9,i)), -sin(X_t(9,i)), 0; 0, 0, 1]'*([X_t(8,i); 0; 0]) +...
+ cross([0; 0; X_t(10,i)], [cos(X_t(9,i)); sin(X_t(9,i)); 0])*X_t(7,i));
+ dVdtheta_dot(:,i) = pars.perifocal2MARSIAU*cross([0; 0; 1], [X_t(7,i)*cos(X_t(9,i)); X_t(7,i)*sin(X_t(9,i)); 0]);
+ end
+
+
+% Coomposizione della covariance P_PhXYZ
+
+ P_PhXYZ = zeros(3,3,size(X_t(9,:),2));
+ P_PhVXYZ = zeros(3,3,size(X_t(9,:),2));
+
+ for i=1:size(X_t(9,:),2)
+ P_PhXYZ(:,:,i) = [dXdR(:,i),dXdtheta(:,i)]*...
+ [P_t(7,7,i), P_t(7,9,i); P_t(9,7,i), P_t(9,9,i)]*...
+ [dXdR(:,i),dXdtheta(:,i)]';
+ P_PhVXYZ(:,:,i) = [dVdR(:,i), dVdR_dot(:,i), dVdtheta(:,i), dVdtheta_dot(:,i)]*...
+ P_t(7:10,7:10,i)*[dVdR(:,i), dVdR_dot(:,i), dVdtheta(:,i), dVdtheta_dot(:,i)]';
+
+ end
+
+% Correlation covariance MMX-Phobos
+ P_PhXYZ_MMX = zeros(3,3,size(X_t(9,:),2));
+ P_VPhXYZ_VMMX = zeros(3,3,size(X_t(9,:),2));
+
+ for i=1:size(X_t(9,:),2)
+ P_PhXYZ_MMX(:,:,i) = [dXdR(:,i),dXdtheta(:,i)]*[P_t(1:3,7,i), P_t(1:3,9,i)]';
+ P_VPhXYZ_VMMX(:,:,i) = [dVdR(:,i), dVdR_dot(:,i), dVdtheta(:,i), dVdtheta_dot(:,i)]*P_t(4:6,7:10,i)';
+
+ end
+
+% Relative distance covariance between MMX and PHobos
+ P_rel = P_PhXYZ + P_t(1:3,1:3,:) - 2*P_PhXYZ_MMX;
+
+% Relative velocity covariance between MMX and PHobos
+ P_Vrel = P_PhVXYZ + P_t(4:6,4:6,:) - 2*P_VPhXYZ_VMMX;
+
+
+end
\ No newline at end of file
diff --git a/Results.m b/Results.m
index f22b8a303b9c4180fc91ada3e759ddc578999708..3fb72386e37685361a7fef71c56e2f2bd564a1c4 100644
--- a/Results.m
+++ b/Results.m
@@ -186,6 +186,81 @@ function Results(Est, YObs_Full, pars, units)
xlabel('time [hour]')
ylabel('$\dot{\Phi}_{Ph}$ [rad/s]')
+
+% Phobos cartesian covaraiance
+
+ [P_PhXYZ, ~, P_rel_MMXPh, P_PhVXYZ, P_Vrel_MMXPh] = Phobos_cartesian_covariance(Est.X_t, Est.P_t, pars, units);
+
+ sqx_Ph = squeeze(P_PhXYZ(1,1,1:end-1));
+ sqy_Ph = squeeze(P_PhXYZ(2,2,1:end-1));
+ sqz_Ph = squeeze(P_PhXYZ(3,3,1:end-1));
+ SQRT_X_Ph = 3*sqrt(sqx_Ph + sqy_Ph + sqz_Ph);
+
+ sqx_VPh = squeeze(P_PhVXYZ(1,1,1:end-1));
+ sqy_VPh = squeeze(P_PhVXYZ(2,2,1:end-1));
+ sqz_VPh = squeeze(P_PhVXYZ(3,3,1:end-1));
+ SQRT_V_Ph = 3*sqrt(sqx_VPh + sqy_VPh + sqz_VPh);
+
+ sqx_rel = squeeze(P_rel_MMXPh(1,1,1:end-1));
+ sqy_rel = squeeze(P_rel_MMXPh(2,2,1:end-1));
+ sqz_rel = squeeze(P_rel_MMXPh(3,3,1:end-1));
+ SQRT_X_rel = 3*sqrt(sqx_rel + sqy_rel + sqz_rel);
+
+ sqx_Vrel = squeeze(P_Vrel_MMXPh(1,1,1:end-1));
+ sqy_Vrel = squeeze(P_Vrel_MMXPh(2,2,1:end-1));
+ sqz_Vrel = squeeze(P_Vrel_MMXPh(3,3,1:end-1));
+ SQRT_V_rel = 3*sqrt(sqx_Vrel + sqy_Vrel + sqz_Vrel);
+
+
+ figure()
+ subplot(2,2,1)
+ semilogy(t_obs(1:end-1)/3600,3.*squeeze(real(sqrt(P_PhXYZ(1,1,1:end-1)))),'Color','b','LineWidth',1)
+ hold on;
+ grid on;
+ semilogy(t_obs(1:end-1)/3600,3*squeeze(real(sqrt(P_PhXYZ(2,2,1:end-1)))),'Color','r','LineWidth',1)
+ semilogy(t_obs(1:end-1)/3600,3*squeeze(real(sqrt(P_PhXYZ(3,3,1:end-1)))),'Color','g','LineWidth',1)
+ semilogy(t_obs(1:end-1)/3600,SQRT_X_Ph,'Color','k','LineWidth',2)
+ xlabel('time $[hour]$')
+ ylabel('$[km]$')
+ title('Phobos position vector $3\sigma$ envelopes','Interpreter','latex','FontSize',16)
+ legend('$3\sigma_{x}$','$3\sigma_{y}$','$3\sigma_{z}$','$3 RMS$','Interpreter','latex','FontSize',14)
+ subplot(2,2,2)
+ semilogy(t_obs(1:end-1)/3600,3.*squeeze(real(sqrt(P_PhVXYZ(1,1,1:end-1)))),'Color','b','LineWidth',1)
+ hold on;
+ grid on;
+ semilogy(t_obs(1:end-1)/3600,3*squeeze(real(sqrt(P_PhVXYZ(2,2,1:end-1)))),'Color','r','LineWidth',1)
+ semilogy(t_obs(1:end-1)/3600,3*squeeze(real(sqrt(P_PhVXYZ(3,3,1:end-1)))),'Color','g','LineWidth',1)
+ semilogy(t_obs(1:end-1)/3600,SQRT_V_Ph,'Color','k','LineWidth',2)
+ xlabel('time $[hour]$')
+ ylabel('$[km]$')
+ title('Phobos velocity vector $3\sigma$ envelopes','Interpreter','latex','FontSize',16)
+ legend('$3\sigma_{\dot{x}}$','$3\sigma_{\dot{y}}$','$3\sigma_{\dot{z}}$','$3 RMS$','Interpreter','latex','FontSize',14)
+ subplot(2,2,3)
+ semilogy(t_obs(1:end-1)/3600,3.*squeeze(real(sqrt(P_rel_MMXPh(1,1,1:end-1)))),'Color','b','LineWidth',1)
+ hold on;
+ grid on;
+ semilogy(t_obs(1:end-1)/3600,3*squeeze(real(sqrt(P_rel_MMXPh(2,2,1:end-1)))),'Color','r','LineWidth',1)
+ semilogy(t_obs(1:end-1)/3600,3*squeeze(real(sqrt(P_rel_MMXPh(3,3,1:end-1)))),'Color','g','LineWidth',1)
+ semilogy(t_obs(1:end-1)/3600,SQRT_X_rel,'Color','k','LineWidth',2)
+ xlabel('time $[hour]$')
+ ylabel('$[km]$')
+ title('MMX-Phobos relative distance $3\sigma$ envelopes','Interpreter','latex','FontSize',16)
+ legend('$3\sigma_{x}$','$3\sigma_{y}$','$3\sigma_{z}$','$3 RMS$','Interpreter','latex','FontSize',14)
+ subplot(2,2,4)
+ semilogy(t_obs(1:end-1)/3600,3.*squeeze(real(sqrt(P_Vrel_MMXPh(1,1,1:end-1)))),'Color','b','LineWidth',1)
+ hold on;
+ grid on;
+ semilogy(t_obs(1:end-1)/3600,3*squeeze(real(sqrt(P_Vrel_MMXPh(2,2,1:end-1)))),'Color','r','LineWidth',1)
+ semilogy(t_obs(1:end-1)/3600,3*squeeze(real(sqrt(P_Vrel_MMXPh(3,3,1:end-1)))),'Color','g','LineWidth',1)
+ semilogy(t_obs(1:end-1)/3600,SQRT_V_rel,'Color','k','LineWidth',2)
+ xlabel('time $[hour]$')
+ ylabel('$[km]$')
+ title('MMX-Phobos relative velocity $3\sigma$ envelopes','Interpreter','latex','FontSize',16)
+ legend('$3\sigma_{x}$','$3\sigma_{y}$','$3\sigma_{z}$','$3 RMS$','Interpreter','latex','FontSize',14)
+
+
+
+
% Harmonics 3sigma Envelopes
% place0 = size(Est.X,1)-pars.nCoeff-pars.nBias;
diff --git a/YObs.mat b/YObs.mat
index 0c417514d7b3e3d30c4efd06924bbede9b577607..e8126f38fbd64912bfcbaf4b51b5000a6ef1fa09 100644
Binary files a/YObs.mat and b/YObs.mat differ
diff --git a/visible_features_model.m b/visible_features_model.m
index 61e08897de5d62de9da9c57248603f84260848d8..41bda2dc625adac0bea7d3081aee3bf46c432916 100644
--- a/visible_features_model.m
+++ b/visible_features_model.m
@@ -105,11 +105,11 @@ function [Y_LOS, Y_pix] = visible_features_model(MMX, Phobos, Sun, Phi, file_fea
z_cam = cross(v, I)/norm(cross(v, I));
y_cam = cross(z_cam, I)/norm(cross(z_cam, I));
- angle_feat = acos(dot(I,-i_feat));
- angle_sun = acos(dot(I,r_PhSun));
+ % angle_feat = acos(dot(I,-i_feat));
+ % angle_sun = acos(dot(I,r_PhSun));
- if (angle_feat<pars.FOV)&&(angle_sun>pars.FOV)&&(dot(r_sb_Ph, point_in_light(1:3,n))>0)
+ % if (angle_feat<pars.FOV)&&(angle_sun>pars.FOV)&&(dot(r_sb_Ph, point_in_light(1:3,n))>0)
k = 1;
visible = [candidate, visible];
LOS_feat = [y_cam'; z_cam'] * r_sf/norm(r_sf) + [random('Normal',0, pars.ObsNoise.camera);...
@@ -121,7 +121,7 @@ function [Y_LOS, Y_pix] = visible_features_model(MMX, Phobos, Sun, Phi, file_fea
pix_feat = pix_feat./repmat(pix_feat(3),3,1);
Y_pix = [[pix_feat(1:2); point_in_light(end,n)], Y_pix];
- end
+ % end
end