diff --git a/Cov_AnalysisHPC.m b/Cov_AnalysisHPC.m index 59368b71fdbccfd8f2137a136b35b9e562372fa0..fa96d59cc8104f963861e238f85df764ae098944 100644 --- a/Cov_AnalysisHPC.m +++ b/Cov_AnalysisHPC.m @@ -16,7 +16,7 @@ clc addpath(genpath('../mice/')); addpath(genpath('../dynamical_model/')); addpath(genpath('../generic_kernels/')); - addpath(genpath('../paper_after_conf/MMX_Product/MMX_BSP_Files_GravLib/')); + addpath(genpath('./MMX_Product/MMX_BSP_Files_GravLib/')); MMX_InitializeSPICE cspice_furnsh('../generic_kernels/mar097.bsp'); cspice_furnsh('./Model_Check/MARPHOFRZ.txt'); @@ -25,10 +25,10 @@ clc %% QSO-H - cspice_furnsh('./MMX_Product/MMX_BSP_Files_GravLib/MMX_SwingQSO_031_011_2x2_826891269_828619269.bsp'); + cspice_furnsh('./MMX_Product/MMX_BSP_Files_GravLib/MMX_QSO_094_2x2_826891269_828619269.bsp'); cspice_furnsh('./MMX_Product/MMX_BSP_Files_GravLib/Phobos_826891269_828619269.bsp'); - load("./ObservationsHPC/Observations-SwingQSOLb_16.mat"); + load("./ObservationsHPC/Observations-QSOM_16_Lidar20Min.mat"); % Model parameters [par, units] = MMX_DefineNewModelParametersAndUnits; @@ -41,6 +41,7 @@ clc % Covariance analysis parameters [par, units] = New_MMX_CovarianceAnalysisParameters(par, units); + par.interval_lidar = 1200; par.sigma = 1e-10/(units.vsf*units.tsf); par.sigmaPh = 0/(units.vsf*units.tsf); @@ -59,19 +60,19 @@ clc Est0.P0 = par.P0; Est0.t0 = data*units.tsf; - par.alpha = 6.6e-1; + par.alpha = 1; - [Est] = UKF_CovarianceAnalysis(Est0,@Dynamics_MPHandMMX_Inertia,... - @Cov_Dynamics_Good, @New_Observables_model,... - par.R,YObs_Full,par,units); - - save('SwingQSOLb_PN10_alpha661_16.mat','Est'); - movefile('SwingQSOLb_PN10_alpha661_16.mat','./ResultsHPC/') +% [Est] = UKF_CovarianceAnalysis(Est0,@Dynamics_MPHandMMX_Inertia,... +% @Cov_Dynamics_Good, @New_Observables_model,... +% par.R,YObs_Full,par,units); +% +% save('QSOM_PN10_alpha1_16_Lidar20Min.mat','Est'); +% movefile('QSOM_PN10_alpha1_16_Lidar20Min.mat','./ResultsHPC/') %% QSO-M - cspice_furnsh('./MMX_Product/MMX_BSP_Files_GravLib/MMX_QSO_049_2x2_826891269_828619269.bsp'); - load("./ObservationsHPC/Observations-QSOLa_16_NORange.mat"); +% cspice_furnsh('./MMX_Product/MMX_BSP_Files_GravLib/MMX_QSO_049_2x2_826891269_828619269.bsp'); + load("./ObservationsHPC/Observations-QSOM_16_Lidar5Min.mat"); % Model parameters [par, units] = MMX_DefineNewModelParametersAndUnits; @@ -84,6 +85,7 @@ clc % Covariance analysis parameters [par, units] = New_MMX_CovarianceAnalysisParameters(par, units); + par.interval_lidar = 300; par.sigma = 1e-10/(units.vsf*units.tsf); par.sigmaPh = 0/(units.vsf*units.tsf); @@ -102,19 +104,19 @@ clc Est0.P0 = par.P0; Est0.t0 = data*units.tsf; - par.alpha = 7e-1; + par.alpha = 1; - [Est] = UKF_CovarianceAnalysis(Est0,@Dynamics_MPHandMMX_Inertia,... - @Cov_Dynamics_Good, @New_Observables_model,... - par.R,YObs_Full,par,units); - - save('QSOLa_PN10_alpha71_16_NORange.mat','Est'); - movefile('QSOLa_PN10_alpha71_16_NORange.mat','./ResultsHPC/') +% [Est] = UKF_CovarianceAnalysis(Est0,@Dynamics_MPHandMMX_Inertia,... +% @Cov_Dynamics_Good, @New_Observables_model,... +% par.R,YObs_Full,par,units); +% +% save('QSOM_PN10_alpha1_16_Lidar5Min.mat','Est'); +% movefile('QSOM_PN10_alpha1_16_Lidar5Min.mat','./ResultsHPC/') %% QSO-La - cspice_furnsh('./MMX_Product/MMX_BSP_Files_GravLib/MMX_QSO_027_2x2_826891269_828619269.bsp'); - load("./ObservationsHPC/Observations-QSOLc_16_NOLOS.mat"); +% cspice_furnsh('./MMX_Product/MMX_BSP_Files_GravLib/MMX_QSO_027_2x2_826891269_828619269.bsp'); + load("./ObservationsHPC/Observations-QSOM_16_Lidar10Min.mat"); % Model parameters [par, units] = MMX_DefineNewModelParametersAndUnits; @@ -127,6 +129,7 @@ clc % Covariance analysis parameters [par, units] = New_MMX_CovarianceAnalysisParameters(par, units); + par.interval_lidar = 600; par.sigma = 1e-10/(units.vsf*units.tsf); par.sigmaPh = 0/(units.vsf*units.tsf); @@ -145,19 +148,19 @@ clc Est0.P0 = par.P0; Est0.t0 = data*units.tsf; - par.alpha = 3e-1; + par.alpha = 1; - [Est] = UKF_CovarianceAnalysis(Est0,@Dynamics_MPHandMMX_Inertia,... - @Cov_Dynamics_Good, @New_Observables_model,... - par.R,YObs_Full,par,units); - - save('QSOLc_PN10_alpha31_16_NOLOS.mat','Est'); - movefile('QSOLc_PN10_alpha31_16_NOLOS.mat','./ResultsHPC/') +% [Est] = UKF_CovarianceAnalysis(Est0,@Dynamics_MPHandMMX_Inertia,... +% @Cov_Dynamics_Good, @New_Observables_model,... +% par.R,YObs_Full,par,units); +% +% save('QSOM_PN10_alpha1_16_Lidar10Min.mat','Est'); +% movefile('QSOM_PN10_alpha1_16_Lidar10Min.mat','./ResultsHPC/') %% QSO-Lb - cspice_furnsh('./MMX_Product/MMX_BSP_Files_GravLib/MMX_QSO_031_2x2_826891269_828619269.bsp'); - load("./ObservationsHPC/Observations-QSOLb_16_NOLOS.mat"); +% cspice_furnsh('./MMX_Product/MMX_BSP_Files_GravLib/MMX_QSO_031_2x2_826891269_828619269.bsp'); + load("./ObservationsHPC/Observations-QSOM_16_Lidar1Min.mat"); % Model parameters [par, units] = MMX_DefineNewModelParametersAndUnits; @@ -170,6 +173,7 @@ clc % Covariance analysis parameters [par, units] = New_MMX_CovarianceAnalysisParameters(par, units); + par.interval_lidar = 60; par.sigma = 1e-10/(units.vsf*units.tsf); par.sigmaPh = 0/(units.vsf*units.tsf); @@ -188,19 +192,19 @@ clc Est0.P0 = par.P0; Est0.t0 = data*units.tsf; - par.alpha = 7.7e-1; + par.alpha = 1; [Est] = UKF_CovarianceAnalysis(Est0,@Dynamics_MPHandMMX_Inertia,... @Cov_Dynamics_Good, @New_Observables_model,... par.R,YObs_Full,par,units); - save('QSOLb_PN10_alpha771_16_NOLOS.mat','Est'); - movefile('QSOLb_PN10_alpha771_16_NOLOS.mat','./ResultsHPC/') + save('QSOM_PN10_alpha1_16_Lidar1Min.mat','Est'); + movefile('QSOM_PN10_alpha1_16_Lidar1Min.mat','./ResultsHPC/') %% QSO-Lc - cspice_furnsh('./MMX_Product/MMX_BSP_Files_GravLib/MMX_QSO_049_2x2_826891269_828619269.bsp'); - load("./ObservationsHPC/Observations-QSOLa_16_NOLOS.mat"); +% cspice_furnsh('./MMX_Product/MMX_BSP_Files_GravLib/MMX_QSO_049_2x2_826891269_828619269.bsp'); + load("./ObservationsHPC/Observations-QSOM_16.mat"); % Model parameters [par, units] = MMX_DefineNewModelParametersAndUnits; @@ -213,6 +217,7 @@ clc % Covariance analysis parameters [par, units] = New_MMX_CovarianceAnalysisParameters(par, units); + par.interval_Range_rate = 600; par.sigma = 1e-10/(units.vsf*units.tsf); par.sigmaPh = 0/(units.vsf*units.tsf); @@ -231,20 +236,20 @@ clc Est0.P0 = par.P0; Est0.t0 = data*units.tsf; - par.alpha = 6e-1; + par.alpha = 1; -% [Est] = UKF_CovarianceAnalysis(Est0,@Dynamics_MPHandMMX_Inertia,... -% @Cov_Dynamics_Good, @New_Observables_model,... -% par.R,YObs_Full,par,units); -% -% save('QSOLa_PN10_alpha61_16_NOLOS.mat','Est'); -% movefile('QSOLa_PN10_alpha61_16_NOLOS.mat','./ResultsHPC/') + [Est] = UKF_CovarianceAnalysis(Est0,@Dynamics_MPHandMMX_Inertia,... + @Cov_Dynamics_Good, @New_Observables_model,... + par.R,YObs_Full,par,units); + + save('QSOM_PN10_alpha1_16_Rate10Min.mat','Est'); + movefile('QSOM_PN10_alpha1_16_Rate10Min.mat','./ResultsHPC/') %% QSO-H - cspice_furnsh('./MMX_Product/MMX_BSP_Files_GravLib/MMX_QSO_198_2x2_826891269_828619269.bsp'); - load("./ObservationsHPC/Observations-QSOH_16_DSN.mat"); +% cspice_furnsh('./MMX_Product/MMX_BSP_Files_GravLib/MMX_QSO_198_2x2_826891269_828619269.bsp'); + load("./ObservationsHPC/Observations-QSOM_16_Rate5Min.mat"); % Model parameters [par, units] = MMX_DefineNewModelParametersAndUnits; @@ -257,6 +262,7 @@ clc % Covariance analysis parameters [par, units] = New_MMX_CovarianceAnalysisParameters(par, units); + par.interval_rate = 300; par.sigma = 1e-10/(units.vsf*units.tsf); par.sigmaPh = 0/(units.vsf*units.tsf); @@ -275,19 +281,19 @@ clc Est0.P0 = par.P0; Est0.t0 = data*units.tsf; - par.alpha = 1e-1; + par.alpha = 1; -% [Est] = UKF_CovarianceAnalysis(Est0,@Dynamics_MPHandMMX_Inertia,... -% @Cov_Dynamics_Good, @New_Observables_model,... -% par.R,YObs_Full,par,units); -% -% save('QSOH_PN10_alpha1_16_DSN.mat','Est'); -% movefile('QSOH_PN10_alpha1_16_DSN.mat','./ResultsHPC/') + [Est] = UKF_CovarianceAnalysis(Est0,@Dynamics_MPHandMMX_Inertia,... + @Cov_Dynamics_Good, @New_Observables_model,... + par.R,YObs_Full,par,units); + + save('QSOM_PN10_alpha1_16_Rate5Min.mat','Est'); + movefile('QSOM_PN10_alpha1_16_Rate5Min.mat','./ResultsHPC/') %% QSO-M - cspice_furnsh('./MMX_Product/MMX_BSP_Files_GravLib/MMX_QSO_094_2x2_826891269_828619269.bsp'); - load("./ObservationsHPC/Observations-QSOM_16_DSN.mat"); +% cspice_furnsh('./MMX_Product/MMX_BSP_Files_GravLib/MMX_QSO_094_2x2_826891269_828619269.bsp'); + load("./ObservationsHPC/Observations-QSOM_16_Rate1Min.mat"); % Model parameters [par, units] = MMX_DefineNewModelParametersAndUnits; @@ -300,6 +306,8 @@ clc % Covariance analysis parameters [par, units] = New_MMX_CovarianceAnalysisParameters(par, units); + par.interval_Range_rate = 60; + par.interval_lidar = 60; par.sigma = 1e-10/(units.vsf*units.tsf); par.sigmaPh = 0/(units.vsf*units.tsf); @@ -318,19 +326,21 @@ clc Est0.P0 = par.P0; Est0.t0 = data*units.tsf; - par.alpha = 2e-1; -% -% [Est] = UKF_CovarianceAnalysis(Est0,@Dynamics_MPHandMMX_Inertia,... -% @Cov_Dynamics_Good, @New_Observables_model,... -% par.R,YObs_Full,par,units); -% -% save('QSOM_PN10_alpha1_16_DSN.mat','Est'); -% movefile('QSOM_PN10_alpha1_16_DSN.mat','./ResultsHPC/') + par.alpha = 1; + + [Est] = UKF_CovarianceAnalysis(Est0,@Dynamics_MPHandMMX_Inertia,... + @Cov_Dynamics_Good, @New_Observables_model,... + par.R,YObs_Full,par,units); + + save('QSOM_PN10_alpha1_16_Rate1Min.mat','Est'); + movefile('QSOM_PN10_alpha1_16_Rate1Min.mat','./ResultsHPC/') + + %% QSO-La - cspice_furnsh('./MMX_Product/MMX_BSP_Files_GravLib/MMX_QSO_049_2x2_826891269_828619269.bsp'); - load("./ObservationsHPC/Observations-QSOLa_16_DSN.mat"); +% cspice_furnsh('./MMX_Product/MMX_BSP_Files_GravLib/MMX_QSO_049_2x2_826891269_828619269.bsp'); + load("./ObservationsHPC/Observations-QSOM_16_Camera20Min.mat"); % Model parameters [par, units] = MMX_DefineNewModelParametersAndUnits; @@ -343,6 +353,7 @@ clc % Covariance analysis parameters [par, units] = New_MMX_CovarianceAnalysisParameters(par, units); + pars.interval_camera = 1200; % s par.sigma = 1e-10/(units.vsf*units.tsf); par.sigmaPh = 0/(units.vsf*units.tsf); @@ -361,19 +372,21 @@ clc Est0.P0 = par.P0; Est0.t0 = data*units.tsf; - par.alpha = 1.5e-1; + par.alpha = .8; + + [Est] = UKF_CovarianceAnalysis(Est0,@Dynamics_MPHandMMX_Inertia,... + @Cov_Dynamics_Good, @New_Observables_model,... + par.R,YObs_Full,par,units); + + save('QSOM_PN10_alpha81_16_Camera20Min.mat','Est'); + movefile('QSOM_PN10_alpha81_16_Camera20Min.mat','./ResultsHPC/') + -% [Est] = UKF_CovarianceAnalysis(Est0,@Dynamics_MPHandMMX_Inertia,... -% @Cov_Dynamics_Good, @New_Observables_model,... -% par.R,YObs_Full,par,units); -% -% save('QSOLa_PN10_alpha1_16_DSN.mat','Est'); -% movefile('QSOLa_PN10_alpha1_16_DSN.mat','./ResultsHPC/') %% QSO-Lb - cspice_furnsh('./MMX_Product/MMX_BSP_Files_GravLib/MMX_QSO_031_2x2_826891269_828619269.bsp'); - load("./ObservationsHPC/Observations-QSOLb_16.mat"); +% cspice_furnsh('./MMX_Product/MMX_BSP_Files_GravLib/MMX_QSO_031_2x2_826891269_828619269.bsp'); + load("./ObservationsHPC/Observations-QSOM_16_Camera20Min.mat"); % Model parameters [par, units] = MMX_DefineNewModelParametersAndUnits; @@ -386,6 +399,7 @@ clc % Covariance analysis parameters [par, units] = New_MMX_CovarianceAnalysisParameters(par, units); + pars.interval_camera = 600; % s par.sigma = 1e-10/(units.vsf*units.tsf); par.sigmaPh = 0/(units.vsf*units.tsf); @@ -404,19 +418,20 @@ clc Est0.P0 = par.P0; Est0.t0 = data*units.tsf; - par.alpha = 7.8e-1; + par.alpha = .8; + + [Est] = UKF_CovarianceAnalysis(Est0,@Dynamics_MPHandMMX_Inertia,... + @Cov_Dynamics_Good, @New_Observables_model,... + par.R,YObs_Full,par,units); + + save('QSOM_PN10_alpha81_16_Camera10Min.mat','Est'); + movefile('QSOM_PN10_alpha81_16_Camera10Min.mat','./ResultsHPC/') -% [Est] = UKF_CovarianceAnalysis(Est0,@Dynamics_MPHandMMX_Inertia,... -% @Cov_Dynamics_Good, @New_Observables_model,... -% par.R,YObs_Full,par,units); -% -% save('QSOLb_PN10_alpha781_16.mat','Est'); -% movefile('QSOLb_PN10_alpha781_16.mat','./ResultsHPC/') %% QSO-Lc - cspice_furnsh('./MMX_Product/MMX_BSP_Files_GravLib/MMX_QSO_027_2x2_826891269_828619269.bsp'); - load("./ObservationsHPC/Observations-QSOLc_16_DSN.mat"); +% cspice_furnsh('./MMX_Product/MMX_BSP_Files_GravLib/MMX_QSO_027_2x2_826891269_828619269.bsp'); + load("./ObservationsHPC/Observations-QSOM_16_Camera2Min.mat"); % Model parameters [par, units] = MMX_DefineNewModelParametersAndUnits; @@ -429,6 +444,7 @@ clc % Covariance analysis parameters [par, units] = New_MMX_CovarianceAnalysisParameters(par, units); + pars.interval_camera = 100; % s par.sigma = 1e-10/(units.vsf*units.tsf); par.sigmaPh = 0/(units.vsf*units.tsf); @@ -447,11 +463,11 @@ clc Est0.P0 = par.P0; Est0.t0 = data*units.tsf; - par.alpha = 2e-1; + par.alpha = .8; -% [Est] = UKF_CovarianceAnalysis(Est0,@Dynamics_MPHandMMX_Inertia,... -% @Cov_Dynamics_Good, @New_Observables_model,... -% par.R,YObs_Full,par,units); -% -% save('QSOLc_PN10_alpha1_16_DSN.mat','Est'); -% movefile('QSOLc_PN10_alpha1_16_DSN.mat','./ResultsHPC/') + [Est] = UKF_CovarianceAnalysis(Est0,@Dynamics_MPHandMMX_Inertia,... + @Cov_Dynamics_Good, @New_Observables_model,... + par.R,YObs_Full,par,units); + + save('QSOM_PN10_alpha81_16_Camera2Min.mat','Est'); + movefile('QSOM_PN10_alpha81_16_Camera2Min.mat','./ResultsHPC/') diff --git a/MMX_Fcn_CovarianceAnalyses/NewCov_ResultsPlot.m b/MMX_Fcn_CovarianceAnalyses/NewCov_ResultsPlot.m index 7881aaa3a98e2ad710cde039cf4fd0230780c032..7ccaa4050b8b0993192fd7604a02bb3f6edff680 100644 --- a/MMX_Fcn_CovarianceAnalyses/NewCov_ResultsPlot.m +++ b/MMX_Fcn_CovarianceAnalyses/NewCov_ResultsPlot.m @@ -250,11 +250,11 @@ function NewCov_ResultsPlot(Est, YObs_Full, pars, units) % Biases -% place0 = size(Est.X,1)-pars.nBias; + place0 = size(Est.X,1)-pars.nBias; % corr_label(end+1) = {'$\rho$'}; % corr_label(end+1) = {'$\dot{\rho}$'}; % corr_label(end+1) = {'$LIDAR_b$'}; -% + % figure() % subplot(1,pars.nBias,1) % semilogy(t_obs(1:end-1)/3600,3.*squeeze(real(sqrt(Est.P_t(place0+1,place0+1,1:end-1)))),'LineWidth',1,... @@ -277,8 +277,8 @@ function NewCov_ResultsPlot(Est, YObs_Full, pars, units) % hold on; % xlabel('time $[hour]$') % ylabel('$[km]$') -% -% + + % % % Harmonics Coefficients % diff --git a/MMX_Fcn_CovarianceAnalyses/New_MMX_CovarianceAnalysisParameters.m b/MMX_Fcn_CovarianceAnalyses/New_MMX_CovarianceAnalysisParameters.m index f1c780d64c1c711da0fe3dcbdb874166f61909c0..33884d83e4475722ae1acca2984bbcab8623d857 100644 --- a/MMX_Fcn_CovarianceAnalyses/New_MMX_CovarianceAnalysisParameters.m +++ b/MMX_Fcn_CovarianceAnalyses/New_MMX_CovarianceAnalysisParameters.m @@ -82,18 +82,23 @@ pars.bias = bias; % Unit vector and parameters number useful for later use pars.nCoeff = size(pars.I2,1); pars.nBias = size(bias,1); +% pars.nBias = 0; units.Dim = [units.sfVec; units.lsf; units.vsf; 1; units.tsf;... 1; units.tsf; 1; units.tsf; ones(pars.nCoeff,1); units.lsf; units.vsf;... units.lsf]; +% units.Dim = [units.sfVec; units.lsf; units.vsf; 1; units.tsf;... +% 1; units.tsf; 1; units.tsf; ones(pars.nCoeff,1)]; units.CovDim = (repmat(units.Dim, 1, length(units.Dim)).*... repmat(units.Dim', length(units.Dim), 1)); % A-priori covariance k = 1e-1; -pars.P0 = diag(([1*ones(3,1); 1e-3*ones(3,1);... - 1e-1; 1e-4; 1e-6; 1e-5; 3e-1; 1e-3; 3e-1; 1e-3; 1e-2*ones(pars.nCoeff,1)/k;... - 1; 1e-3; 1]*k./units.Dim).^2); +pars.P0 = diag(([.3*ones(3,1); .3e-3*ones(3,1);... + 1e-1; 1e-4; 1e-5; 1e-7; 3e-1; 1e-3; 3e-1; 1e-3; 1e-2*ones(pars.nCoeff,1);... + 1; 1e-3; 1]./units.Dim).^2); +% pars.P0 = diag(([.3*ones(3,1); .3e-3*ones(3,1);... +% 1e-1; 1e-4; 1e-6; 1e-8; 3e-1; 1e-3; 3e-1; 1e-3; 1e-2*ones(pars.nCoeff,1)]./units.Dim).^2); % PN coefficients diff --git a/MMX_Fcn_CovarianceAnalyses/New_Observables_model.m b/MMX_Fcn_CovarianceAnalyses/New_Observables_model.m index c134fcf4a423fada16940c28fa432a58b2113079..8b35543c7a900c901f727edc542fabb986f70262 100644 --- a/MMX_Fcn_CovarianceAnalyses/New_Observables_model.m +++ b/MMX_Fcn_CovarianceAnalyses/New_Observables_model.m @@ -58,6 +58,8 @@ function G = New_Observables_model(et,Obs,X,par,units) r_Phobos = par.perifocal2MARSIAU*X(d+1)*[cos(X(d+3)); sin(X(d+3)); 0]; switch size(X,1) + case 17 + bias = par.bias; case 18 bias = X(16:end); case 20 diff --git a/MMX_Fcn_CovarianceAnalyses/UKF_CovarianceAnalysis.m b/MMX_Fcn_CovarianceAnalyses/UKF_CovarianceAnalysis.m index 00f5e4f785214e300657c1bc70f082d2a7972db4..5865ebb38ef7d2b7da86e5292e1b4fbe6cbd528a 100644 --- a/MMX_Fcn_CovarianceAnalyses/UKF_CovarianceAnalysis.m +++ b/MMX_Fcn_CovarianceAnalyses/UKF_CovarianceAnalysis.m @@ -91,6 +91,8 @@ function [Est] = UKF_CovarianceAnalysis(Est0,f,fsigma,G,O,Measures,par,units) switch n case 14 St0 = [Xold; par.IPhx_bar; par.IPhy_bar; par.IPhz_bar]; + case 17 + St0 = Xold(1:17); case 18 St0 = par.X0_reference; case 20 @@ -137,6 +139,8 @@ function [Est] = UKF_CovarianceAnalysis(Est0,f,fsigma,G,O,Measures,par,units) switch n case 14 Xmean_bar = X(idx,1:size(Est0.X,1))'; + case 17 + Xmean_bar = X(idx,:)'; case 18 Xmean_bar = [X(idx,1:10)'; X(idx,13:end)'; 0; 0; 0]; par.PhiM = X(idx,11); @@ -148,6 +152,8 @@ function [Est] = UKF_CovarianceAnalysis(Est0,f,fsigma,G,O,Measures,par,units) switch n case 14 Xmean_bar = X(idx,1:size(Est0.X,1))'; + case 17 + Xmean_bar = X(idx,:)'; case 18 Xmean_bar = [X(idx,1:10)'; X(idx,13:end)'; 0; 0; 0]; case 20 @@ -213,6 +219,10 @@ function [Est] = UKF_CovarianceAnalysis(Est0,f,fsigma,G,O,Measures,par,units) Gamma = [eye(q)*deltaT^2/2; eye(q)*deltaT; zeros(n-6,q)]; Q = Gamma*Q*Gamma' + [zeros(6,n); zeros(n-6,6), par.sigmaPh^2*... diag([deltaT^2/2, deltaT, deltaT^2/2, deltaT, deltaT^2/2, deltaT, deltaT^2/2, deltaT])]; + case 17 + Gamma = [eye(q)*deltaT^2/2; eye(q)*deltaT; zeros(n-6,q)]; + Q = Gamma*Q*Gamma' + [zeros(6,n); zeros(n-6,6), par.sigmaPh^2*... + diag([deltaT^2/2, deltaT, deltaT^2/2, deltaT, deltaT^2/2, deltaT, deltaT^2/2, deltaT, zeros(1,n-14)])]; case 18 Gamma = [eye(q)*deltaT^2/2; eye(q)*deltaT; zeros(n-6,q)]; Q = Gamma*Q*Gamma' + [zeros(6,n); zeros(n-6,6), par.sigmaPh^2*... @@ -270,8 +280,8 @@ function [Est] = UKF_CovarianceAnalysis(Est0,f,fsigma,G,O,Measures,par,units) Pold = P; - err(IDX,i) = (Y_obs - G(et,obs(:,i),Xstar,par,units)').*units.Observations(IDX); - prefit(IDX,i) = (Y_obs - G(et,obs(:,i),Xmean_bar,par,units)').*units.Observations(IDX); + err(IDX,i) = (Y_obs - G(et,obs(:,idx),Xstar,par,units)').*units.Observations(IDX); + prefit(IDX,i) = (Y_obs - G(et,obs(:,idx),Xmean_bar,par,units)').*units.Observations(IDX); y_t(IDX,i) = y.*units.Observations(IDX); diff --git a/Model_Check/.DS_Store b/Model_Check/.DS_Store new file mode 100644 index 0000000000000000000000000000000000000000..5008ddfcf53c02e82d7eee2e57c38e5672ef89f6 Binary files /dev/null and b/Model_Check/.DS_Store differ diff --git a/Model_Check/TestIntegratorNew.m b/Model_Check/TestIntegratorNew.m index 9c44ae6c37bee4d7ddb70a9747dc43055f8e9a3a..593553e791be724d0425e90213595f774eb3ba16 100644 --- a/Model_Check/TestIntegratorNew.m +++ b/Model_Check/TestIntegratorNew.m @@ -193,14 +193,14 @@ clc grid on; hold on; plot(t/3600,r_Ph_t); - ylabel('radius [km]') - xlabel('[hours]') - legend('Integrated','Spice') + ylabel('$r_{Ph}$ [km]','Fontsize',26) + xlabel('[hours]','Fontsize',26) + legend('Integrated','Spice','Fontsize',22) subplot(2,4,5) plot(t/3600,X(:,1)*units.lsf-r_Ph_t(:)); grid on; - ylabel('err radius [km]') - xlabel('[hours]') + ylabel('$\Delta r_{Ph}$ [km]','Fontsize',26) + xlabel('[hours]','Fontsize',26) % R_dot Spice and integration subplot(2,4,2) @@ -208,14 +208,14 @@ clc grid on; hold on; plot(t/3600,rdot_Ph_t); - ylabel('Rdot [km/s]') - xlabel('[hours]') - legend('Integrated','Spice') + ylabel('$\dot{r}_{Ph}$ [km/s]','Fontsize',26) + xlabel('[hours]','Fontsize',26) + legend('Integrated','Spice','Fontsize',22) subplot(2,4,6) plot(t/3600,X(:,2)*units.vsf-rdot_Ph_t(:)); grid on; - ylabel('err Rdot [km/s]') - xlabel('[hours]') + ylabel('$\Delta \dot{r}_{Ph}$ [km/s]','Fontsize',26) + xlabel('[hours]','Fontsize',26) % Theta Spice and integration subplot(2,4,3) @@ -223,14 +223,14 @@ clc grid on; hold on; plot(t/3600,mod(theta_t,2*pi)); - ylabel('Theta [rad]') - xlabel('[hours]') - legend('Integrated','Spice') + ylabel('$\theta$ [rad]','Fontsize',26) + xlabel('[hours]','Fontsize',26) + legend('Integrated','Spice','Fontsize',22) subplot(2,4,7) plot(t/3600,angdiff(X(:,3)',theta_t)); grid on; - ylabel('err Theta [rad]') - xlabel('[hours]') + ylabel('$\Delta \theta$ [rad]','Fontsize',26) + xlabel('[hours]','Fontsize',26) % Theta_dot Spice and integration subplot(2,4,4) @@ -238,14 +238,14 @@ clc grid on; hold on; plot(t/3600,thetadot_t); - ylabel('Thetadot [rad/s]') - xlabel('[hours]') - legend('Integrated','Spice') + ylabel('$\dot{\theta}$ [rad/s]','Fontsize',26) + xlabel('[hours]','Fontsize',26) + legend('Integrated','Spice','Fontsize',22) subplot(2,4,8) plot(t/3600,X(:,4)*units.tsf-thetadot_t(:)); grid on; - ylabel('err Thetadot [rad/s]') - xlabel('[hours]') + ylabel('$\Delta \dot{\theta}$ [rad/s]','Fontsize',26) + xlabel('[hours]','Fontsize',26) % Phi1 Spice and integration figure() @@ -254,14 +254,14 @@ clc grid on; hold on; plot(t/3600,mod(Phi1_t,2*pi)); - ylabel('Phi1 [rad]') - xlabel('[hours]') - legend('Integrated','Spice') + ylabel('$\Phi_{M}$ [rad]','Fontsize',26) + xlabel('[hours]','Fontsize',26) + legend('Integrated','Spice','Fontsize',22) subplot(2,4,5) plot(t/3600,angdiff(X(:,5)',Phi1_t)); grid on; - ylabel('err Phi1 [rad]') - xlabel('[hours]') + ylabel('$\Delta \Phi_{M}$ [rad]','Fontsize',26) + xlabel('[hours]','Fontsize',26) % Phi1_dot Spice and integration subplot(2,4,2) @@ -269,14 +269,14 @@ clc grid on; hold on; plot(t/3600,Phi1dot_t); - ylabel('Phi1dot [rad/s]') - xlabel('[hours]') - legend('Integrated','Spice') + ylabel('$\dot{\Phi}_{M}$ [rad/s]','Fontsize',26) + xlabel('[hours]','Fontsize',26) + legend('Integrated','Spice','Fontsize',22) subplot(2,4,6) plot(t/3600,X(:,6)*units.tsf-Phi1dot_t(:)); grid on; - ylabel('err Phi1dot [rad/s]') - xlabel('[hours]') + ylabel('$\Delta \dot{\Phi}_{M}$ [rad/s]','Fontsize',26) + xlabel('[hours]','Fontsize',26) % Phi2 Spice and integration subplot(2,4,3) @@ -284,14 +284,14 @@ clc grid on; hold on; plot(t/3600,mod(Phi2_t,2*pi)); - ylabel('Phi2 [rad]') - xlabel('[hours]') - legend('Integrated','Spice') + ylabel('$\dot{\Phi}_{Ph}$ [rad]','Fontsize',26) + xlabel('[hours]','Fontsize',26) + legend('Integrated','Spice','Fontsize',22) subplot(2,4,7) plot(t/3600,angdiff(X(:,7)',Phi2_t)); grid on; - ylabel('err Phi2 [rad]') - xlabel('[hours]') + ylabel('$\Delta \dot{\Phi}_{Ph}$ [rad]','Fontsize',26) + xlabel('[hours]','Fontsize',26) % Phi2_dot Spice and integration subplot(2,4,4) @@ -299,14 +299,14 @@ clc grid on; hold on; plot(t/3600,Phi2dot_t); - ylabel('Phi2dot') - xlabel('[hours]') - legend('Integrated','Spice') + ylabel('$\dot{\Phi}_{Ph}$ [rad/s]','Fontsize',26) + xlabel('[hours]','Fontsize',26) + legend('Integrated','Spice','Fontsize',22) subplot(2,4,8) plot(t/3600,X(:,8)*units.tsf-Phi2dot_t(:)); grid on; - ylabel('err Phi2dot') - xlabel('[hours]') + ylabel('$\Delta \dot{\Phi}_{Ph}$ [rad/s]','Fontsize',26) + xlabel('[hours]','Fontsize',26) % Comparison wrt x,y,z in MarsIAU reference frame diff --git a/New_MMXCov_Analysis.m b/New_MMXCov_Analysis.m index 3da930b66f4266d511651ab1375f136a3163e640..fdb21e6655192656fba72da5a3d38f4591e7f04f 100644 --- a/New_MMXCov_Analysis.m +++ b/New_MMXCov_Analysis.m @@ -16,12 +16,12 @@ clc addpath(genpath('../mice/')); addpath(genpath('../dynamical_model/')); addpath(genpath('../generic_kernels/')); - addpath(genpath('../paper_after_conf/MMX_Product/MMX_BSP_Files_GravLib/')); + addpath(genpath('./MMX_Product/MMX_BSP_Files_GravLib/')); MMX_InitializeSPICE cspice_furnsh(which('mar097.bsp')); cspice_furnsh(which('MARPHOFRZ.txt')); - cspice_furnsh(which('MMX_QSO_027_2x2_826891269_828619269.bsp')); -% cspice_furnsh(which('MMX_3DQSO_031_009_2x2_826891269_828619269.bsp')); +% cspice_furnsh(which('MMX_QSO_027_2x2_826891269_828619269.bsp')); + cspice_furnsh(which('MMX_3DQSO_031_009_2x2_826891269_828619269.bsp')); % cspice_furnsh(which('MMX_SwingQSO_031_011_2x2_826891269_831211269.bsp')); cspice_furnsh(which('Phobos_826891269_828619269.bsp')); @@ -29,8 +29,8 @@ clc %% Load observations and path for syntethic pictures -% load("./ObservationsHPC/Observations-3DQSOLb_16.mat"); - load("./Observations.mat"); + load("./ObservationsHPC/Observations-3DQSOLb_16.mat"); +% load("./Observations.mat"); %% Initial conditions for the analysis diff --git a/New_MMXCov_ObsGen.m b/New_MMXCov_ObsGen.m index 26d917ea36ff486ccae74886fa49029df7871d45..04f78970fc6f5f71985b3782e313c3a28c99c195 100644 --- a/New_MMXCov_ObsGen.m +++ b/New_MMXCov_ObsGen.m @@ -18,16 +18,16 @@ clc MMX_InitializeSPICE cspice_furnsh(which('mar097.bsp')); cspice_furnsh(which('MARPHOFRZ.txt')); - cspice_furnsh(which('MMX_QSO_027_2x2_826891269_828619269.bsp')); + cspice_furnsh(which('MMX_QSO_094_2x2_826891269_828619269.bsp')); % cspice_furnsh(which('MMX_3DQSO_094_009_2x2_826891269_828619269.bsp')); -% cspice_furnsh(which('MMX_SwingQSO_094_006_2x2_826891269_828619269.bsp')); +% cspice_furnsh(which('MMX_SwingQSO_031_011_2x2_826891269_828619269.bsp')); cspice_furnsh(which('Phobos_826891269_828619269.bsp')); % Time of the analysis data = '2026-03-16 00:00:00 (UTC)'; data = cspice_str2et(data); day = 86400; - PropTime = 2*day; + PropTime = 7*day; date_end = data+PropTime; % Model parameters @@ -47,5 +47,5 @@ clc % YObs_Full(7,:) = NaN; % YObs_Full(end-2:end,:) = NaN; -% save('./ObservationsHPC/Observations-SwingQSOM_16.mat','YObs_Full'); - save('Observations.mat','YObs_Full'); + save('./ObservationsHPC/Observations-QSOM_16_Lidar5Min.mat','YObs_Full'); +% save('Observations.mat','YObs_Full'); diff --git a/Observations.mat b/Observations.mat index 77280646c5f95fb762c8f31fe2ec2a5f3a7a9991..c7e8cef8a60edcb382b07c9351f2571d04f615eb 100644 Binary files a/Observations.mat and b/Observations.mat differ diff --git a/PlotComparison.m b/PlotComparison.m index db549a1f329395a7b0db90c82a567787c43d0ad4..db89692be1b6725d5d12393c35f6fc56261ad0fa 100644 --- a/PlotComparison.m +++ b/PlotComparison.m @@ -25,7 +25,7 @@ semilogy(t_obs(1:end-1)/3600,SQRT_X5,'LineWidth',1.2) title('Position vector $3\sigma$ RMS','FontSize',30) xlabel('Time $[hour]$','Fontsize',26) ylabel('$[km]$','Fontsize',26) -legend('QSO-Lc','QSO-Lb','QSO-La','QSO-H','QSO-H', 'Interpreter','latex','Fontsize',26) +legend('QSO-H','QSO-M','QSO-La','QSO-Lb','QSO-Lc', 'Interpreter','latex','Fontsize',26) subplot(1,2,2) semilogy(t_obs(1:end-1)/3600,SQRT_V1,'LineWidth',1.2) grid on; @@ -37,7 +37,7 @@ semilogy(t_obs(1:end-1)/3600,SQRT_V5,'LineWidth',1.2) title('Velocity vector $3\sigma$ RMS','FontSize',30) xlabel('Time $[hour]$','Fontsize',26) ylabel('$[km/s]$','Fontsize',26) -legend('QSO-Lc','QSO-Lb','QSO-La','QSO-H','QSO-H', 'Interpreter','latex','Fontsize',26) +legend('QSO-H','QSO-M','QSO-La','QSO-Lb','QSO-Lc', 'Interpreter','latex','Fontsize',26) figure() subplot(1,2,1) @@ -50,7 +50,7 @@ semilogy(t_obs(1:end-1)/3600,3.*real(sqrt(P_t_C_20_4)),'LineWidth',1.2); semilogy(t_obs(1:end-1)/3600,3.*real(sqrt(P_t_C_20_5)),'LineWidth',1.2); xlabel('Time $[hour]$','Fontsize',26) ylabel('$C_{20}$','Fontsize',26) -legend('QSO-Lc','QSO-Lb','QSO-La','QSO-H','QSO-H', 'Interpreter','latex','Fontsize',26) +legend('QSO-H','QSO-M','QSO-La','QSO-Lb','QSO-Lc', 'Interpreter','latex','Fontsize',26) subplot(1,2,2) semilogy(t_obs(1:end-1)/3600,3.*real(sqrt(P_t_C_22_1)),'LineWidth',1.2); grid on; @@ -61,7 +61,7 @@ semilogy(t_obs(1:end-1)/3600,3.*real(sqrt(P_t_C_22_4)),'LineWidth',1.2); semilogy(t_obs(1:end-1)/3600,3.*real(sqrt(P_t_C_22_5)),'LineWidth',1.2); xlabel('Time $[hour]$','Fontsize',26) ylabel('$C_{22}$','Fontsize',26) -legend('QSO-Lc','QSO-Lb','QSO-La','QSO-H','QSO-H', 'Interpreter','latex','Fontsize',26) +legend('QSO-H','QSO-M','QSO-La','QSO-Lb','QSO-Lc', 'Interpreter','latex','Fontsize',26) figure() @@ -75,7 +75,7 @@ semilogy(t_obs(1:end-1)/3600,180/pi*3.*squeeze(real(sqrt(Est5.P_t(13,13,1:end-1) title('Libration angle $3\sigma$ envelopes','Interpreter','latex','Fontsize',30) xlabel('Time [hour]', 'Interpreter','latex','Fontsize',26) ylabel('$\Phi_{Ph}$ [deg]', 'Interpreter','latex','Fontsize',26) -legend('QSO-Lc','QSO-Lb','QSO-La','QSO-H','QSO-H', 'Interpreter','latex','Fontsize',26) +legend('QSO-H','QSO-M','QSO-La','QSO-Lb','QSO-Lc', 'Interpreter','latex','Fontsize',26) % Phobos & Mars states diff --git a/PlotFinali_old.m b/PlotFinali_old.m index eb3f724b5614d54447ba638f9016eef38820a3b8..3fa4c98c29be9b4c48fcbf5f8c058360fbd5a463 100644 --- a/PlotFinali_old.m +++ b/PlotFinali_old.m @@ -3,7 +3,8 @@ close all clc addpath('./Results/'); -addpath('./ResultsHPC/'); +addpath('./ResultsHPC_old/'); +addpath('./ObservationsHPC/'); addpath('../useful_functions/'); addpath(genpath('../mice/')); addpath(genpath('../generic_kernels/')); diff --git a/Plot_DifferentFrequencies.m b/Plot_DifferentFrequencies.m new file mode 100644 index 0000000000000000000000000000000000000000..9bbe15f8acab1ff2816cef1ea8fa68cb246a2f9b --- /dev/null +++ b/Plot_DifferentFrequencies.m @@ -0,0 +1,261 @@ +clear +close all +clc + +addpath('./Results/'); +addpath('./ResultsHPC/'); +addpath('./ObservationsHPC/'); +addpath('../useful_functions/'); +addpath(genpath('../mice/')); +addpath(genpath('../generic_kernels/')); +addpath('MMX_Fcn_CovarianceAnalyses/'); +set(0,'DefaultTextInterpreter','latex'); + +%% Confronto fra diverse Lidar frequencies + +load('QSOM_PN10_alpha1_16_Lidar20Min.mat'); +Est20Min = Est; +load('QSOM_PN10_alpha1_16_Lidar10Min.mat'); +Est10Min = Est; +load('QSOM_PN10_alpha1_16_Lidar5Min.mat'); +Est5Min = Est; +load('QSOM_PN10_alpha1_16.mat'); +Est2Min = Est; +load('QSOM_PN10_alpha1_16_Lidar1Min.mat'); +Est1Min = Est; + + +[SQRT_X99, SQRT_V99, P_t_C_20_99, P_t_C_22_99, mean99, std99, meanV99, stdV99,... + meanLib99, stdLib99, meanC2099, stdC2099, meanC2299, stdC2299] = Envelopes_Plot(Est20Min, 500, 500); +[SQRT_X98, SQRT_V98, P_t_C_20_98, P_t_C_22_98, mean98, std98, meanV98, stdV98,... + meanLib98, stdLib98, meanC2098, stdC2098, meanC2298, stdC2298] = Envelopes_Plot(Est10Min, 500, 500); +[SQRT_X97, SQRT_V97, P_t_C_20_97, P_t_C_22_97, mean97, std97, meanV97, stdV97,... + meanLib97, stdLib97, meanC2097, stdC2097, meanC2297, stdC2297] = Envelopes_Plot(Est5Min, 1000, 1000); +[SQRT_X96, SQRT_V96, P_t_C_20_96, P_t_C_22_96, mean96, std96, meanV96, stdV96,... + meanLib96, stdLib96, meanC2096, stdC2096, meanC2296, stdC2296] = Envelopes_Plot(Est2Min, 4000, 4000); +[SQRT_X95, SQRT_V95, P_t_C_20_95, P_t_C_22_95, mean95, std95, meanV95,... + stdV95, meanLib95, stdLib95, meanC2095, stdC2095, meanC2295, stdC2295] = Envelopes_Plot(Est1Min, 8000, 8000); + +xFreq = (1:5); +yFreq = [mean99; mean98; mean97; mean96; mean95]; +yposFreq = [std99; std98; std97; std96; std95]; +ynegFreq = [std99; std98; std97; std96; std95]; +yVFreq = [meanV99; meanV98; meanV97; meanV96; meanV95]; +yposVFreq = [stdV99; stdV98; stdV97; stdV96; stdV95]; +ynegVFreq = [stdV99; stdV98; stdV97; stdV96; stdV95]; + +labelsFreq = {'$20 Min$','$10 Min$','$5 Min$','$2 Min$','$1 Min$'}; + + +figure() +subplot(1,2,1) +errorbar(xFreq,yFreq,ynegFreq,yposFreq,'-s','LineWidth',1.2,'MarkerSize',5) +set(gca,'YScale','log') +set(gca,'TickLabelInterpreter','latex') +set(gca,'XTick',(1:5)); +set(gca,'XTickLabel',labelsFreq,'Fontsize',26); +set(gca,'XTickLabelRotation',45) +ylabel('$[km]$','Fontsize',26) +xlim([0,6]) +grid on +subplot(1,2,2) +errorbar(xFreq,yVFreq,ynegVFreq,yposVFreq,'-s',... + 'LineWidth',1.2,'MarkerSize',5) +grid on +set(gca,'YScale','log') +set(gca,'TickLabelInterpreter','latex') +set(gca,'XTick',(1:5)); +set(gca,'XTickLabel',labelsFreq,'Fontsize',26); +set(gca,'XTickLabelRotation',45) +ylabel('$[km/s]$','Fontsize',26) +xlim([0,6]) + + +%% Confronto fra diverse Range rate frequencies + +load('QSOM_PN10_alpha1_16_Rate10Min.mat'); +Est10Min_rate = Est; +load('QSOM_PN10_alpha1_16_Rate5Min.mat'); +Est5Min_rate = Est; +load('QSOM_PN10_alpha81_16_Rate2Min.mat'); +Est2Min_rate = Est; +load('QSOM_PN10_alpha1_16_Rate1Min.mat'); +Est1Min_rate = Est; + +t_obs = Est10Min_rate.t(1,:); +coeff = 0.8; +fine = round(coeff*size(t_obs,2)); +fineCss = round(coeff*size(t_obs,2)); + +[SQRT_X98, SQRT_V98, P_t_C_20_98, P_t_C_22_98, mean98, std98, meanV98, stdV98,... + meanLib98, stdLib98, meanC2098, stdC2098, meanC2298, stdC2298] = Envelopes_Plot(Est10Min_rate, fine, fineCss); +[SQRT_X97, SQRT_V97, P_t_C_20_97, P_t_C_22_97, mean97, std97, meanV97, stdV97,... + meanLib97, stdLib97, meanC2097, stdC2097, meanC2297, stdC2297] = Envelopes_Plot(Est5Min_rate, fine, fineCss); +[SQRT_X96, SQRT_V96, P_t_C_20_96, P_t_C_22_96, mean96, std96, meanV96, stdV96,... + meanLib96, stdLib96, meanC2096, stdC2096, meanC2296, stdC2296] = Envelopes_Plot(Est2Min_rate, fine, fineCss); +[SQRT_X95, SQRT_V95, P_t_C_20_95, P_t_C_22_95, mean95, std95, meanV95, stdV95,... + meanLib95, stdLib95, meanC2095, stdC2095, meanC2295, stdC2295] = Envelopes_Plot(Est1Min_rate, 7000, 7000); + +xFreq = (1:4); +yFreq = [mean98; mean97; mean96; mean95]; +yposFreq = [std98; std97; std96; std95]; +ynegFreq = [std98; std97; std96; std95]; +yVFreq = [meanV98; meanV97; meanV96; meanV95]; +yposVFreq = [stdV98; stdV97; stdV96; stdV95]; +ynegVFreq = [stdV98; stdV97; stdV96; stdV95]; + +labelsFreq = {'$10 Min$','$5 Min$','$2 Min$','$1 Min$'}; + + +figure() +subplot(1,2,1) +errorbar(xFreq,yFreq,ynegFreq,yposFreq,'-s','LineWidth',1.2,'MarkerSize',5) +set(gca,'YScale','log') +set(gca,'TickLabelInterpreter','latex') +set(gca,'XTick',(1:4)); +set(gca,'XTickLabel',labelsFreq,'Fontsize',26); +set(gca,'XTickLabelRotation',45) +ylabel('$[km]$','Fontsize',26) +xlim([0,5]) +grid on +subplot(1,2,2) +errorbar(xFreq,yVFreq,ynegVFreq,yposVFreq,'-s',... + 'LineWidth',1.2,'MarkerSize',5) +grid on +set(gca,'YScale','log') +set(gca,'TickLabelInterpreter','latex') +set(gca,'XTick',(1:4)); +set(gca,'XTickLabel',labelsFreq,'Fontsize',26); +set(gca,'XTickLabelRotation',45) +ylabel('$[km/s]$','Fontsize',26) +xlim([0,5]) + + +%% Confronto fra diverse Range frequencies + +load('QSOM_PN10_alpha81_16.mat'); +Est1hour_range = Est; +load('QSOM_PN10_alpha81_16_Range20Min.mat'); +Est20Min_range = Est; +load('QSOM_PN10_alpha81_16_Range10Min.mat'); +Est10Min_range = Est; +load('QSOM_PN10_alpha81_16_Range5Min.mat'); +Est5Min_range = Est; +load('QSOM_PN10_alpha81_16_Range2Min.mat'); +Est2Min_range = Est; + +t_obs = Est1hour_range.t(1,:); +coeff = 0.8; +fine = round(coeff*size(t_obs,2)); +fineCss = round(coeff*size(t_obs,2)); + +[SQRT_X99, SQRT_V99, P_t_C_20_99, P_t_C_22_99, mean99, std99, meanV99, stdV99,... + meanLib99, stdLib99, meanC2099, stdC2099, meanC2299, stdC2299] = Envelopes_Plot(Est1hour_range, fine, fineCss); +[SQRT_X98, SQRT_V98, P_t_C_20_98, P_t_C_22_98, mean98, std98, meanV98, stdV98,... + meanLib98, stdLib98, meanC2098, stdC2098, meanC2298, stdC2298] = Envelopes_Plot(Est20Min_range, fine, fineCss); +[SQRT_X97, SQRT_V97, P_t_C_20_97, P_t_C_22_97, mean97, std97, meanV97, stdV97,... + meanLib97, stdLib97, meanC2097, stdC2097, meanC2297, stdC2297] = Envelopes_Plot(Est10Min_range, fine, fineCss); +[SQRT_X96, SQRT_V96, P_t_C_20_96, P_t_C_22_96, mean96, std96, meanV96, stdV96,... + meanLib96, stdLib96, meanC2096, stdC2096, meanC2296, stdC2296] = Envelopes_Plot(Est5Min_range, fine, fineCss); +[SQRT_X95, SQRT_V95, P_t_C_20_95, P_t_C_22_95, mean95, std95, meanV95,... + stdV95, meanLib95, stdLib95, meanC2095, stdC2095, meanC2295, stdC2295] = Envelopes_Plot(Est2Min_range, fine, fineCss); + +xFreq = (1:5); +yFreq = [mean99; mean98; mean97; mean96; mean95]; +yposFreq = [std99; std98; std97; std96; std95]; +ynegFreq = [std99; std98; std97; std96; std95]; +yVFreq = [meanV99; meanV98; meanV97; meanV96; meanV95]; +yposVFreq = [stdV99; stdV98; stdV97; stdV96; stdV95]; +ynegVFreq = [stdV99; stdV98; stdV97; stdV96; stdV95]; + + PlotComparison(Est1hour_range, Est20Min_range, Est10Min_range, Est5Min_range, Est2Min_range, fine, fineCss); + + +labelsFreq = {'$1 hour$','$20 Min$','$10 Min$','$5 Min$','$2 Min$'}; + + +figure() +subplot(1,2,1) +errorbar(xFreq,yFreq,ynegFreq,yposFreq,'-s','LineWidth',1.2,'MarkerSize',5) +set(gca,'YScale','log') +set(gca,'TickLabelInterpreter','latex') +set(gca,'XTick',(1:5)); +set(gca,'XTickLabel',labelsFreq,'Fontsize',26); +set(gca,'XTickLabelRotation',45) +ylabel('$[km]$','Fontsize',26) +xlim([0,6]) +grid on +subplot(1,2,2) +errorbar(xFreq,yVFreq,ynegVFreq,yposVFreq,'-s',... + 'LineWidth',1.2,'MarkerSize',5) +grid on +set(gca,'YScale','log') +set(gca,'TickLabelInterpreter','latex') +set(gca,'XTick',(1:5)); +set(gca,'XTickLabel',labelsFreq,'Fontsize',26); +set(gca,'XTickLabelRotation',45) +ylabel('$[km/s]$','Fontsize',26) +xlim([0,6]) + + +%% Confronto fra diverse Camera frequencies + +load('QSOM_PN10_alpha81_16.mat'); +Est1hour_camera = Est; +load('QSOM_PN10_alpha81_16_Camera20Min.mat'); +Est20Min_camera = Est; +load('QSOM_PN10_alpha81_16_Camera10Min.mat'); +Est10Min_camera = Est; +load('QSOM_PN10_alpha81_16_Camera2Min.mat'); +Est2Min_camera = Est; + +t_obs = Est1hour_camera.t(1,:); +coeff = 0.8; +fine = round(coeff*size(t_obs,2)); +fineCss = round(coeff*size(t_obs,2)); + +[SQRT_X99, SQRT_V99, P_t_C_20_99, P_t_C_22_99, mean99, std99, meanV99, stdV99,... + meanLib99, stdLib99, meanC2099, stdC2099, meanC2299, stdC2299] = Envelopes_Plot(Est1hour_camera, fine, fineCss); +[SQRT_X98, SQRT_V98, P_t_C_20_98, P_t_C_22_98, mean98, std98, meanV98, stdV98,... + meanLib98, stdLib98, meanC2098, stdC2098, meanC2298, stdC2298] = Envelopes_Plot(Est20Min_camera, fine, fineCss); +[SQRT_X97, SQRT_V97, P_t_C_20_97, P_t_C_22_97, mean97, std97, meanV97, stdV97,... + meanLib97, stdLib97, meanC2097, stdC2097, meanC2297, stdC2297] = Envelopes_Plot(Est10Min_camera, fine, fineCss); +[SQRT_X95, SQRT_V95, P_t_C_20_95, P_t_C_22_95, mean95, std95, meanV95,... + stdV95, meanLib95, stdLib95, meanC2095, stdC2095, meanC2295, stdC2295] = Envelopes_Plot(Est2Min_camera, fine, fineCss); + +xFreq = (1:4); +yFreq = [mean99; mean98; mean97;mean95]; +yposFreq = [std99; std98; std97; std95]; +ynegFreq = [std99; std98; std97; std95]; +yVFreq = [meanV99; meanV98; meanV97; meanV95]; +yposVFreq = [stdV99; stdV98; stdV97; stdV95]; +ynegVFreq = [stdV99; stdV98; stdV97; stdV95]; + +% PlotComparison(Est1hour_camera, Est20Min_camera, Est10Min_camera, Est5Min_camera, Est2Min_camera, fine, fineCss); + + +labelsFreq = {'$1 hour$','$20 Min$','$10 Min$','$2 Min$'}; + + +figure() +subplot(1,2,1) +errorbar(xFreq,yFreq,ynegFreq,yposFreq,'-s','LineWidth',1.2,'MarkerSize',5) +set(gca,'YScale','log') +set(gca,'TickLabelInterpreter','latex') +set(gca,'XTick',(1:4)); +set(gca,'XTickLabel',labelsFreq,'Fontsize',26); +set(gca,'XTickLabelRotation',45) +ylabel('$[km]$','Fontsize',26) +xlim([0,5]) +grid on +subplot(1,2,2) +errorbar(xFreq,yVFreq,ynegVFreq,yposVFreq,'-s',... + 'LineWidth',1.2,'MarkerSize',5) +grid on +set(gca,'YScale','log') +set(gca,'TickLabelInterpreter','latex') +set(gca,'XTick',(1:4)); +set(gca,'XTickLabel',labelsFreq,'Fontsize',26); +set(gca,'XTickLabelRotation',45) +ylabel('$[km/s]$','Fontsize',26) +xlim([0,5]) \ No newline at end of file diff --git a/launchMMX.sh b/launchMMX.sh index ec8abeb2d96220a856e3aeb3d82fce5b3534a3c5..a24c9fcbfae64ae90141202347f7e77fe7cc4a9f 100644 --- a/launchMMX.sh +++ b/launchMMX.sh @@ -14,7 +14,7 @@ #SBATCH --mail-type=ALL #SBATCH --mail-user=ec01259@surrey.ac.uk -cd /users/ec01259/MMX_Edo/paper_after_conf +cd /users/ec01259/MMX_Edo/paper module load matlab