Code owners
Assign users and groups as approvers for specific file changes. Learn more.
PlotFinali.m 27.15 KiB
clear
close all
clc
addpath('./ResultsHPC/');
addpath('../useful_functions/');
addpath(genpath('../mice/'));
addpath(genpath('../generic_kernels/'));
addpath('MMX_Fcn_CovarianceAnalyses/');
MMX_InitializeSPICE
cspice_furnsh(which('mar097.bsp'));
set(0,'DefaultTextInterpreter','latex');
set(0,'DefaultAxesFontSize', 16);
%% Confronto altezze
name10 = 'QSOH_PN10_alpha51_16.mat';
load(name10);
Est10 = Est;
name9 = 'QSOM_PN10_alpha1_16.mat';
load(name9);
Est9 = Est;
name8 = 'QSOLa_PN10_alpha71_16.mat';
load(name8);
Est8 = Est;
name7 = 'QSOLb_PN10_alpha781_16.mat';
load(name7);
Est7 = Est;
name6 = 'QSOLc_PN10_alpha1_16.mat';
load(name6);
Est6 = Est;
t_obs = Est6.t(1,:);
coeff = 0.8;
fine = round(coeff*size(t_obs,2));
fineCss = round(coeff*size(t_obs,2));
% Con Process noice
[SQRT_X6, SQRT_V6, P_t_C_20_6, P_t_C_22_6, mean6, std6, meanV6, stdV6,...
meanLib6, stdLib6, meanC206, stdC206, meanC226, stdC226] = Envelopes_Plot(Est6, fine, fineCss);
[SQRT_X7, SQRT_V7, P_t_C_20_7, P_t_C_22_7, mean7, std7, meanV7, stdV7,...
meanLib7, stdLib7, meanC207, stdC207, meanC227, stdC227] = Envelopes_Plot(Est7, fine, fineCss);
[SQRT_X8, SQRT_V8, P_t_C_20_8, P_t_C_22_8, mean8, std8, meanV8, stdV8,...
meanLib8, stdLib8, meanC208, stdC208, meanC228, stdC228] = Envelopes_Plot(Est8, fine, fineCss);
[SQRT_X9, SQRT_V9, P_t_C_20_9, P_t_C_22_9, mean9, std9, meanV9, stdV9,...
meanLib9, stdLib9, meanC209, stdC209, meanC229, stdC229] = Envelopes_Plot(Est9, fine, fineCss);
[SQRT_X10, SQRT_V10, P_t_C_20_10, P_t_C_22_10, mean10, std10, meanV10,...
stdV10, meanLib10, stdLib10, meanC2010, stdC2010, meanC2210, stdC2210] = Envelopes_Plot(Est10, fine, fineCss);
% Tradeoff
x = [27; 31; 49; 94; 198];
y = [mean6; mean7; mean8; mean9; mean10];
ypos = [std6; std7; std8; std9; std10];
yneg = [std6; std7; std8; std9; std10];
yV = [meanV6; meanV7; meanV8; meanV9; meanV10];
yposV = [stdV6; stdV7; stdV8; stdV9; stdV10];
ynegV = [stdV6; stdV7; stdV8; stdV9; stdV10];
xC20 = [27; 31; 49; 94; 198];
yC20 = [meanC206; meanC207; meanC208; meanC209; meanC2010];
yposC20 = [stdC206; stdC207; stdC208; stdC209; stdC2010];
ynegC20 = [stdC206; stdC207; stdC208; stdC209; stdC2010];
yC22 = [meanC226; meanC227; meanC228; meanC229; meanC2210];
yposC22 = [stdC226; stdC227; stdC228; stdC229; stdC2210];
ynegC22 = [stdC226; stdC227; stdC228; stdC229; stdC2210];
labels = {'$QSO-Lc$','$QSO-Lb$','$QSO-La$','$QSO-M$','$QSO-H$'};
% Same thing with the libration angle
xLib = [27; 31; 49; 94; 198];
yLib = [meanLib6; meanLib7; meanLib8; meanLib9; meanLib10].*180/pi;
yposLib = [stdLib6; stdLib7; stdLib8; stdLib9; stdLib10].*180/pi;
ynegLib = [stdLib6; stdLib7; stdLib8; stdLib9; stdLib10].*180/pi;
labelsLib = {'$QSO-Lc$','$QSO-Lb$','$QSO-La$','$QSO-M$','$QSO-H$'};
PlotComparison(Est10, Est9, Est8, Est7, Est6, fine, fineCss);
%% Senza Lidar
name5 = 'QSOH_PN10_alpha1_16_NOLidar.mat';
load(name5);
Est5 = Est;
name4 = 'QSOM_PN10_alpha1_16_NOLidar.mat';
load(name4);
Est4 = Est;
name3 = 'QSOLa_PN10_alpha1_16_NOLidar.mat';
load(name3);
Est3 = Est;
name2 = 'QSOLb_PN10_alpha1_16_NOLidar.mat';
load(name2);
Est2 = Est;
name1 = 'QSOLc_PN10_alpha1_16_NOLidar.mat';
load(name1);
Est1 = Est;
t_obs1 = Est1.t(1,:);
fine1 = round(coeff*size(t_obs1,2));
fineCss1 = round(coeff*size(t_obs1,2));
[SQRT_X1, SQRT_V1, P_t_C_20_1, P_t_C_22_1, mean1, std1, meanV1, stdV1,...
meanLib1, stdLib1, meanC201, stdC201, meanC221, stdC221] = Envelopes_Plot(Est1, fine1, fineCss1);
[SQRT_X2, SQRT_V2, P_t_C_20_2, P_t_C_22_2, mean2, std2, meanV2, stdV2,...
meanLib2, stdLib2, meanC202, stdC202, meanC222, stdC222] = Envelopes_Plot(Est2, fine1, fineCss1);
[SQRT_X3, SQRT_V3, P_t_C_20_3, P_t_C_22_3, mean3, std3, meanV3, stdV3,...
meanLib3, stdLib3, meanC203, stdC203, meanC223, stdC223] = Envelopes_Plot(Est3, fine1, fineCss1);
[SQRT_X4, SQRT_V4, P_t_C_20_4, P_t_C_22_4, mean4, std4, meanV4, stdV4,...
meanLib4, stdLib4, meanC204, stdC204, meanC224, stdC224] = Envelopes_Plot(Est4, fine1, fineCss1);
[SQRT_X5, SQRT_V5, P_t_C_20_5, P_t_C_22_5, mean5, std5, meanV5,...
stdV5, meanLib5, stdLib5, meanC205, stdC205, meanC225, stdC225] = Envelopes_Plot(Est5, fine1, fineCss1);
xNo_lidar = [27; 31; 49; 94; 198];
yNo_lidar = [mean1; mean2; mean3; mean4; mean5];
yposNo_lidar = [std1; std2; std3; std4; std5];
ynegNo_lidar = [std1; std2; std3; std4; std5];
yVNo_lidar = [meanV1; meanV2; meanV3; meanV4; meanV5];
yposVNo_lidar = [stdV1; stdV2; stdV3; stdV4; stdV5];
ynegVNo_lidar = [stdV1; stdV2; stdV3; stdV4; stdV5];
xC20No_lidar = [27; 31; 49; 94; 198];
yC20No_lidar = [meanC201; meanC202; meanC203; meanC204; meanC205];
yposC20No_lidar = [stdC201; stdC202; stdC203; stdC204; stdC205];
ynegC20No_lidar = [stdC201; stdC202; stdC203; stdC204; stdC205];
yC22No_lidar = [meanC221; meanC222; meanC223; meanC224; meanC225];
yposC22No_lidar = [stdC221; stdC222; stdC223; stdC224; stdC225];
ynegC22No_lidar = [stdC221; stdC222; stdC223; stdC224; stdC225];
% Same thing with the libration angle
xLib_NoLidar = [27; 31; 49; 94; 198];
yLib_NoLidar = [mean1; mean2; mean3; mean4; mean5].*180/pi;
yposLib_NoLidar = [std1; std2; std3; std4; std5].*180/pi;
ynegLib_NoLidar = [std1; std2; std3; std4; std5].*180/pi;
% PlotComparison(Est5, Est4, Est3, Est2, Est1, fine1, fineCss1)
%% Solo DSN
name15 = 'QSOH_PN10_alpha1_16_DSN.mat';
load(name15);
Est15 = Est;
name14 = 'QSOM_PN10_alpha1_16_DSN.mat';
load(name14);
Est14 = Est;
name13 = 'QSOLa_PN10_alpha1_16_DSN.mat';
load(name13);
Est13 = Est;
name12 = 'QSOLb_PN10_alpha1_16_DSN.mat';
load(name12);
Est12 = Est;
name11 = 'QSOLc_PN10_alpha1_16_DSN.mat';
load(name11);
Est11 = Est;
[SQRT_X11, SQRT_V11, P_t_C_20_11, P_t_C_22_11, mean11, std11, meanV11, stdV11,...
meanLib11, stdLib11, meanC2011, stdC2011, meanC2211, stdC2211] = Envelopes_Plot(Est11, fine1, fineCss1);
[SQRT_X12, SQRT_V12, P_t_C_20_12, P_t_C_22_12, mean12, std12, meanV12, stdV12,...
meanLib12, stdLib12, meanC2012, stdC2012, meanC2212, stdC2212] = Envelopes_Plot(Est12, fine1, fineCss1);
[SQRT_X13, SQRT_V13, P_t_C_20_13, P_t_C_22_13, mean13, std13, meanV13, stdV13,...
meanLib13, stdLib13, meanC2013, stdC2013, meanC2213, stdC2213] = Envelopes_Plot(Est13, fine1, fineCss1);
[SQRT_X14, SQRT_V14, P_t_C_20_14, P_t_C_22_14, mean14, std14, meanV14, stdV14,...
meanLib14, stdLib14, meanC2014, stdC2014, meanC2214, stdC2214] = Envelopes_Plot(Est14, fine1, fineCss1);
[SQRT_X15, SQRT_V15, P_t_C_20_15, P_t_C_22_15, mean15, std15, meanV15,...
stdV15, meanLib15, stdLib15, meanC2015, stdC2015, meanC2215, stdC2215] = Envelopes_Plot(Est15, fine1, fineCss1);
xDSN = [27; 31; 49; 94; 198];
yDSN = [mean11; mean12; mean13; mean14; mean15];
yposDSN = [std11; std12; std13; std14; std15];
ynegDSN = [std11; std12; std13; std14; std15];
yVDSN = [meanV11; meanV12; meanV13; meanV14; meanV15];
yposVDSN = [stdV11; stdV12; stdV13; stdV14; stdV15];
ynegVDSN = [stdV11; stdV12; stdV13; stdV14; stdV15];
xC20DSN = [27; 31; 49; 94; 198];
yC20DSN = [meanC2011; meanC2012; meanC2013; meanC2014; meanC2015];
yposC20DSN = [stdC2011; stdC2012; stdC2013; stdC2014; stdC2015];
ynegC20DSN = [stdC2011; stdC2012; stdC2013; stdC2014; stdC2015];
yC22DSN = [meanC2211; meanC2212; meanC2213; meanC2214; meanC2215];
yposC22DSN = [stdC2211; stdC2212; stdC2213; stdC2214; stdC2215];
ynegC22DSN = [stdC2211; stdC2212; stdC2213; stdC2214; stdC2215];
% Same thing with the libration angle
xLib_DSN = [27; 31; 49; 94; 198];
yLib_DSN = [mean11; mean12; mean13; mean14; mean15].*180/pi;
yposLib_DSN = [std11; std12; std13; std14; std15].*180/pi;
ynegLib_DSN = [std11; std12; std13; std14; std15].*180/pi;
% PlotComparison(Est15, Est14, Est13, Est12, Est11, fine1, fineCss1)
%% No range
name20 = 'QSOH_PN10_alpha71_16_NORange.mat';
load(name20);
Est20 = Est;
name19 = 'QSOM_PN10_alpha1_16_NORange.mat';
load(name19);
Est19 = Est;
name18 = 'QSOLa_PN10_alpha71_16_NORange.mat';
load(name18);
Est18 = Est;
name17 = 'QSOLb_PN10_alpha1_16_NORange.mat';
load(name17);
Est17 = Est;
name16 = 'QSOLc_PN10_alpha1_16_NORange.mat';
load(name16);
Est16 = Est;
[SQRT_X16, SQRT_V16, P_t_C_20_16, P_t_C_22_16, mean16, std16, meanV16, stdV16,...
meanLib16, stdLib16, meanC2016, stdC2016, meanC2216, stdC2216] = Envelopes_Plot(Est16, fine, fineCss);
[SQRT_X17, SQRT_V17, P_t_C_20_17, P_t_C_22_17, mean17, std17, meanV17, stdV17,...
meanLib17, stdLib17, meanC2017, stdC2017, meanC2217, stdC2217] = Envelopes_Plot(Est17, fine, fineCss);
[SQRT_X18, SQRT_V18, P_t_C_20_18, P_t_C_22_18, mean18, std18, meanV18, stdV18,...
meanLib18, stdLib18, meanC2018, stdC2018, meanC2218, stdC2218] = Envelopes_Plot(Est18, fine, fineCss);
[SQRT_X19, SQRT_V19, P_t_C_20_19, P_t_C_22_19, mean19, std19, meanV19, stdV19,...
meanLib19, stdLib19, meanC2019, stdC2019, meanC2219, stdC2219] = Envelopes_Plot(Est19, fine, fineCss);
[SQRT_X20, SQRT_V20, P_t_C_20_20, P_t_C_22_20, mean20, std20, meanV20,...
stdV20, meanLib20, stdLib20, meanC2020, stdC2020, meanC2220, stdC2220] = Envelopes_Plot(Est20, fine, fineCss);
xNorho = [27; 31; 49; 94; 198];
yNorho = [mean16; mean17; mean18; mean19; mean20];
yposNorho = [std16; std17; std18; std19; std20];
ynegNorho = [std16; std17; std18; std19; std20];
yVNorho = [meanV16; meanV17; meanV18; meanV19; meanV20];
yposVNorho = [stdV16; stdV17; stdV18; stdV19; stdV20];
ynegVNorho = [stdV16; stdV17; stdV18; stdV19; stdV20];
xC20Norho = [27; 31; 49; 94; 198];
yC20Norho = [meanC2016; meanC2017; meanC208; meanC2019; meanC2020];
yposC20Norho = [stdC2016; stdC2017; stdC2018; stdC2019; stdC2020];
ynegC20Norho = [stdC2016; stdC2017; stdC2018; stdC2019; stdC2020];
yC22Norho = [meanC2216; meanC2217; meanC2218; meanC2219; meanC2220];
yposC22Norho = [stdC2216; stdC2217; stdC2218; stdC2219; stdC2220];
ynegC22Norho = [stdC2216; stdC2217; stdC2218; stdC2219; stdC2220];
% Libration
xLib_Norho = [27; 31; 49; 94; 198];
yLib_Norho = [meanLib16; meanLib17; meanLib18; meanLib19; meanLib20].*180/pi;
yposLib_Norho = [stdLib16; stdLib17; stdLib18; stdLib19; stdLib20].*180/pi;
ynegLib_Norho = [stdLib16; stdLib17; stdLib18; stdLib19; stdLib20].*180/pi;
% PlotComparison(Est20, Est19, Est18, Est17, Est16, fine, fineCss)
%% No range rate
name25 = 'QSOH_PN10_alpha1_16_NORate.mat';
load(name25);
Est25 = Est;
name24 = 'QSOM_PN10_alpha1_16_NORate.mat';
load(name24);
Est24 = Est;
name23 = 'QSOLa_PN10_alpha1_16_NORate.mat';
load(name23);Est23 = Est;
name22 = 'QSOLb_PN10_alpha41_16_NORate.mat';
load(name22);
Est22 = Est;
name21 = 'QSOLc_PN10_alpha1_16_NORate.mat';
load(name21);
Est21 = Est;
[SQRT_X21, SQRT_V21, P_t_C_20_21, P_t_C_22_21, mean21, std21, meanV21, stdV21,...
meanLib21, stdLib21, meanC2021, stdC2021, meanC2221, stdC2221] = Envelopes_Plot(Est21, fine, fineCss);
[SQRT_X22, SQRT_V22, P_t_C_20_22, P_t_C_22_22, mean22, std22, meanV22, stdV22,...
meanLib22, stdLib22, meanC2022, stdC2022, meanC2222, stdC2222] = Envelopes_Plot(Est22, fine, fineCss);
[SQRT_X23, SQRT_V23, P_t_C_20_23, P_t_C_22_23, mean23, std23, meanV23, stdV23,...
meanLib23, stdLib23, meanC2023, stdC2023, meanC2223, stdC2223] = Envelopes_Plot(Est23, fine, fineCss);
[SQRT_X24, SQRT_V24, P_t_C_20_24, P_t_C_22_24, mean24, std24, meanV24, stdV24,...
meanLib24, stdLib24, meanC2024, stdC2024, meanC2224, stdC2224] = Envelopes_Plot(Est24, fine, fineCss);
[SQRT_X25, SQRT_V25, P_t_C_20_25, P_t_C_22_25, mean25, std25, meanV25,...
stdV25, meanLib25, stdLib25, meanC2025, stdC2025, meanC2225, stdC2225] = Envelopes_Plot(Est25, fine, fineCss);
xNorate = [27; 31; 49; 94; 198];
yNorate = [mean21,; mean22; mean23; mean24; mean25];
yposNorate = [std21,; std22; std23; std24; std25];
ynegNorate = [std21,; std22; std23; std24; std25];
yVNorate = [meanV21,; meanV22; meanV23; meanV24; meanV25];
yposVNorate = [stdV21,; stdV22; stdV23; stdV24; stdV25];
ynegVNorate = [stdV21,; stdV22; stdV23; stdV24; stdV25];
xC20Norate = [27; 31; 49; 94; 198];
yC20Norate = [meanC2021,; meanC2022; meanC2023; meanC2024; meanC2025];
yposC20Norate = [stdC2021,; stdC2022; stdC2023; stdC2024; stdC2025];
ynegC20Norate = [stdC2021,; stdC2022; stdC2023; stdC2024; stdC2025];
yC22Norate = [meanC2221,; meanC2222; meanC2223; meanC2224; meanC2225];
yposC22Norate = [stdC2221,; stdC2222; stdC2223; stdC2224; stdC2225];
ynegC22Norate = [stdC2221,; stdC2222; stdC2223; stdC2224; stdC2225];
% Libration
xLib_Norate = [27; 31; 49; 94; 198];
yLib_Norate = [meanLib21; meanLib22; meanLib23; meanLib24; meanLib25].*180/pi;
yposLib_Norate = [stdLib21; stdLib22; stdLib23; stdLib24; stdLib25].*180/pi;
ynegLib_Norate = [stdLib21; stdLib22; stdLib23; stdLib24; stdLib25].*180/pi;
% PlotComparison(Est25, Est24, Est23, Est22, Est21, fine, fineCss)
%% No LOS
name30 = 'QSOH_PN10_alpha51_16_NOLOS.mat';
load(name30);
Est30 = Est;
name29 = 'QSOM_PN10_alpha1_16_NOLOS.mat';
load(name29);
Est29 = Est;
name28 = 'QSOLa_PN10_alpha71_16_NOLOS.mat';
load(name28);
Est28 = Est;
name27 = 'QSOLb_PN10_alpha771_16_NOLOS.mat';
load(name27);
Est27 = Est;
name26 = 'QSOLc_PN10_alpha551_16_NOLOS.mat';
load(name26);
Est26 = Est;
[SQRT_X26, SQRT_V26, P_t_C_20_26, P_t_C_22_26, mean26, std26, meanV26, stdV26,...
meanLib26, stdLib26, meanC2026, stdC2026, meanC2226, stdC2226] = Envelopes_Plot(Est26, fine, fineCss);[SQRT_X27, SQRT_V27, P_t_C_20_27, P_t_C_22_27, mean27, std27, meanV27, stdV27,...
meanLib27, stdLib27, meanC2027, stdC2027, meanC2227, stdC2227] = Envelopes_Plot(Est27, fine, fineCss);
[SQRT_X28, SQRT_V28, P_t_C_20_28, P_t_C_22_28, mean28, std28, meanV28, stdV28,...
meanLib28, stdLib28, meanC2028, stdC2028, meanC2228, stdC2228] = Envelopes_Plot(Est28, fine, fineCss);
[SQRT_X29, SQRT_V29, P_t_C_20_29, P_t_C_22_29, mean29, std29, meanV29, stdV29,...
meanLib29, stdLib29, meanC2029, stdC2029, meanC2229, stdC2229] = Envelopes_Plot(Est29, fine, fineCss);
[SQRT_X30, SQRT_V30, P_t_C_20_30, P_t_C_22_30, mean30, std30, meanV30,...
stdV30, meanLib30, stdLib30, meanC2030, stdC2030, meanC2230, stdC2230] = Envelopes_Plot(Est30, fine-300, fineCss-300);
% MMX state vector RMS
xNoLOS = [27; 31; 49; 94; 198];
yNoLOS = [mean26; mean27; mean8; mean29; mean30];
yposNoLOS = [std26; std27; std28; std29; std30];
ynegNoLOS = [std26; std27; std28; std29; std30];
yVNoLOS = [meanV26; meanV27; meanV8; meanV29; meanV30];
yposVNoLOS = [stdV26; stdV27; stdV28; stdV29; stdV30];
ynegVNoLOS = [stdV26; stdV27; stdV28; stdV29; stdV30];
% Css
xC20NoLOS = [27; 31; 49; 94; 198];
yC20NoLOS = [meanC206; meanC207; meanC2028; meanC2029; meanC2030];
yposC20NoLOS = [stdC2026; stdC2027; stdC2028; stdC2029; stdC2030];
ynegC20NoLOS = [stdC2026; stdC2027; stdC2028; stdC2029; stdC2030];
yC22NoLOS = [meanC2216; meanC2227; meanC2228; meanC2229; meanC2230];
yposC22NoLOS = [stdC2226; stdC2227; stdC2228; stdC2229; stdC2230];
ynegC22NoLOS = [stdC2226; stdC2227; stdC2228; stdC2229; stdC2230];
% Libration
xLib_NoLOS = [27; 31; 49; 94; 198];
yLib_NoLOS = [meanLib26; meanLib27; meanLib28; meanLib29; meanLib30].*180/pi;
yposLib_NoLOS = [stdLib26; stdLib27; stdLib28; stdLib29; stdLib30].*180/pi;
ynegLib_NoLOS = [stdLib26; stdLib27; stdLib28; stdLib29; stdLib30].*180/pi;
% PlotComparison(Est30, Est29, Est28, Est27, Est26, fine, fineCss)
%% Plots
figure()
subplot(1,2,1)
errorbar(x,y,yneg,ypos,'-s','LineWidth',2.2,'MarkerSize',5)
hold on
errorbar(xNorho,yNorho,ynegNorho,yposNorho,'-s',...
'LineWidth',1.2,'MarkerSize',5)
errorbar(xNorate,yNorate,ynegNorate,yposNorate,'-s',...
'LineWidth',1.2,'MarkerSize',5)
errorbar(xNo_lidar,yNo_lidar,ynegNo_lidar,yposNo_lidar,'-s',...
'LineWidth',1.2,'MarkerSize',5)
errorbar(xNoLOS,yNoLOS,ynegNoLOS,yposNoLOS,'-s',...
'LineWidth',1.2,'MarkerSize',5)
errorbar(xDSN,yDSN,ynegDSN,yposDSN,'-s',...
'LineWidth',1.2,'MarkerSize',5)
set(gca,'YScale','log')
set(gca,'TickLabelInterpreter','latex')
set(gca,'XScale','log');
set(gca,'XTick',[27; 31; 49; 94; 198]);
set(gca,'XTickLabel',labels,'Fontsize',26);
set(gca,'XTickLabelRotation',45)
legend('Full','No $\rho$','No $\dot{\rho}$','No Lidar','No LOS','DSN only','Interpreter','latex','Fontsize',26)
ylabel('$[km]$','Fontsize',26)
xlim([2e1,4e2])
grid on
subplot(1,2,2)
errorbar(x,yV,ynegV,yposV,'-s','LineWidth',2.2,'MarkerSize',5)
hold on
errorbar(xNorho,yVNorho,ynegVNorho,yposVNorho,'-s',... 'LineWidth',1.2,'MarkerSize',5)
errorbar(xNorate,yVNorate,ynegVNorate,yposVNorate,'-s',...
'LineWidth',1.2,'MarkerSize',5)
errorbar(xNo_lidar,yVNo_lidar,ynegVNo_lidar,yposVNo_lidar,'-s',...
'LineWidth',1.2,'MarkerSize',5)
errorbar(xNoLOS,yVNoLOS,ynegVNoLOS,yposVNoLOS,'-s',...
'LineWidth',1.2,'MarkerSize',5)
errorbar(xDSN,yVDSN,ynegVDSN,yposVDSN,'-s',...
'LineWidth',1.2,'MarkerSize',5)
grid on
set(gca,'YScale','log')
set(gca,'TickLabelInterpreter','latex')
set(gca,'XScale','log')
set(gca,'XTick',[27; 31; 49; 94; 198]);
set(gca,'XTickLabel',labels,'Fontsize',26);
set(gca,'XTickLabelRotation',45)
legend('Full','No $\rho$','No $\dot{\rho}$','No Lidar','No LOS','DSN only','Interpreter','latex','Fontsize',26)
ylabel('$[km/s]$','Fontsize',26)
xlim([2e1,4e2])
% Css
figure()
subplot(1,2,1)
errorbar(xC20,yC20,ynegC20,yposC20,'-s','LineWidth',2.2,'MarkerSize',5)
hold on
errorbar(xC20Norho,yC20Norho,ynegC20Norho,yposC20Norho,'-s',...
'LineWidth',1.2,'MarkerSize',5)
errorbar(xC20Norate,yC20Norate,ynegC20Norate,yposC20Norate,'-s',...
'LineWidth',1.2,'MarkerSize',5)
errorbar(xC20No_lidar,yC20No_lidar,ynegC20No_lidar,yposC20No_lidar,'-s',...
'LineWidth',1.2,'MarkerSize',5)
errorbar(xC20NoLOS,yC20NoLOS,ynegC20NoLOS,yposC20NoLOS,'-s',...
'LineWidth',1.2,'MarkerSize',5)
errorbar(xC20DSN,yC20DSN,ynegC20DSN,yposC20DSN,'-s',...
'LineWidth',1.2,'MarkerSize',5)
set(gca,'YScale','log')
set(gca,'TickLabelInterpreter','latex')
set(gca,'XScale','log')
set(gca,'XTick',[27; 31; 49; 94; 198]);
set(gca,'XTickLabel',labels,'Fontsize',26);
set(gca,'XTickLabelRotation',45)
legend('Full','No $\rho$','No $\dot{\rho}$','No Lidar','No LOS','DSN only','Interpreter','latex','Fontsize',26)
ylabel('$3\sigma_{C_{20}}$','Fontsize',30)
xlim([2e1,4e2])
grid on
subplot(1,2,2)
errorbar(x,yC22,ynegC22,yposC22,'-s','LineWidth',2.2,'MarkerSize',5)
hold on
errorbar(xNorho,yC22Norho,ynegC22Norho,yposC22Norho,'-s',...
'LineWidth',1.2,'MarkerSize',5)
errorbar(xNorate,yC22Norate,ynegC22Norate,yposC22Norate,'-s',...
'LineWidth',1.2,'MarkerSize',5)
errorbar(xNo_lidar,yC22No_lidar,ynegC22No_lidar,yposC22No_lidar,'-s',...
'LineWidth',1.2,'MarkerSize',5)
errorbar(xNoLOS,yC22NoLOS,ynegC22NoLOS,yposC22NoLOS,'-s',...
'LineWidth',1.2,'MarkerSize',5)
errorbar(xDSN,yC22DSN,ynegC22DSN,yposC22DSN,'-s',...
'LineWidth',1.2,'MarkerSize',5)
grid on
set(gca,'YScale','log')
set(gca,'TickLabelInterpreter','latex')
set(gca,'XScale','log')
set(gca,'XTick',[27; 31; 49; 94; 198]);
set(gca,'XTickLabel',labels,'Fontsize',26);
set(gca,'XTickLabelRotation',45)
legend('Full','No $\rho$','No $\dot{\rho}$','No Lidar','No LOS','DSN only','Interpreter','latex','Fontsize',26)
ylabel('$3\sigma_{C_{22}}$','Fontsize',30)
xlim([2e1,4e2])
% Libration
figure()
errorbar(xLib,yLib,ynegLib,yposLib,'-s','LineWidth',2.5,'MarkerSize',5)
hold on
errorbar(xLib_Norho,yLib_Norho,ynegLib_Norho,yposLib_Norho,'-s',...
'LineWidth',1.5,'MarkerSize',5)
errorbar(xLib_Norate,yLib_Norate,ynegLib_Norate,yposLib_Norate,'-s',...
'LineWidth',1.5,'MarkerSize',5)
errorbar(xLib_NoLidar,yLib_NoLidar,ynegLib_NoLidar,yposLib_NoLidar,'-s',...
'LineWidth',1.5,'MarkerSize',5)
errorbar(xLib_NoLOS,yLib_NoLOS,ynegLib_NoLOS,yposLib_NoLOS,'-s',...
'LineWidth',1.5,'MarkerSize',5)
errorbar(xLib_DSN,yLib_DSN,ynegLib_DSN,yposLib_DSN,'-s',...
'LineWidth',1.5,'MarkerSize',5)
set(gca,'YScale','log')
set(gca,'TickLabelInterpreter','latex')
set(gca,'XScale','log')
set(gca,'XTick',[27; 31; 49; 94; 198]);
set(gca,'XTickLabel',labelsLib,'Fontsize',30);
set(gca,'XTickLabelRotation',45)
ylabel('$[deg]$','Fontsize',30)
legend('Full','No $\rho$','No $\dot{\rho}$','No Lidar','No LOS','DSN only','Interpreter','latex','Fontsize',30)
xlim([2e1,4e2])
grid on
%% Confronto geometrie sulla Lb
% load('QSOM_PN10_alpha1_16.mat');
% EstQSO = Est;
% load('3DQSOM_PN10_alpha91_16.mat');
% Est3D = Est;
% load('SwingQSOM_PN10_alpha91_16.mat');
% EstSwing = Est;
% load('QSOLa_PN10_alpha71_16.mat');
% EstQSO = Est;
% load('3DQSOLa_PN10_alpha71_16.mat');
% Est3D = Est;
% load('SwingQSOLa_PN10_alpha71_16.mat');
% EstSwing = Est;
load('QSOLb_PN10_alpha61_16.mat');
EstQSO = Est;
load('3DQSOLb_PN10_alpha61_16.mat');
Est3D = Est;
load('SwingQSOLb_PN10_alpha61_16.mat');
EstSwing = Est;
coeff = .99;
t_obs = EstQSO.t(1,:);
fine = round(coeff*size(t_obs,2));
fineCss = round(coeff*size(t_obs,2));
[par,~] = MMX_DefineNewModelParametersAndUnits;
C_20 = par.SH.Clm(3,1);
C_22 = par.SH.Clm(3,3);
t_obs = EstQSO.t(1,:);
[SQRT_X_QSO, SQRT_V_QSO, P_t_C_20_QSO, P_t_C_22_QSO, meanQSO, stdQSO, meanVQSO, stdVQSO,...
meanLibQSO, stdLibQSO,] = Envelopes_Plot(EstQSO, fine, fineCss);
[SQRT_X_3D, SQRT_V_3D, P_t_C_20_3D, P_t_C_22_3D, mean3D, std3D, meanV3D, stdV3D,...
meanLib3D, stdLib3D] = Envelopes_Plot(Est3D, fine, fineCss);[SQRT_X_Swing, SQRT_V_Swing, P_t_C_20_Swing, P_t_C_22_Swing, meanSwing, stdSwing, meanVSwing, stdVSwing,...
meanLibSwing, stdLibSwing] = Envelopes_Plot(EstSwing, fine, fineCss);
C_20_QSO_perc= abs(3*real(sqrt(P_t_C_20_QSO))/C_20);
C_22_QSO_perc= abs(3*real(sqrt(P_t_C_22_QSO))/C_22);
C_20_3D_perc= abs(3*real(sqrt(P_t_C_20_3D))/C_20);
C_22_3D_perc= abs(3*real(sqrt(P_t_C_22_3D))/C_22);
C_20_Swing_perc= abs(3*real(sqrt(P_t_C_20_Swing))/C_20);
C_22_Swing_perc= abs(3*real(sqrt(P_t_C_22_Swing))/C_22);
figure()
subplot(1,2,1)
semilogy(t_obs(1:end-1)/3600,3.*real(sqrt(P_t_C_20_QSO)),'LineWidth',1.2);
grid on;
hold on;
semilogy(t_obs(1:end-1)/3600,3.*real(sqrt(P_t_C_20_3D)),'LineWidth',1.2);
semilogy(t_obs(1:end-1)/3600,3.*real(sqrt(P_t_C_20_Swing)),'LineWidth',1.2);
xlabel('Time $[hour]$','Fontsize',26)
ylabel('$3\sigma_{C_{20}}$','Fontsize',30)
legend('QSO','3D-QSO','Swing-QSO','Interpreter','latex','Fontsize',26)
subplot(1,2,2)
semilogy(t_obs(1:end-1)/3600,3.*real(sqrt(P_t_C_22_QSO)),'LineWidth',1.2);
grid on;
hold on;
semilogy(t_obs(1:end-1)/3600,3.*real(sqrt(P_t_C_22_3D)),'LineWidth',1.2);
semilogy(t_obs(1:end-1)/3600,3.*real(sqrt(P_t_C_22_Swing)),'LineWidth',1.2);
xlabel('Time $[hour]$','Fontsize',26)
ylabel('$3\sigma_{C_{22}}$','Fontsize',30)
legend('QSO','3D-QSO','Swing-QSO','Interpreter','latex','Fontsize',26)
% figure()
% subplot(1,2,1)
% semilogy(t_obs(1:end-1)/3600,C_20_QSO_perc,'LineWidth',1.2);
% grid on;
% hold on;
% semilogy(t_obs(1:end-1)/3600,C_20_3D_perc,'LineWidth',1.2);
% semilogy(t_obs(1:end-1)/3600,C_20_Swing_perc,'LineWidth',1.2);
% xlabel('Time $[hour]$','Fontsize',26)
% ylabel('%','Fontsize',26)
% legend('QSO','3D-QSO','Swing-QSO','Interpreter','latex','Fontsize',26)
% subplot(1,2,2)
% semilogy(t_obs(1:end-1)/3600,C_22_QSO_perc,'LineWidth',1.2);
% grid on;
% hold on;
% semilogy(t_obs(1:end-1)/3600,C_22_3D_perc,'LineWidth',1.2);
% semilogy(t_obs(1:end-1)/3600,C_22_Swing_perc,'LineWidth',1.2);
% xlabel('Time $[hour]$','Fontsize',26)
% ylabel('%','Fontsize',26)
% legend('QSO','3D-QSO','Swing-QSO','Interpreter','latex','Fontsize',26)
figure()
semilogy(t_obs(1:end-1)/3600,180/pi*3.*squeeze(real(sqrt(EstQSO.P_t(13,13,1:end-1)))),'LineWidth',1.2)
grid on;
hold on;
semilogy(t_obs(1:end-1)/3600,180/pi*3.*squeeze(real(sqrt(Est3D.P_t(13,13,1:end-1)))),'LineWidth',1.2)
semilogy(t_obs(1:end-1)/3600,180/pi*3.*squeeze(real(sqrt(EstSwing.P_t(13,13,1:end-1)))),'LineWidth',1.2)
title('Libration angle $3\sigma$ envelopes','FontSize',26)
xlabel('Time [hour]','Interpreter','latex','Fontsize',26)
ylabel('$3\sigma_{\Phi_{Ph}}$ [deg]','Interpreter','latex','Fontsize',30)
legend('QSO','3D-QSO','Swing-QSO','Interpreter','latex','Fontsize',26)
% MMX Position and velocity
figure()
subplot(1,2,1)
semilogy(t_obs(1:end-1)/3600,SQRT_X_QSO)hold on;
grid on;
semilogy(t_obs(1:end-1)/3600,SQRT_X_3D)
semilogy(t_obs(1:end-1)/3600,SQRT_X_Swing)
xlabel('time $[hour]$')
ylabel('$[km]$')
title('MMX position vector $3\sigma$ envelopes','Interpreter','latex','FontSize',16)
legend('QSO','3D-QSO','Swing-QSO','Interpreter','latex','FontSize',14)
subplot(1,2,2)
semilogy(t_obs(1:end-1)/3600,SQRT_V_QSO)
hold on;
grid on;
semilogy(t_obs(1:end-1)/3600,SQRT_V_3D)
semilogy(t_obs(1:end-1)/3600,SQRT_V_Swing)
xlabel('time $[hour]$')
ylabel('$[km/s]$')
title('MMX velocity vector $3\sigma$ envelopes','Interpreter','latex','FontSize',16)
legend('QSO','3D-QSO','Swing-QSO','Interpreter','latex','FontSize',14)
% Phobos states
figure()
subplot(2,4,1)
semilogy(t_obs(1:end-1)/3600,3.*squeeze(real(sqrt(EstQSO.P_t(7,7,1:end-1)))))
grid on;
hold on;
semilogy(t_obs(1:end-1)/3600,3.*squeeze(real(sqrt(Est3D.P_t(7,7,1:end-1)))))
semilogy(t_obs(1:end-1)/3600,3.*squeeze(real(sqrt(EstSwing.P_t(7,7,1:end-1)))))
xlabel('time [hour]')
ylabel('$R_{Ph}$ [km]')
subplot(2,4,2)
semilogy(t_obs(1:end-1)/3600,180/pi*3.*squeeze(real(sqrt(EstQSO.P_t(9,9,1:end-1)))))
grid on;
hold on;
semilogy(t_obs(1:end-1)/3600,180/pi*3.*squeeze(real(sqrt(Est3D.P_t(9,9,1:end-1)))))
semilogy(t_obs(1:end-1)/3600,180/pi*3.*squeeze(real(sqrt(EstSwing.P_t(9,9,1:end-1)))))
xlabel('time [hour]')
ylabel('$\theta_{Ph}$ [deg]')
subplot(2,4,3)
semilogy(t_obs(1:end-1)/3600,180/pi*3.*squeeze(real(sqrt(EstQSO.P_t(11,11,1:end-1)))))
grid on;
hold on;
semilogy(t_obs(1:end-1)/3600,180/pi*3.*squeeze(real(sqrt(Est3D.P_t(11,11,1:end-1)))))
semilogy(t_obs(1:end-1)/3600,180/pi*3.*squeeze(real(sqrt(EstSwing.P_t(11,11,1:end-1)))))
xlabel('time [hour]')
ylabel('$\Phi_{M}$ [deg]')
subplot(2,4,4)
semilogy(t_obs(1:end-1)/3600,180/pi*3.*squeeze(real(sqrt(EstQSO.P_t(13,13,1:end-1)))))
grid on;
hold on;
semilogy(t_obs(1:end-1)/3600,180/pi*3.*squeeze(real(sqrt(Est3D.P_t(13,13,1:end-1)))))
semilogy(t_obs(1:end-1)/3600,180/pi*3.*squeeze(real(sqrt(EstSwing.P_t(13,13,1:end-1)))))
xlabel('time [hour]')
ylabel('$\Phi_{Ph}$ [deg]')
subplot(2,4,5)
semilogy(t_obs(1:end-1)/3600,3.*squeeze(real(sqrt(EstQSO.P_t(8,8,1:end-1)))))
grid on;
hold on;
semilogy(t_obs(1:end-1)/3600,3.*squeeze(real(sqrt(Est3D.P_t(8,8,1:end-1)))))
semilogy(t_obs(1:end-1)/3600,3.*squeeze(real(sqrt(EstSwing.P_t(8,8,1:end-1)))))
xlabel('time [hour]')
ylabel('$\dot{R}_{Ph}$ [km/s]')
subplot(2,4,6)
semilogy(t_obs(1:end-1)/3600,180/pi*3.*squeeze(real(sqrt(EstQSO.P_t(10,10,1:end-1)))))
grid on;
hold on;
semilogy(t_obs(1:end-1)/3600,180/pi*3.*squeeze(real(sqrt(Est3D.P_t(10,10,1:end-1)))))
semilogy(t_obs(1:end-1)/3600,180/pi*3.*squeeze(real(sqrt(EstSwing.P_t(10,10,1:end-1)))))
xlabel('time [hour]')ylabel('$\dot{\theta}_{Ph}$ [deg/s]')
subplot(2,4,7)
semilogy(t_obs(1:end-1)/3600,180/pi*3.*squeeze(real(sqrt(EstQSO.P_t(12,12,1:end-1)))))
grid on;
hold on;
semilogy(t_obs(1:end-1)/3600,180/pi*3.*squeeze(real(sqrt(Est3D.P_t(12,12,1:end-1)))))
semilogy(t_obs(1:end-1)/3600,180/pi*3.*squeeze(real(sqrt(EstSwing.P_t(12,12,1:end-1)))))
xlabel('time [hour]')
ylabel('$\dot{\Phi}_{M}$ [deg/s]')
subplot(2,4,8)
semilogy(t_obs(1:end-1)/3600,180/pi*3.*squeeze(real(sqrt(EstQSO.P_t(14,14,1:end-1)))))
grid on;
hold on;
semilogy(t_obs(1:end-1)/3600,180/pi*3.*squeeze(real(sqrt(Est3D.P_t(14,14,1:end-1)))))
semilogy(t_obs(1:end-1)/3600,180/pi*3.*squeeze(real(sqrt(EstSwing.P_t(14,14,1:end-1)))))
xlabel('time [hour]')
ylabel('$\dot{\Phi}_{Ph}$ [deg/s]')