aggiunta covariance cartesian Phobos and relativa con MMX

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')

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;
 

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
 

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);
 
 
 

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;
 

Phobos_cartesian_covariance.m

+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

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;

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