<< Jacobian and covariance matrices Cookbook Basic orbit determination >>

celestlab >> - Introduction - > Cookbook > Comparison of orbit propagation models

Comparison of orbit propagation models

Comparison of propagation models

1. Comparison of eckhech and lydlp - case 1/2

The following code compares the 2 analytical models "eckhech" and "lydlp" using the same values for the initial mean elements.

The same physical constants are used by the 2 models (5 zonal terms only).

// -------------------------
// eckhech / lydlp
// Same initial mean elements
// -------------------------
J1Jn = CL_dataGet("j1jn");

// Initial state vector
// Epoch (cjd) + circular mean orbital elements 
cjd0 = 20000;
mean_cir0 = [7.e6; 0; 1.e-3; 98*(%pi/180); 0; 0];

// Final dates
cjd = cjd0 + (0 : 60 : 86400) / 86400;

// Propagation using "eckhech"
[osc1, mean1] = CL_ex_propagate("eckhech", "cir", cjd0, mean_cir0, cjd, res="om", ...
                                j1jn=J1Jn(1:5));
[pos1, vel1] = CL_oe_cir2car(osc1); 
[posm1, velm1] = CL_oe_cir2car(mean1); 

// Propagation using "lydlp"
[osc2, mean2] = CL_ex_propagate("lydlp", "cir", cjd0, mean_cir0, cjd, res="om", ...
                                j1jn=J1Jn(1:5));
[pos2, vel2] = CL_oe_cir2car(osc2); 
[posm2, velm2] = CL_oe_cir2car(mean2); 

scf(); 
plot(cjd-cjd0, CL_norm(pos1-pos2), "b");
plot(cjd-cjd0, CL_norm(posm1-posm2), "r");
CL_g_stdaxes(); 

 
// The main effect is on the mean anomaly. 
// A simple correction on the initial mean orbital elements is then: 

alpha = CL_unMod(osc1(6,:), 2 * %pi); 
dalpha = CL_unMod(osc2(6,:) - osc1(6,:), 2 * %pi); 
n = reglin((cjd-cjd0) * 86400, alpha);
[dn, dalpha0] = reglin((cjd-cjd0) * 86400, dalpha);
dsma = (-2/3) * (dn/n) * mean_cir0(1); 
dmean_cir0 = [dsma; 0; 0; 0; 0; dalpha0]; 

[osc2, mean2] = CL_ex_propagate("lydlp", "cir", cjd0, mean_cir0 - dmean_cir0, ...
                                cjd, res="om", j1jn=J1Jn(1:5));
[pos2, vel2] = CL_oe_cir2car(osc2); 
[posm2, velm2] = CL_oe_cir2car(mean2); 

scf(); 
plot(cjd-cjd0, CL_norm(pos1-pos2), "b");
plot(cjd-cjd0, CL_norm(posm1-posm2), "r");
CL_g_stdaxes();

2. Comparison of eckhech and lydlp - case 2/2

The following code compares the 2 analytical models "eckhech" and "lydlp" using the same values for the initial osculating elements.

The same physical constants are used by the 2 models (5 zonal terms only).

The distance between the positions resulting from the 2 propagation models depends on the initial argument of latitude. It is in fact the consequence of the difference between the osculating elements.

// -------------------------
// eckhech / lydlp
// Same initial osculating elements
// -------------------------
J1Jn = CL_dataGet("j1jn");

// Initial state vector
// Epoch (cjd) + circular mean orbital elements 
cjd0 = 20000;
mean_cir0 = [7.e6; 0; 1.e-3; 98*(%pi/180); 0; 0];

// Final dates
cjd = cjd0 + (0 : 30 : 1*86400)/86400;

f = scf(); 

for pso = (0 : 30 : 360) * %pi/180
  mean_cir0(6) = pso; 
  osc_cir0 = CL_ex_mean2osc("eckhech", "cir", mean_cir0, j1jn=J1Jn(1:5)); 

  // Propagation using "eckhech"
  m = CL_ex_osc2mean("eckhech", "cir", osc_cir0, j1jn=J1Jn(1:5));
  osc_cir1 = CL_ex_propagate("eckhech", "cir", cjd0, m, cjd, res="o", ...
                             j1jn=J1Jn(1:5));
  [pos1, vel1] = CL_oe_cir2car(osc_cir1); 

  // Propagation using "lydlp"
  m = CL_ex_osc2mean("lydlp", "cir", osc_cir0, j1jn=J1Jn(1:5));
  osc_cir2 = CL_ex_propagate("lydlp", "cir", cjd0, m, cjd, res="o", ...
                             j1jn=J1Jn(1:5));
  [pos2, vel2] = CL_oe_cir2car(osc_cir2); 

  scf(f); 
  plot(cjd-cjd0, CL_norm(pos1-pos2));
  CL_g_stdaxes(); 
end

3. Comparison of lydlp and STELA

The following code compares the 2 models "lydlp" and "STELA" using the same values for the initial mean elements.

The model constants used are the same: zonal terms up to degree 5.

The conclusion is that the theories behing lydlp and STELA are very close to each other.

// -------------------------
// lydlp / STELA
// Same mean initial elements
// -------------------------
// Initial date (cjd, time scale: TREF)
cjd0 = 20000;

// Keplerian mean orbital elements (frame: ECI)
mean_cir0 = [7.e6; 0; 1.e-3; 98*(%pi/180); 0; %pi/2];

// Final dates
cjd = cjd0 + (0 : 0.5 : 30);

// The constants used come from STELA
cst = CL_stela_getConst(); 

// Propagated mean elements using lydlp 
mean_cir = CL_ex_propagate("lydlp", "cir", cjd0, mean_cir0, cjd, res="m", ...
                           er=cst.er, mu=cst.mu, j1jn=cst.j1jn(1:5));
[posm, velm] = CL_oe_cir2car(mean_cir, mu=cst.mu); 

// Propagated mean elements using STELA 
params = CL_stela_params("none");
params.central_enabled = %t; 
params.zonal_enabled = %t; 
params.zonal_maxDeg = 5; 
params.integrator_step = 0.5 * 86400; 

// Propagation
mean_cir2 = CL_stela_extrap("cir", cjd0, mean_cir0, cjd, params, "m");
[posm2, velm2] = CL_oe_cir2car(mean_cir2, mu=cst.mu); 

// Plot gap between (mean) positions
scf(); 
plot(cjd-cjd0, CL_norm(posm-posm2));

4. Comparison of lydsec and SGP4 (TLE propagation)

The following code compares the 2 models "lydsec" and "SGP4" using the same values for the initial mean elements.

The model used are not exactly the same because SGP4 (potentially) includes the effect of the third body perturbation.

However the results obtained by the two models are very close to each other. This means that the underlying theories are essentially the same.

// -------------------------
// lydsec / SGP4
// Same mean initial elements
// -------------------------
cjd0 = CL_dat_cal2cjd(2015,1,1); 
mean_kep0 = [7.e6; 1.e-3; 98 * %pi/180; %pi/2; 0; 0];

cjd = cjd0 + (0 : 300 : 30*86400)/86400;

// The constants used come from SGP4
cst = CL_tle_getConst(); 
mu = cst.mu; 
j1jn = [0; cst.j2; cst.j3; cst.j4]; 
er = cst.er; 

// Propagated osculating elements using lydsec 
[osc_kep, mean_kep] = CL_ex_propagate("lydsec", "kep", cjd0, mean_kep0, cjd, "om", ...
                                      mu=mu, er=er, j1jn=j1jn);
[pos, vel] = CL_oe_kep2car(osc_kep, mu=mu); 

// Propagated osculating elements using SGP4 
tle = CL_tle_new(); 
tle = CL_tle_setElem(tle, "kep", cjd0, mean_kep0, frame="native"); 
[pos2, vel2] = CL_tle_genEphem(tle, cjd, frame="native"); 

// Plot gap between (osculating) positions
scf(); 
plot(cjd-cjd0, CL_norm(pos2 - pos));


Report an issue
<< Jacobian and covariance matrices Cookbook Basic orbit determination >>