List of usage examples for org.apache.commons.math3.ode.nonstiff DormandPrince853Integrator DormandPrince853Integrator
public DormandPrince853Integrator(final double minStep, final double maxStep, final double[] vecAbsoluteTolerance, final double[] vecRelativeTolerance)
From source file:Ephemeris.mySlavePropagator.java
public static void propagate(double a, double e, double i, double omega, double raan, double lM, double stepS, double duration, boolean solarBL, boolean aeroBL, JTextArea TF) { try {//from w w w.ja va 2s . c om // configure Orekit Autoconfiguration.configureOrekit(); // Inertial frame Frame inertialFrame = FramesFactory.getEME2000(); // Initial date in UTC time scale TimeScale utc = TimeScalesFactory.getUTC(); AbsoluteDate initialDate = new AbsoluteDate(2016, 01, 01, 12, 00, 00.000, utc); // gravitation coefficient double mu = 3.986004415e+14; // Orbit construction as Keplerian Orbit initialOrbit = new KeplerianOrbit(a, e, i, omega, raan, lM, PositionAngle.MEAN, inertialFrame, initialDate, mu); // Initial state definition SpacecraftState initialState = new SpacecraftState(initialOrbit); // Adaptive step integrator with a minimum step of 0.001 and a maximum step of 1000 final double minStep = 0.001; final double maxstep = 1000.0; final double positionTolerance = 0.1; final double initStep = stepS; final OrbitType propagationType = OrbitType.KEPLERIAN; final double[][] tolerances = NumericalPropagator.tolerances(positionTolerance, initialOrbit, propagationType); AdaptiveStepsizeIntegrator integrator = new DormandPrince853Integrator(minStep, maxstep, tolerances[0], tolerances[1]); integrator.setInitialStepSize(initStep); NumericalPropagator propagator = new NumericalPropagator(integrator); propagator.setOrbitType(propagationType); // Force Model (reduced to perturbing gravity field) final NormalizedSphericalHarmonicsProvider provider = GravityFieldFactory.getNormalizedProvider(10, 10); ForceModel holmesFeatherstone = new HolmesFeatherstoneAttractionModel( FramesFactory.getITRF(IERSConventions.IERS_2010, true), provider); // Aero drag final Frame earthFrame = CelestialBodyFactory.getEarth().getBodyOrientedFrame(); final double ae = provider.getAe(); final SphericalSpacecraft ssc = new SphericalSpacecraft(1, 0.47, 0., 1.2); final OneAxisEllipsoid earth = new OneAxisEllipsoid(ae, Constants.WGS84_EARTH_FLATTENING, earthFrame); Atmosphere atmosphere = new HarrisPriester(CelestialBodyFactory.getSun(), earth); ForceModel aerodrag = new DragForce(atmosphere, ssc); //Solar radiation pressure ForceModel solarPressure = new SolarRadiationPressure(CelestialBodyFactory.getSun(), ae, ssc); // Add force model to the propagator propagator.addForceModel(holmesFeatherstone); if (aeroBL) propagator.addForceModel(aerodrag); if (solarBL) propagator.addForceModel(solarPressure); // Set up initial state in the propagator propagator.setInitialState(initialState); // Set the propagator to slave mode (could be omitted as it is the default mode) propagator.setSlaveMode(); // Overall duration in seconds for extrapolation //double duration = 630.; // Stop date final AbsoluteDate finalDate = initialDate.shiftedBy(duration); // Step duration in seconds double stepT = stepS; // Extrapolation loop //int cpt = 1; for (AbsoluteDate extrapDate = initialDate; extrapDate .compareTo(finalDate) <= 0; extrapDate = extrapDate.shiftedBy(60 * stepT)) { SpacecraftState currentState = propagator.propagate(extrapDate); //System.out.println(/*"step " + */cpt++); //System.out.println(/*" time : " + */currentState.getDate()); //System.out.println(" " + currentState.getOrbit()); KeplerianOrbit o = (KeplerianOrbit) currentState.getOrbit(); String tempSTR; tempSTR = String.format(Locale.US, "%s %12.3f %10.8f %10.6f %10.6f %10.6f %10.6f%n", currentState.getDate(), o.getA(), o.getE(), FastMath.toDegrees(o.getI()), FastMath.toDegrees(o.getPerigeeArgument()), FastMath.toDegrees(o.getRightAscensionOfAscendingNode()), FastMath.toDegrees(o.getTrueAnomaly())); TF.append(tempSTR); } } catch (OrekitException oe) { System.err.println(oe.getMessage()); } }
From source file:fr.cs.examples.bodies.Phasing.java
private void run(final File input) throws IOException, IllegalArgumentException, ParseException, OrekitException { // read input parameters KeyValueFileParser<ParameterKey> parser = new KeyValueFileParser<ParameterKey>(ParameterKey.class); parser.parseInput(new FileInputStream(input)); TimeScale utc = TimeScalesFactory.getUTC(); // simulation properties AbsoluteDate date = parser.getDate(ParameterKey.ORBIT_DATE, utc); int nbOrbits = parser.getInt(ParameterKey.PHASING_ORBITS_NUMBER); int nbDays = parser.getInt(ParameterKey.PHASING_DAYS_NUMBER); double latitude = parser.getAngle(ParameterKey.SUN_SYNCHRONOUS_REFERENCE_LATITUDE); boolean ascending = parser.getBoolean(ParameterKey.SUN_SYNCHRONOUS_REFERENCE_ASCENDING); double mst = parser.getTime(ParameterKey.SUN_SYNCHRONOUS_MEAN_SOLAR_TIME).getSecondsInDay() / 3600; int degree = parser.getInt(ParameterKey.GRAVITY_FIELD_DEGREE); int order = parser.getInt(ParameterKey.GRAVITY_FIELD_ORDER); String gridOutput = parser.getString(ParameterKey.GRID_OUTPUT); double[] gridLatitudes = new double[] { parser.getAngle(ParameterKey.GRID_LATITUDE_1), parser.getAngle(ParameterKey.GRID_LATITUDE_2), parser.getAngle(ParameterKey.GRID_LATITUDE_3), parser.getAngle(ParameterKey.GRID_LATITUDE_4), parser.getAngle(ParameterKey.GRID_LATITUDE_5) }; boolean[] gridAscending = new boolean[] { parser.getBoolean(ParameterKey.GRID_ASCENDING_1), parser.getBoolean(ParameterKey.GRID_ASCENDING_2), parser.getBoolean(ParameterKey.GRID_ASCENDING_3), parser.getBoolean(ParameterKey.GRID_ASCENDING_4), parser.getBoolean(ParameterKey.GRID_ASCENDING_5) }; gravityField = GravityFieldFactory.getNormalizedProvider(degree, order); // initial guess for orbit CircularOrbit orbit = guessOrbit(date, FramesFactory.getEME2000(), nbOrbits, nbDays, latitude, ascending, mst);//from w ww . ja v a2 s . c o m System.out.println("initial orbit: " + orbit); System.out.println("please wait while orbit is adjusted..."); System.out.println(); // numerical model for improving orbit double[][] tolerances = NumericalPropagator.tolerances(0.1, orbit, OrbitType.CIRCULAR); DormandPrince853Integrator integrator = new DormandPrince853Integrator(1.0e-4 * orbit.getKeplerianPeriod(), 1.0e-1 * orbit.getKeplerianPeriod(), tolerances[0], tolerances[1]); integrator.setInitialStepSize(1.0e-2 * orbit.getKeplerianPeriod()); NumericalPropagator propagator = new NumericalPropagator(integrator); propagator.addForceModel(new HolmesFeatherstoneAttractionModel( FramesFactory.getGTOD(IERSConventions.IERS_2010, true), gravityField)); propagator.addForceModel(new ThirdBodyAttraction(CelestialBodyFactory.getSun())); propagator.addForceModel(new ThirdBodyAttraction(CelestialBodyFactory.getMoon())); double deltaP = Double.POSITIVE_INFINITY; double deltaV = Double.POSITIVE_INFINITY; int counter = 0; DecimalFormat f = new DecimalFormat("0.000E00", new DecimalFormatSymbols(Locale.US)); while (deltaP > 3.0e-1 || deltaV > 3.0e-4) { CircularOrbit previous = orbit; CircularOrbit tmp1 = improveEarthPhasing(previous, nbOrbits, nbDays, propagator); CircularOrbit tmp2 = improveSunSynchronization(tmp1, nbOrbits * tmp1.getKeplerianPeriod(), latitude, ascending, mst, propagator); orbit = improveFrozenEccentricity(tmp2, nbOrbits * tmp2.getKeplerianPeriod(), propagator); double da = orbit.getA() - previous.getA(); double dex = orbit.getCircularEx() - previous.getCircularEx(); double dey = orbit.getCircularEy() - previous.getCircularEy(); double di = FastMath.toDegrees(orbit.getI() - previous.getI()); double dr = FastMath.toDegrees( orbit.getRightAscensionOfAscendingNode() - previous.getRightAscensionOfAscendingNode()); System.out.println(" iteration " + (++counter) + ": deltaA = " + f.format(da) + " m, deltaEx = " + f.format(dex) + ", deltaEy = " + f.format(dey) + ", deltaI = " + f.format(di) + " deg, deltaRAAN = " + f.format(dr) + " deg"); PVCoordinates delta = new PVCoordinates(previous.getPVCoordinates(), orbit.getPVCoordinates()); deltaP = delta.getPosition().getNorm(); deltaV = delta.getVelocity().getNorm(); } // final orbit System.out.println(); System.out.println("final orbit (osculating): " + orbit); // generate the ground track grid file PrintStream output = new PrintStream(new File(input.getParent(), gridOutput)); for (int i = 0; i < gridLatitudes.length; ++i) { printGridPoints(output, gridLatitudes[i], gridAscending[i], orbit, propagator, nbOrbits); } output.close(); }
From source file:gbt.ubt.tool.Orthorectifier.java
public static BoundedPropagator generateEphemeris(InputParameters params) throws IOException, OrekitException { /*//from w w w.j av a2s .co m This function generates the satellite ephemeris for the duration of the acquisition using the orbital state vector that is embedded in the product header (MPH). This state vector is normally the restituted orbit from the Flight Operations Segment. The open source Orekit library (https://www.orekit.org/) is used to perform the numerical propagation. The generated ephemeris compares well against the available DORIS Precise Orbit files (maximum 10 metre displacement vector magnitude during propagation duration). */ if (params.orthorectify) { System.out.println("Orthorectification is selected, performing orbit propagation"); /* Get the acquisition duration and state vector from the L1b product */ DataProvidersManager.getInstance().addProvider(new ZipJarCrawler(new File("orekit-data.zip"))); Product readProduct = ProductIO.readProduct(params.inputFileLocation); Date startTime = readProduct.getStartTime().getAsDate(); Date stopTime = readProduct.getEndTime().getAsDate(); MetadataElement metadataRoot = readProduct.getMetadataRoot(); Date vectorTime = metadataRoot.getElementAt(0).getAttributeUTC("STATE_VECTOR_TIME").getAsDate(); Double xPos = metadataRoot.getElementAt(0).getAttributeDouble("X_POSITION"); Double yPos = metadataRoot.getElementAt(0).getAttributeDouble("Y_POSITION"); Double zPos = metadataRoot.getElementAt(0).getAttributeDouble("Z_POSITION"); Double xVel = metadataRoot.getElementAt(0).getAttributeDouble("X_VELOCITY"); Double yVel = metadataRoot.getElementAt(0).getAttributeDouble("Y_VELOCITY"); Double zVel = metadataRoot.getElementAt(0).getAttributeDouble("Z_VELOCITY"); /* Product state vector in fixed frame and UTC Note state vector and product/acquisition start time are not coincident */ Frame fixedFrame = FramesFactory.getITRF(IERSConventions.IERS_2010, false); TimeScale utc = TimeScalesFactory.getUTC(); GregorianCalendar cal = new GregorianCalendar(TimeZone.getTimeZone("UTC")); cal.setTime(vectorTime); AbsoluteDate initialDate = new AbsoluteDate(cal.get(Calendar.YEAR), cal.get(Calendar.MONTH) + 1, cal.get(Calendar.DATE), cal.get(Calendar.HOUR_OF_DAY), cal.get(Calendar.MINUTE), cal.get(Calendar.SECOND), utc); Vector3D position = new Vector3D(xPos, yPos, zPos); Vector3D velocity = new Vector3D(xVel, yVel, zVel); /* Orekit can only propate orbits defined in inertial reference frame Convert state vector to J2000 (EME2000) frame. */ Frame inertialFrame = FramesFactory.getEME2000(); Transform frameTransform = fixedFrame.getTransformTo(inertialFrame, initialDate); PVCoordinates transformPVCoordinates = frameTransform .transformPVCoordinates(new PVCoordinates(position, velocity)); /* Set initial spacecraft state */ double mu = 3.986004415e+14; // Central attraction coefficient Orbit initialOrbit = new CartesianOrbit(transformPVCoordinates, inertialFrame, initialDate, mu); SpacecraftState initialState = new SpacecraftState(initialOrbit, 7892.0); // Orbital parameters and Envisat dry mass /* Set up numerical integrator for propagation */ double minStep = 1; double maxStep = 1000; double initialStep = 60; OrbitType propagationType = OrbitType.CARTESIAN; double[][] tolerance = NumericalPropagator.tolerances(0.001, initialOrbit, propagationType); AdaptiveStepsizeIntegrator integrator = new DormandPrince853Integrator(minStep, maxStep, tolerance[0], tolerance[1]); integrator.setInitialStepSize(initialStep); NumericalPropagator propagator = new NumericalPropagator(integrator); propagator.setOrbitType(propagationType); /* Add force models to the propagator: Consider central body (Earth) gravity pertubation Consider third body (Sun & Moon) gravity pertubation Other pertubations are not considered (computation speed/accuracy trade-off) */ NormalizedSphericalHarmonicsProvider provider = GravityFieldFactory.getNormalizedProvider(21, 21); ForceModel holmesFeatherstone = new HolmesFeatherstoneAttractionModel( FramesFactory.getITRF(IERSConventions.IERS_2010, true), provider); propagator.addForceModel(holmesFeatherstone); propagator.addForceModel(new ThirdBodyAttraction(CelestialBodyFactory.getSun())); propagator.addForceModel(new ThirdBodyAttraction(CelestialBodyFactory.getMoon())); /* Run the propagator to generate the ephemeris */ propagator.setInitialState(initialState); propagator.setEphemerisMode(); cal.setTime(startTime); AbsoluteDate startDate = new AbsoluteDate(cal.get(Calendar.YEAR), cal.get(Calendar.MONTH) + 1, cal.get(Calendar.DATE), cal.get(Calendar.HOUR_OF_DAY), cal.get(Calendar.MINUTE), cal.get(Calendar.SECOND), utc); cal.setTime(stopTime); AbsoluteDate stopDate = new AbsoluteDate(cal.get(Calendar.YEAR), cal.get(Calendar.MONTH) + 1, cal.get(Calendar.DATE), cal.get(Calendar.HOUR_OF_DAY), cal.get(Calendar.MINUTE), cal.get(Calendar.SECOND), utc); SpacecraftState finalState = propagator.propagate(startDate.shiftedBy(-600.0), stopDate.shiftedBy(600.0)); BoundedPropagator ephemeris = propagator.getGeneratedEphemeris(); return ephemeris; } else { return null; } }
From source file:fr.cs.examples.propagation.DSSTPropagation.java
/** Set up the DSST Propagator * * @param orbit initial orbit/*from w w w.j a va 2 s . com*/ * @param mass S/C mass (kg) * @param isOsculating if orbital elements are osculating * @param fixedStepSize step size for fixed step integrator (s) * @throws OrekitException */ private DSSTPropagator createDSSTProp(final Orbit orbit, final double mass, final boolean isOsculating, final double fixedStepSize) throws OrekitException { AbstractIntegrator integrator; if (fixedStepSize > 0.) { integrator = new ClassicalRungeKuttaIntegrator(fixedStepSize); } else { final double minStep = orbit.getKeplerianPeriod(); final double maxStep = minStep * 100.; final double[][] tol = DSSTPropagator.tolerances(1.0, orbit); integrator = new DormandPrince853Integrator(minStep, maxStep, tol[0], tol[1]); ((AdaptiveStepsizeIntegrator) integrator).setInitialStepSize(10. * minStep); } DSSTPropagator dsstProp = new DSSTPropagator(integrator); dsstProp.setInitialState(new SpacecraftState(orbit, mass), isOsculating); return dsstProp; }
From source file:fr.cs.examples.propagation.DSSTPropagation.java
/** Create the numerical propagator * * @param orbit initial orbit//from ww w . j a v a 2 s . c om * @param mass S/C mass (kg) * @throws OrekitException */ private NumericalPropagator createNumProp(final Orbit orbit, final double mass) throws OrekitException { final double[][] tol = NumericalPropagator.tolerances(1.0, orbit, orbit.getType()); final double minStep = 1.e-3; final double maxStep = 1.e+3; AdaptiveStepsizeIntegrator integrator = new DormandPrince853Integrator(minStep, maxStep, tol[0], tol[1]); integrator.setInitialStepSize(100.); NumericalPropagator numProp = new NumericalPropagator(integrator); numProp.setInitialState(new SpacecraftState(orbit, mass)); return numProp; }
From source file:nl.rivm.cib.episim.model.disease.infection.MSEIRSPlot.java
@Override public void start(final Stage stage) { final SIRConfig conf = ConfigFactory.create(SIRConfig.class); final double[] t = conf.t(); final long[] pop = conf.population(); final double n0 = Arrays.stream(pop).sum(); final String[] colors = conf.colors(), colors2 = conf.colors2(); final Pane plot = new Pane(); plot.setPrefSize(400, 300);/*w ww. j a v a 2 s .c om*/ plot.setMinSize(50, 50); final NumberAxis xAxis = new NumberAxis(t[0], t[1], (t[1] - t[0]) / 10); final NumberAxis yAxis = new NumberAxis(0, n0, n0 / 10); final Pane axes = new Pane(); axes.prefHeightProperty().bind(plot.heightProperty()); axes.prefWidthProperty().bind(plot.widthProperty()); xAxis.setSide(Side.BOTTOM); xAxis.setMinorTickVisible(false); xAxis.setPrefWidth(axes.getPrefWidth()); xAxis.prefWidthProperty().bind(axes.widthProperty()); xAxis.layoutYProperty().bind(axes.heightProperty()); yAxis.setSide(Side.LEFT); yAxis.setMinorTickVisible(false); yAxis.setPrefHeight(axes.getPrefHeight()); yAxis.prefHeightProperty().bind(axes.heightProperty()); yAxis.layoutXProperty().bind(Bindings.subtract(1, yAxis.widthProperty())); axes.getChildren().setAll(xAxis, yAxis); final Label lbl = new Label(String.format("R0=%.1f, recovery=%.1ft\nSIR(0)=%s", conf.reproduction(), conf.recovery(), Arrays.toString(pop))); lbl.setTextAlignment(TextAlignment.CENTER); lbl.setTextFill(Color.WHITE); final Path[] deterministic = { new Path(), new Path(), new Path() }; IntStream.range(0, pop.length).forEach(i -> { final Color color = Color.valueOf(colors[i]); final Path path = deterministic[i]; path.setStroke(color.deriveColor(0, 1, 1, 0.6)); path.setStrokeWidth(2); path.setClip(new Rectangle(0, 0, plot.getPrefWidth(), plot.getPrefHeight())); }); plot.getChildren().setAll(axes); // fill paths with integration estimates final double xl = xAxis.getLowerBound(), sx = plot.getPrefWidth() / (xAxis.getUpperBound() - xl), yh = plot.getPrefHeight(), sy = yh / (yAxis.getUpperBound() - yAxis.getLowerBound()); final TreeMap<Double, Integer> iDeterministic = new TreeMap<>(); MSEIRSTest.deterministic(conf, () -> new DormandPrince853Integrator(1.0E-8, 10, 1.0E-20, 1.0E-20)) .subscribe(yt -> { iDeterministic.put(yt.getKey(), deterministic[0].getElements().size()); final double[] y = yt.getValue(); final double x = (yt.getKey() - xl) * sx; for (int i = 0; i < y.length; i++) { final double yi = yh - y[i] * sy; final PathElement di = deterministic[i].getElements().isEmpty() ? new MoveTo(x, yi) : new LineTo(x, yi); deterministic[i].getElements().add(di); } }, e -> LOG.error("Problem", e), () -> plot.getChildren().addAll(deterministic)); final Path[] stochasticTau = { new Path(), new Path(), new Path() }; IntStream.range(0, pop.length).forEach(i -> { final Color color = Color.valueOf(colors[i]); final Path path = stochasticTau[i]; path.setStroke(color); path.setStrokeWidth(1); path.setClip(new Rectangle(0, 0, plot.getPrefWidth(), plot.getPrefHeight())); }); final TreeMap<Double, Integer> iStochasticTau = new TreeMap<>(); MSEIRSTest.stochasticGillespie(conf).subscribe(yt -> { final double x = (yt.getKey() - xl) * sx; iStochasticTau.put(yt.getKey(), stochasticTau[0].getElements().size()); final long[] y = yt.getValue(); for (int i = 0; i < y.length; i++) { final double yi = yh - y[i] * sy; final ObservableList<PathElement> path = stochasticTau[i].getElements(); if (path.isEmpty()) { path.add(new MoveTo(x, yi)); // first } else { final PathElement last = path.get(path.size() - 1); final double y_prev = last instanceof MoveTo ? ((MoveTo) last).getY() : ((LineTo) last).getY(); path.add(new LineTo(x, y_prev)); path.add(new LineTo(x, yi)); } } }, e -> LOG.error("Problem", e), () -> plot.getChildren().addAll(stochasticTau)); final Path[] stochasticRes = { new Path(), new Path(), new Path() }; IntStream.range(0, pop.length).forEach(i -> { final Color color = Color.valueOf(colors2[i]); final Path path = stochasticRes[i]; path.setStroke(color); path.setStrokeWidth(1); path.setClip(new Rectangle(0, 0, plot.getPrefWidth(), plot.getPrefHeight())); }); final TreeMap<Double, Integer> iStochasticRes = new TreeMap<>(); MSEIRSTest.stochasticSellke(conf).subscribe(yt -> { final double x = (yt.getKey() - xl) * sx; iStochasticRes.put(yt.getKey(), stochasticRes[0].getElements().size()); final long[] y = yt.getValue(); for (int i = 0; i < y.length; i++) { final double yi = yh - y[i] * sy; final ObservableList<PathElement> path = stochasticRes[i].getElements(); if (path.isEmpty()) { path.add(new MoveTo(x, yi)); // first } else { final PathElement last = path.get(path.size() - 1); final double y_prev = last instanceof MoveTo ? ((MoveTo) last).getY() : ((LineTo) last).getY(); path.add(new LineTo(x, y_prev)); path.add(new LineTo(x, yi)); } } }, e -> LOG.error("Problem", e), () -> plot.getChildren().addAll(stochasticRes)); // auto-scale on stage/plot resize // FIXME scaling around wrong origin, use ScatterChart? // xAxis.widthProperty() // .addListener( (ChangeListener<Number>) ( observable, // oldValue, newValue ) -> // { // final double scale = ((Double) newValue) // / plot.getPrefWidth(); // plot.getChildren().filtered( n -> n instanceof Path ) // .forEach( n -> // { // final Path path = (Path) n; // path.setScaleX( scale ); // path.setTranslateX( (path // .getBoundsInParent().getWidth() // - path.getLayoutBounds().getWidth()) // / 2 ); // } ); // } ); // plot.heightProperty() // .addListener( (ChangeListener<Number>) ( observable, // oldValue, newValue ) -> // { // final double scale = ((Double) newValue) // / plot.getPrefHeight(); // plot.getChildren().filtered( n -> n instanceof Path ) // .forEach( n -> // { // final Path path = (Path) n; // path.setScaleY( scale ); // path.setTranslateY( // (path.getBoundsInParent() // .getHeight() * (scale - 1)) // / 2 ); // } ); // } ); final StackPane layout = new StackPane(lbl, plot); layout.setAlignment(Pos.TOP_CENTER); layout.setPadding(new Insets(50)); layout.setStyle("-fx-background-color: rgb(35, 39, 50);"); final Line vertiCross = new Line(); vertiCross.setStroke(Color.SILVER); vertiCross.setStrokeWidth(1); vertiCross.setVisible(false); axes.getChildren().add(vertiCross); final Tooltip tip = new Tooltip(""); tip.setAutoHide(false); tip.hide(); axes.setOnMouseExited(ev -> tip.hide()); axes.setOnMouseMoved(ev -> { final Double x = (Double) xAxis.getValueForDisplay(ev.getX()); if (x > xAxis.getUpperBound() || x < xAxis.getLowerBound()) { tip.hide(); vertiCross.setVisible(false); return; } final Double y = (Double) yAxis.getValueForDisplay(ev.getY()); if (y > yAxis.getUpperBound() || y < yAxis.getLowerBound()) { tip.hide(); vertiCross.setVisible(false); return; } final double xs = xAxis.getDisplayPosition(x); vertiCross.setStartX(xs); vertiCross.setStartY(yAxis.getDisplayPosition(0)); vertiCross.setEndX(xs); vertiCross.setEndY(yAxis.getDisplayPosition(yAxis.getUpperBound())); vertiCross.setVisible(true); final int i = (iDeterministic.firstKey() > x ? iDeterministic.firstEntry() : iDeterministic.floorEntry(x)).getValue(); final Object[] yi = Arrays.stream(deterministic).mapToDouble(p -> getY(p, i)) .mapToObj(yAxis::getValueForDisplay).map(n -> DecimalUtil.toScale(n, 1)).toArray(); final int j = (iStochasticTau.firstKey() > x ? iStochasticTau.firstEntry() : iStochasticTau.floorEntry(x)).getValue(); final Object[] yj = Arrays.stream(stochasticTau).mapToDouble(p -> getY(p, j)) .mapToObj(yAxis::getValueForDisplay).map(n -> DecimalUtil.toScale(n, 0)).toArray(); final int k = (iStochasticRes.firstKey() > x ? iStochasticRes.firstEntry() : iStochasticRes.floorEntry(x)).getValue(); final Object[] yk = Arrays.stream(stochasticRes).mapToDouble(p -> getY(p, k)) .mapToObj(yAxis::getValueForDisplay).map(n -> DecimalUtil.toScale(n, 0)).toArray(); final String txt = String.format("SIR(t=%.1f)\n" + "~det%s\n" + "~tau%s\n" + "~res%s", x, Arrays.toString(yi), Arrays.toString(yj), Arrays.toString(yk)); tip.setText(txt); tip.show(axes, ev.getScreenX() - ev.getSceneX() + xs, ev.getScreenY() + 15); }); try { stage.getIcons().add(new Image(FileUtil.toInputStream("icon.jpg"))); } catch (final IOException e) { LOG.error("Problem", e); } stage.setTitle("Deterministic vs. Stochastic"); stage.setScene(new Scene(layout, Color.rgb(35, 39, 50))); // stage.setOnHidden( ev -> tip.hide() ); stage.show(); }
From source file:nl.rivm.cib.episim.model.disease.infection.MSEIRSTest.java
@Test public void compareMSEIRS() throws TimeoutException { final SIRConfig conf = ConfigFactory.create(SIRConfig.class); final double[] T = conf.t(); final long[] pop = conf.population(); LOG.info("Starting with SIR: {}, t in {}", pop, T); final int t = (int) T[0], N = (int) T[1], C = pop.length; final MatrixBuilder result = MatrixBuilder.sparse(N - t + 1, pop.length * 4).label("Results") .labelColumns("~dS.dt", "~dI.dt", "~dR.dt", "S_gil(t)", "I_gil(t)", "R_gil(t)", "S_sel(t)", "I_sel(t)", "R_sel(t)", "S_dyn(t)", "I_dyn(t)", "R_dyn(t)") .labelRows(i -> "max(t=<" + i + ")"); // the Dormand-Prince (embedded Runge-Kutta) ODE integrator // see https://www.wikiwand.com/en/Runge%E2%80%93Kutta_methods deterministic(conf, () -> new DormandPrince853Integrator(1.0E-8, 10, 1.0E-20, 1.0E-20)) .blockingForEach(yt -> result.withContent(Arrays.stream(yt.getValue()).mapToObj(d -> d), DecimalUtil.ceil(yt.getKey()).longValue(), 0 * C)); final int n = 10; final Waiter waiter = new Waiter(); LOG.trace("Running Gillespie x {}...", n); // see https://www.wikiwand.com/en/Gillespie_algorithm averages(() -> stochasticGillespie(conf, 1d), d -> DecimalUtil.ceil(d).longValue(), n) .forEach(stats -> result.withContent(stats.getValue(), stats.getKey(), 1 * C)); LOG.trace("Running Sellke x {}...", n); // see e.g. Cook et al. 2008, equation 3.5: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3227033/#__sec2title averages(() -> stochasticSellke(conf, 1d), d -> DecimalUtil.ceil(d).longValue(), n) .forEach(stats -> result.withContent(stats.getValue(), stats.getKey(), 2 * C)); LOG.trace("Running Epidemes x {}...", 1); averages(() -> stochasticLocal(conf, 1d), d -> DecimalUtil.ceil(d).longValue(), 1).subscribe( stats -> result.withContent(stats.getValue(), stats.getKey(), 3 * C), waiter::fail, waiter::resume); //waiter.await(); // Application.launch( CartesianPlot.class ); LOG.trace("result: \n{}", result.build().toString()); LOG.info("Complete"); }
From source file:org.orekit.bodies.SolarBodyTest.java
@Test public void testPropagationVsEphemeris() throws OrekitException { Utils.setDataRoot("regular-data"); //Creation of the celestial bodies of the solar system final CelestialBody sun = CelestialBodyFactory.getSun(); final CelestialBody mercury = CelestialBodyFactory.getMercury(); final CelestialBody venus = CelestialBodyFactory.getVenus(); final CelestialBody earth = CelestialBodyFactory.getEarth(); final CelestialBody mars = CelestialBodyFactory.getMars(); final CelestialBody jupiter = CelestialBodyFactory.getJupiter(); final CelestialBody saturn = CelestialBodyFactory.getSaturn(); final CelestialBody uranus = CelestialBodyFactory.getUranus(); final CelestialBody neptune = CelestialBodyFactory.getNeptune(); final CelestialBody pluto = CelestialBodyFactory.getPluto(); //Starting and end dates final AbsoluteDate startingDate = new AbsoluteDate(2000, 1, 2, TimeScalesFactory.getUTC()); AbsoluteDate endDate = startingDate.shiftedBy(30 * Constants.JULIAN_DAY); final Frame icrf = FramesFactory.getICRF(); // fake orbit around negligible point mass at solar system barycenter double negligibleMu = 1.0e-3; SpacecraftState initialState = new SpacecraftState( new CartesianOrbit(venus.getPVCoordinates(startingDate, icrf), icrf, startingDate, negligibleMu)); //Creation of the numerical propagator final double[][] tol = NumericalPropagator.tolerances(1000, initialState.getOrbit(), OrbitType.CARTESIAN); AbstractIntegrator dop1 = new DormandPrince853Integrator(1.0, 1.0e5, tol[0], tol[1]); NumericalPropagator propag = new NumericalPropagator(dop1); propag.setOrbitType(OrbitType.CARTESIAN); propag.setInitialState(initialState); propag.setMu(negligibleMu);//from w w w .j a va2 s. c o m //Creation of the ForceModels propag.addForceModel(new BodyAttraction(sun)); propag.addForceModel(new BodyAttraction(mercury)); propag.addForceModel(new BodyAttraction(earth)); propag.addForceModel(new BodyAttraction(mars)); propag.addForceModel(new BodyAttraction(jupiter)); propag.addForceModel(new BodyAttraction(saturn)); propag.addForceModel(new BodyAttraction(uranus)); propag.addForceModel(new BodyAttraction(neptune)); propag.addForceModel(new BodyAttraction(pluto)); // checks are done within the step handler propag.setMasterMode(1000.0, new OrekitFixedStepHandler() { public void init(SpacecraftState s0, AbsoluteDate t) { } public void handleStep(SpacecraftState currentState, boolean isLast) throws PropagationException { try { // propagated position should remain within 1400m of ephemeris for one month Vector3D propagatedP = currentState.getPVCoordinates(icrf).getPosition(); Vector3D ephemerisP = venus.getPVCoordinates(currentState.getDate(), icrf).getPosition(); Assert.assertEquals(0, Vector3D.distance(propagatedP, ephemerisP), 1400.0); } catch (OrekitException oe) { throw new PropagationException(oe); } } }); propag.propagate(startingDate, endDate); }
From source file:org.orekit.forces.drag.DragForceTest.java
@Test public void testStateJacobianSphere() throws OrekitException { // initialization AbsoluteDate date = new AbsoluteDate(new DateComponents(2003, 03, 01), new TimeComponents(13, 59, 27.816), TimeScalesFactory.getUTC()); double i = FastMath.toRadians(98.7); double omega = FastMath.toRadians(93.0); double OMEGA = FastMath.toRadians(15.0 * 22.5); Orbit orbit = new KeplerianOrbit(7201009.7124401, 1e-3, i, omega, OMEGA, 0, PositionAngle.MEAN, FramesFactory.getEME2000(), date, Constants.EIGEN5C_EARTH_MU); OrbitType integrationType = OrbitType.CARTESIAN; double[][] tolerances = NumericalPropagator.tolerances(0.01, orbit, integrationType); NumericalPropagator propagator = new NumericalPropagator( new DormandPrince853Integrator(1.0e-3, 120, tolerances[0], tolerances[1])); propagator.setOrbitType(integrationType); final DragForce forceModel = new DragForce(new HarrisPriester(CelestialBodyFactory.getSun(), new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, Constants.WGS84_EARTH_FLATTENING, FramesFactory.getITRF(IERSConventions.IERS_2010, true))), new SphericalSpacecraft(2.5, 1.2, 0.7, 0.2)); propagator.addForceModel(forceModel); SpacecraftState state0 = new SpacecraftState(orbit); checkStateJacobian(propagator, state0, date.shiftedBy(3.5 * 3600.0), 1e3, tolerances[0], 3.0e-4); }
From source file:org.orekit.forces.drag.DragForceTest.java
@Test public void testStateJacobianBox() throws OrekitException { // initialization AbsoluteDate date = new AbsoluteDate(new DateComponents(2003, 03, 01), new TimeComponents(13, 59, 27.816), TimeScalesFactory.getUTC()); double i = FastMath.toRadians(98.7); double omega = FastMath.toRadians(93.0); double OMEGA = FastMath.toRadians(15.0 * 22.5); Orbit orbit = new KeplerianOrbit(7201009.7124401, 1e-3, i, omega, OMEGA, 0, PositionAngle.MEAN, FramesFactory.getEME2000(), date, Constants.EIGEN5C_EARTH_MU); OrbitType integrationType = OrbitType.CARTESIAN; double[][] tolerances = NumericalPropagator.tolerances(0.01, orbit, integrationType); NumericalPropagator propagator = new NumericalPropagator( new DormandPrince853Integrator(1.0e-3, 120, tolerances[0], tolerances[1])); propagator.setOrbitType(integrationType); final DragForce forceModel = new DragForce(new HarrisPriester(CelestialBodyFactory.getSun(), new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, Constants.WGS84_EARTH_FLATTENING, FramesFactory.getITRF(IERSConventions.IERS_2010, true))), new BoxAndSolarArraySpacecraft(1.5, 2.0, 1.8, CelestialBodyFactory.getSun(), 20.0, Vector3D.PLUS_J, 1.2, 0.7, 0.2));//from w w w . ja va2 s. c om propagator.addForceModel(forceModel); SpacecraftState state0 = new SpacecraftState(orbit); checkStateJacobian(propagator, state0, date.shiftedBy(3.5 * 3600.0), 1e3, tolerances[0], 2.0e-3); }