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");
 
    }
 
}