//============================================================
// 
//             Bureau d'études :  ISAE 321
//             =========================== 
//
//============================================================  
// Script scilab pour :
//  -calculer les paramètres orbitaux pour une orbite héliosynchrosne 
//   phasée et gelée (orbites phasées 1 et 2 jour). Un groupe peut tracer 
//   l orbite un jour l autre 2 jours (orbites Mistigri)
//  -calculer les caracteristiques des orbites
//  -extrapoler les orbites
//  -tracer les orbites
//  -estmimer la consommation d ergol
//
//
//              CNES - DCT/SB/MS
//
//============================================================
// Historique
// V1   Version initiale 03/03/2015
//============================================================
//----------------------------------------------------------
// Chargement des constantes et des fonctions 
//----------------------------------------------------------

CL_init();  // Chargement des constantes de CelestLab (cf help CL_init()

//============================================================
//    Caractéristiques des orbites phasées hélio gelées
//============================================================

// Données d'entrée pour une mission heliosynchrone phasee :
t0 = CL_dat_cal2cjd(2016,7,9); 
// Recherche des orbites héliosynchrones phasées 1 ou 2 jours entre 500 et 800km
alt_min = 500.e3; //Zone de recherche 
alt_max = 800.e3;
Qmin = 1; //jour
Qmax = 2; //jours
ecc = 0.001; //Excentricité gelée
sso = 1; //Sun Synchronous Orbit
phasage = CL_op_searchRepeatOrbits(%CL_eqRad+alt_min,%CL_eqRad+alt_max,Qmin,Qmax,ecc,sso);

// Calcul de l'ascension droite pour que le plan de l'orbite ait la bonne heure locale
hloc = 13; // heure locale désirée
gom = CL_op_locTime(t0, 'mlh', hloc, 'ra'); // hloc en heure

//Calcul de l'excentricité gelée (optionnel)
dga = phasage(:,1)';
inc = phasage(:,3)';
[ecc, pom] = CL_op_frozenOrbit(dga, inc);



//anomalie moyenne
anm = [-%pi/2, -%pi/2];


//============================================================
//  Calcul des caractéristiques de l'orbite keplerienne
//============================================================

//Calcul de la période képlerienne
per = CL_kp_params('per',dga);

//Transformation anomalie moyenne vers anomalie vraie 
v = CL_kp_M2v(ecc,anm);

//Caractéristiques de l'orbite képlerienne
res = CL_kp_characteristics(dga, ecc,v);

//============================================================
//  Calcul des caractéristiques de l'orbite kepler + J2
//============================================================

//Intertrace orbitale (J2) et periode nodale (J2) 
[lgap, nodper] = CL_op_paramsJ2(['lgap', 'nodper'],dga, ecc, inc);

//============================================================
//         Extrapolation des orbites
//============================================================

step = 30 / 86400; // pas de propagation 
t = t0:step:t0+2;
//repère inertiel
kep_t1 = CL_ex_secularJ2(t0, [dga(1), ecc(1), inc(1), pom(1), gom(1), anm(1)]', t); // propagation d orbite

//calcul des coordonnées en cartesien.
pos_eci1 = CL_oe_kep2car(kep_t1);      // position cartesienne en repere inertiel
pos_ecf1 = CL_fr_convert("ECI", "ECF", t, pos_eci1); // position en repere terrestre 


kep_t2 = CL_ex_secularJ2(t0, [dga(2), ecc(2), inc(2), pom(2), gom(1), anm(2)]', t); // orbit propagation

pos_eci2 = CL_oe_kep2car(kep_t2);      // position cartesienne en repere inertie
pos_ecf2 = CL_fr_convert("ECI", "ECF", t, pos_eci2); // position en repere terrestre
						     
//============================================================
//         Tracé des orbites
//============================================================						     
						     						     
scf();
CL_plot_ephem(pos_ecf1, color_id = 3);
CL_plot_earthMap();
a = gca();
a.title.text = "Orbite phasée 2 jours";

scf();
CL_plot_ephem(pos_ecf2, color_id = 5);
CL_plot_earthMap();
a = gca();
a.title.text = "Orbite phasée 1 jour";

//============================================================
//            Ecriture des résultats
//============================================================

entete = [  "Eléments orbitaux [km, deg]"
            "dga, exx, inc, argp, adna anom"  ];
params =[];
//Convertion de format pour écriture 
params(1,:) = dga/1e3; //km
params(2,:) = ecc; //km  
params(3,:) = inc * %CL_rad2deg; // deg 
params(4,:) = pom * %CL_rad2deg; // deg  
params(5,:) = gom * %CL_rad2deg; // deg  
params(6,:) = anm * %CL_rad2deg; // deg  
str = [ msprintf("\n%6.6f %6.6f\n", params) ];
  
str = [entete; str];
nomfic = "D:\DATA_IDM_CIC\orbites.txt";  
mputl(str, nomfic);
//disp(str)

//============================================================
//           Budget d'ergol
//============================================================
//Correction des erreurs d injection
Ddga = [10e3, 10e3];
[deltav_mip, dv_mip, anv] = CL_man_dvSma(dga, ecc ,dga + Ddga);

//Cout du maintien à poste (demi grand axe)
alt = dga - %CL_eqRad;
vel = [sqrt(%CL_mu./ dga); zeros(alt); zeros(alt)];
rho = CL_mod_atmUS76Density(alt);
surface = 0.95; //m2
mass = 200; //kg
coefd = 2.7 * surface ./ mass;
acc = CL_fo_dragAcc(vel, rho, [coefd, coefd]);
//Incrément de vitesse = accélération * temps (3 ans)
deltav_map = CL_norm(acc) * 86400 * 365 * 3;

//Manoeuvre de désorbitation
//Orbite rentrant en 25 ans (calcul STELA)
//6991.82750145 0.00120477 98.2693 90.0005859052 122.152359 270
[deltav_desorb, dv_desorb, anv] = CL_man_dvSma(dga(1), 0, 6991.82750145*1000);
deltav_desorb = [deltav_desorb, 0];

//Delta V total (m/s)
deltav_tot = deltav_mip + deltav_map + deltav_desorb;

//Calcul de la consommation (kg)
conso_tot = CL_man_consumption('dm', deltav_tot, mass, 200);

//============================================================
//             Couverture et ouverture senseur
//============================================================
dlon = CL_op_paramsJ2('lgap',dga,ecc,inc); //Ecart entre 2 traces consécutives
dlon(1) = dlon(1)/2; //Orbite phasée 2 jours
//Angle au centre pour couvrir dlon/2
cen = CL_op_equatorialSwath("dlon2cen", dga, inc, dlon/2);
//ouverture instrument
sat = CL_gm_visiParams(dga, ones(dga).*%CL_eqRad,'cen', cen,'sat');
//Vérifier avec la démo
CL_rad2deg(sat)