NumericalPropagationWithSRP

NumericalPropagationWithSRP

/**

* Example allowing to numerically propagate an orbit with a SRP model force.
* @author goesterjf
*
*/

public class NumericalPropagationWithSRP {

   public static void main(String[] args) throws PatriusException, IOException, ParseException {
       
       // Patrius Dataset initialization (needed for example to get the UTC time)
       PatriusDataset.addResourcesFromPatriusDataset() ;

       // Recovery of the UTC time scale using a "factory" (not to duplicate such unique object)
       final TimeScale TUC = TimeScalesFactory.getUTC();
       
       // Date of the orbit given in UTC time scale)
       final AbsoluteDate date = new AbsoluteDate("2010-01-01T12:00:00.000", TUC);
       
       // Getting the frame with wich will defined the orbit parameters
       // As for time scale, we will use also a "factory".
       final Frame GCRF = FramesFactory.getGCRF();

       // Initial orbit
       final double sma = 7000.e+3;
       final double exc = 0.;
       final double per = sma*(1.-exc);
       final double apo = sma*(1.+exc);
       final double inc = FastMath.toRadians(98.);
       final double pa = FastMath.toRadians(0.);
       final double raan = FastMath.toRadians(0.);
       final double anm = FastMath.toRadians(0.);
       final double MU = Constants.WGS84_EARTH_MU;
       
       final ApsisRadiusParameters par = new ApsisRadiusParameters(per, apo, inc, pa, raan, anm, PositionAngle.MEAN, MU);
       final Orbit iniOrbit = new ApsisOrbit(par, GCRF, date);
       
       // Mass model using an Assembly
       
       final AssemblyBuilder builder = new AssemblyBuilder();
       
       // Initial mass (mandatory to take into account mass for atmospheric force computation)
       final double dryMass = 100.;
       builder.addMainPart("MAIN");
       builder.addProperty(new MassProperty(dryMass), "MAIN");
       

//SPECIFIC

       final double ka = 1.0;
       final double ks = 0.0;
       final double kd = 0.0;
       builder.addProperty(new RadiativeProperty(ka, ks, kd), "MAIN");
       builder.addProperty(new RadiativeIRProperty(ka, ks, kd), "MAIN");
       final double radius = 10.;
       builder.addProperty(new RadiativeSphereProperty(radius), "MAIN");

//SPECIFIC

       final Assembly assembly = builder.returnAssembly();
       

//SPECIFIC

       // Sun ephemeris
       CelestialBody sun = new MeeusSun();
       
       // Definition of the Earth ellipsoid for later SRP computation
       final Frame ITRF = FramesFactory.getITRF();
       final double AE = Constants.WGS84_EARTH_EQUATORIAL_RADIUS;
       final GeometricBodyShape EARTH = new ExtendedOneAxisEllipsoid(AE, Constants.WGS84_EARTH_FLATTENING, ITRF, "EARTH");
       
       // Direct SRP data
       final double dRef = 1.4959787E11;
       final double pRef = 4.5605E-6;
       DirectRadiativeModel rm = new DirectRadiativeModel(assembly);
       SolarRadiationPressureEllipsoid radPres = new SolarRadiationPressureEllipsoid(dRef, pRef, sun, EARTH, rm);
       
       // Rediffused DRP data
       final int inCorona = 1;
       final int inMeridian = 10;
       final IEmissivityModel inEmissivityModel = new KnockeRiesModel();
       final boolean inAlbedo = true;
       final boolean inIr = true;
       final double coefAlbedo = 1.;
       final double coefIr = 1.;
       RediffusedRadiativeModel rdm = new RediffusedRadiativeModel(inAlbedo, inIr, coefAlbedo, coefIr, assembly);
       RediffusedRadiationPressure reDiff =  new RediffusedRadiationPressure(sun, GCRF, inCorona, inMeridian, inEmissivityModel, rdm);

//SPECIFIC

       final MassProvider mm = new MassModel(assembly);

       // We create a spacecratftstate
       final SpacecraftState iniState = new SpacecraftState(iniOrbit, mm);
       
       // Initialization of the Runge Kutta integrator with a 2 s step
       final double pasRk = 2.;
       final FirstOrderIntegrator integrator = new ClassicalRungeKuttaIntegrator(pasRk);

       // Initialization of the propagator
       final NumericalPropagator propagator = new NumericalPropagator(integrator);
       propagator.resetInitialState(iniState);
       
       // Adding additional state (change name add to set for V3.3)
       propagator.setMassProviderEquation(mm);
       
       // Forcing integration using cartesian equations
       propagator.setOrbitType(OrbitType.CARTESIAN);

//SPECIFIC

       // Adding SRP forces
       propagator.addForceModel(radPres);        
       propagator.addForceModel(reDiff);        

//SPECIFIC

       // Propagating 5 periods
       final double dt = 5.*iniOrbit.getKeplerianPeriod();
       final AbsoluteDate finalDate = date.shiftedBy(dt);
       final SpacecraftState finalState = propagator.propagate(finalDate);
       final Orbit finalOrbit = finalState.getOrbit();
       
       // Printing new date and semi major axis
       System.out.println();
       System.out.println("Initial semi major axis = "+iniOrbit.getA()/1000.+" km");
       System.out.println("New date = "+finalOrbit.getDate().toString(TUC)+" deg");
       System.out.println("Final semi major axis = "+finalOrbit.getA()/1000.+" km");
       // Printing mass
       System.out.println();
       System.out.println("Mass = "+finalState.getMass("MAIN")+" kg");

   }
   

}