List of usage examples for org.apache.commons.math3.analysis.differentiation DerivativeStructure getValue
public double getValue()
From source file:com.bwc.ora.models.Lrp.java
/** * Get the local maximums from a collection of Points. * * @param lrpSeries/* ww w .j a v a2s . c om*/ * @param seriesTitle * @return */ public static XYSeries getMaximumsWithHiddenPeaks(XYSeries lrpSeries, String seriesTitle) { XYSeries maxPoints = new XYSeries(seriesTitle); //convert to x and y coordinate arrays double[][] xyline = lrpSeries.toArray(); //use a spline interpolator to converts points into an equation UnivariateInterpolator interpolator = new SplineInterpolator(); UnivariateFunction function = interpolator.interpolate(xyline[0], xyline[1]); // create a differentiator using 5 points and 0.01 step FiniteDifferencesDifferentiator differentiator = new FiniteDifferencesDifferentiator(5, 0.01); // create a new function that computes both the value and the derivatives // using DerivativeStructure UnivariateDifferentiableFunction completeF = differentiator.differentiate(function); // now we can compute the value and its derivatives // here we decided to display up to second order derivatives, // because we feed completeF with order 2 DerivativeStructure instances //find local minima in second derivative, these indicate the peaks (and hidden peaks) //of the input for (double x = xyline[0][0] + 1; x < xyline[0][xyline[0].length - 1] - 1; x += 0.5) { DerivativeStructure xDSc = new DerivativeStructure(1, 2, 0, x); DerivativeStructure xDSl = new DerivativeStructure(1, 2, 0, x - 0.5); DerivativeStructure xDSr = new DerivativeStructure(1, 2, 0, x + 0.5); DerivativeStructure yDSc = completeF.value(xDSc); DerivativeStructure yDSl = completeF.value(xDSl); DerivativeStructure yDSr = completeF.value(xDSr); double c2d = yDSc.getPartialDerivative(2); if (c2d < yDSl.getPartialDerivative(2) && c2d < yDSr.getPartialDerivative(2)) { maxPoints.add((int) Math.round(x), yDSc.getValue()); } } return maxPoints; }
From source file:de.tuberlin.uebb.jdae.examples.FactorialEquation.java
private final DerivativeStructure dfactorial(final DerivativeStructure ds) { if (ds.getValue() <= 1.0) { final DerivativeStructure divide = ds.divide(ds); return divide; } else// www. ja v a 2 s . com return ds.multiply(dfactorial(ds.subtract(1.0))); }
From source file:com.itemanalysis.psychometrics.irt.estimation.IrtExaminee.java
/** * Returns value of loglikelihood using DerivativeStructure per interface requirements * @param t//from w ww .j a v a 2 s.c o m * @return */ public DerivativeStructure value(DerivativeStructure t) { return new DerivativeStructure(1, 0, 0, logLikelihood(t.getValue())); }
From source file:com.itemanalysis.psychometrics.irt.estimation.IrtExaminee.java
/** * Returns first derivative of loglikelihood using DerivativeStructure per interface requirements * @param t/*from w w w . java 2 s .co m*/ * @return */ public DerivativeStructure derivLogLikelihood(DerivativeStructure t) { return new DerivativeStructure(1, 1, 0, derivLogLikelihood(t.getValue())); }
From source file:org.orekit.bodies.OneAxisEllipsoidTest.java
@Test public void testMovingGeodeticPoint() throws OrekitException { final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, Constants.WGS84_EARTH_FLATTENING, FramesFactory.getITRF(IERSConventions.IERS_2010, true)); double lat0 = FastMath.toRadians(60.0); double lon0 = FastMath.toRadians(25.0); double alt0 = 100.0; double lat1 = 1.0e-3; double lon1 = -2.0e-3; double alt1 = 1.2; double lat2 = -1.0e-5; double lon2 = -3.0e-5; double alt2 = -0.01; final DerivativeStructure latDS = new DerivativeStructure(1, 2, lat0, lat1, lat2); final DerivativeStructure lonDS = new DerivativeStructure(1, 2, lon0, lon1, lon2); final DerivativeStructure altDS = new DerivativeStructure(1, 2, alt0, alt1, alt2); // direct computation of position, velocity and acceleration PVCoordinates pv = earth.transform(new FieldGeodeticPoint<DerivativeStructure>(latDS, lonDS, altDS)); // finite differences computation FiniteDifferencesDifferentiator differentiator = new FiniteDifferencesDifferentiator(5, 0.1); UnivariateDifferentiableFunction fx = differentiator.differentiate(new UnivariateFunction() { public double value(double dt) { GeodeticPoint gp = new GeodeticPoint(latDS.taylor(dt), lonDS.taylor(dt), altDS.taylor(dt)); return earth.transform(gp).getX(); }/*from w w w . java2s . co m*/ }); UnivariateDifferentiableFunction fy = differentiator.differentiate(new UnivariateFunction() { public double value(double dt) { GeodeticPoint gp = new GeodeticPoint(latDS.taylor(dt), lonDS.taylor(dt), altDS.taylor(dt)); return earth.transform(gp).getY(); } }); UnivariateDifferentiableFunction fz = differentiator.differentiate(new UnivariateFunction() { public double value(double dt) { GeodeticPoint gp = new GeodeticPoint(latDS.taylor(dt), lonDS.taylor(dt), altDS.taylor(dt)); return earth.transform(gp).getZ(); } }); DerivativeStructure dtZero = new DerivativeStructure(1, 2, 0, 0.0); DerivativeStructure xDS = fx.value(dtZero); DerivativeStructure yDS = fy.value(dtZero); DerivativeStructure zDS = fz.value(dtZero); Assert.assertEquals(xDS.getValue(), pv.getPosition().getX(), 2.0e-20 * FastMath.abs(xDS.getValue())); Assert.assertEquals(xDS.getPartialDerivative(1), pv.getVelocity().getX(), 2.0e-12 * FastMath.abs(xDS.getPartialDerivative(1))); Assert.assertEquals(xDS.getPartialDerivative(2), pv.getAcceleration().getX(), 2.0e-9 * FastMath.abs(xDS.getPartialDerivative(2))); Assert.assertEquals(yDS.getValue(), pv.getPosition().getY(), 2.0e-20 * FastMath.abs(yDS.getValue())); Assert.assertEquals(yDS.getPartialDerivative(1), pv.getVelocity().getY(), 2.0e-12 * FastMath.abs(yDS.getPartialDerivative(1))); Assert.assertEquals(yDS.getPartialDerivative(2), pv.getAcceleration().getY(), 2.0e-9 * FastMath.abs(yDS.getPartialDerivative(2))); Assert.assertEquals(zDS.getValue(), pv.getPosition().getZ(), 2.0e-20 * FastMath.abs(zDS.getValue())); Assert.assertEquals(zDS.getPartialDerivative(1), pv.getVelocity().getZ(), 2.0e-12 * FastMath.abs(zDS.getPartialDerivative(1))); Assert.assertEquals(zDS.getPartialDerivative(2), pv.getAcceleration().getZ(), 2.0e-9 * FastMath.abs(zDS.getPartialDerivative(2))); }
From source file:org.orekit.forces.BoxAndSolarArraySpacecraft.java
/** Get solar array normal in spacecraft frame. * @param date current date// w ww . jav a 2s. c o m * @param frame inertial reference frame for state (both orbit and attitude) * @param position position of spacecraft in reference frame * @param rotation orientation (attitude) of the spacecraft with respect to reference frame * @return solar array normal in spacecraft frame * @exception OrekitException if sun direction cannot be computed in best lightning * configuration */ public synchronized FieldVector3D<DerivativeStructure> getNormal(final AbsoluteDate date, final Frame frame, final FieldVector3D<DerivativeStructure> position, final FieldRotation<DerivativeStructure> rotation) throws OrekitException { final DerivativeStructure one = position.getX().getField().getOne(); if (referenceDate != null) { // use a simple rotation at fixed rate final DerivativeStructure alpha = one.multiply(rotationRate * date.durationFrom(referenceDate)); return new FieldVector3D<DerivativeStructure>(alpha.cos(), saX, alpha.sin(), saY); } // compute orientation for best lightning final FieldVector3D<DerivativeStructure> sunInert = position .subtract(sun.getPVCoordinates(date, frame).getPosition()).negate().normalize(); final FieldVector3D<DerivativeStructure> sunSpacecraft = rotation.applyTo(sunInert); final DerivativeStructure d = FieldVector3D.dotProduct(sunSpacecraft, saZ); final DerivativeStructure f = d.multiply(d).subtract(1).negate(); if (f.getValue() < Precision.EPSILON) { // extremely rare case: the sun is along solar array rotation axis // (there will not be much output power ...) // we set up an arbitrary normal return new FieldVector3D<DerivativeStructure>(one, saZ.orthogonal()); } final DerivativeStructure s = f.sqrt().reciprocal(); return new FieldVector3D<DerivativeStructure>(s, sunSpacecraft) .subtract(new FieldVector3D<DerivativeStructure>(s.multiply(d), saZ)); }
From source file:org.orekit.forces.BoxAndSolarArraySpacecraft.java
/** {@inheritDoc} */ public FieldVector3D<DerivativeStructure> dragAcceleration(final AbsoluteDate date, final Frame frame, final FieldVector3D<DerivativeStructure> position, final FieldRotation<DerivativeStructure> rotation, final DerivativeStructure mass, final double density, final FieldVector3D<DerivativeStructure> relativeVelocity) throws OrekitException { // relative velocity in spacecraft frame final FieldVector3D<DerivativeStructure> v = rotation.applyTo(relativeVelocity); // solar array contribution final FieldVector3D<DerivativeStructure> solarArrayFacet = new FieldVector3D<DerivativeStructure>( solarArrayArea, getNormal(date, frame, position, rotation)); DerivativeStructure sv = FieldVector3D.dotProduct(v, solarArrayFacet).abs(); // body facets contribution for (final Facet facet : facets) { final DerivativeStructure dot = FieldVector3D.dotProduct(v, facet.getNormal()); if (dot.getValue() < 0) { // the facet intercepts the incoming flux sv = sv.subtract(dot.multiply(facet.getArea())); }//from ww w .ja v a2 s . c om } return new FieldVector3D<DerivativeStructure>(sv.multiply(density * dragCoeff / 2.0).divide(mass), relativeVelocity); }
From source file:org.orekit.forces.BoxAndSolarArraySpacecraft.java
/** {@inheritDoc} */ public FieldVector3D<DerivativeStructure> radiationPressureAcceleration(final AbsoluteDate date, final Frame frame, final FieldVector3D<DerivativeStructure> position, final FieldRotation<DerivativeStructure> rotation, final DerivativeStructure mass, final FieldVector3D<DerivativeStructure> flux) throws OrekitException { if (flux.getNormSq().getValue() < Precision.SAFE_MIN) { // null illumination (we are probably in umbra) return new FieldVector3D<DerivativeStructure>(0.0, flux); }/* w w w .j a v a 2s . c o m*/ // radiation flux in spacecraft frame final FieldVector3D<DerivativeStructure> fluxSat = rotation.applyTo(flux); // solar array contribution FieldVector3D<DerivativeStructure> normal = getNormal(date, frame, position, rotation); DerivativeStructure dot = FieldVector3D.dotProduct(normal, fluxSat); if (dot.getValue() > 0) { // the solar array is illuminated backward, // fix signs to compute contribution correctly dot = dot.negate(); normal = normal.negate(); } FieldVector3D<DerivativeStructure> force = facetRadiationAcceleration(normal, solarArrayArea, fluxSat, dot); // body facets contribution for (final Facet bodyFacet : facets) { normal = new FieldVector3D<DerivativeStructure>(mass.getField().getOne(), bodyFacet.getNormal()); dot = FieldVector3D.dotProduct(normal, fluxSat); if (dot.getValue() < 0) { // the facet intercepts the incoming flux force = force.add(facetRadiationAcceleration(normal, bodyFacet.getArea(), fluxSat, dot)); } } // convert to inertial frame return rotation.applyInverseTo(new FieldVector3D<DerivativeStructure>(mass.reciprocal(), force)); }
From source file:org.orekit.propagation.analytical.EcksteinHechlerPropagator.java
/** Computes the eccentric latitude argument from the mean latitude argument. * @param alphaM = M + mean latitude argument (rad) * @param ex e cos(), first component of circular eccentricity vector * @param ey e sin(), second component of circular eccentricity vector * @return the eccentric latitude argument. */// w w w .j av a 2 s. c o m private DerivativeStructure meanToEccentric(final DerivativeStructure alphaM, final DerivativeStructure ex, final DerivativeStructure ey) { // Generalization of Kepler equation to circular parameters // with alphaE = PA + E and // alphaM = PA + M = alphaE - ex.sin(alphaE) + ey.cos(alphaE) DerivativeStructure alphaE = alphaM; DerivativeStructure shift = alphaM.getField().getZero(); DerivativeStructure alphaEMalphaM = alphaM.getField().getZero(); DerivativeStructure cosAlphaE = alphaE.cos(); DerivativeStructure sinAlphaE = alphaE.sin(); int iter = 0; do { final DerivativeStructure f2 = ex.multiply(sinAlphaE).subtract(ey.multiply(cosAlphaE)); final DerivativeStructure f1 = alphaM.getField().getOne().subtract(ex.multiply(cosAlphaE)) .subtract(ey.multiply(sinAlphaE)); final DerivativeStructure f0 = alphaEMalphaM.subtract(f2); final DerivativeStructure f12 = f1.multiply(2); shift = f0.multiply(f12).divide(f1.multiply(f12).subtract(f0.multiply(f2))); alphaEMalphaM = alphaEMalphaM.subtract(shift); alphaE = alphaM.add(alphaEMalphaM); cosAlphaE = alphaE.cos(); sinAlphaE = alphaE.sin(); } while ((++iter < 50) && (FastMath.abs(shift.getValue()) > 1.0e-12)); return alphaE; }
From source file:org.orekit.propagation.numerical.Jacobianizer.java
/** Compute acceleration and derivatives with respect to state. * @param date current date// w w w .jav a 2 s . c o m * @param frame inertial reference frame for state (both orbit and attitude) * @param position position of spacecraft in reference frame * @param velocity velocity of spacecraft in reference frame * @param rotation orientation (attitude) of the spacecraft with respect to reference frame * @param mass spacecraft mass * @return acceleration with derivatives * @exception OrekitException if the underlying force models cannot compute the acceleration */ public FieldVector3D<DerivativeStructure> accelerationDerivatives(final AbsoluteDate date, final Frame frame, final FieldVector3D<DerivativeStructure> position, final FieldVector3D<DerivativeStructure> velocity, final FieldRotation<DerivativeStructure> rotation, final DerivativeStructure mass) throws OrekitException { final int parameters = mass.getFreeParameters(); final int order = mass.getOrder(); // estimate the scalar velocity step, assuming energy conservation // and hence differentiating equation V = sqrt(mu (2/r - 1/a)) final Vector3D p0 = position.toVector3D(); final Vector3D v0 = velocity.toVector3D(); final double r2 = p0.getNormSq(); final double hVel = mu * hPos / (v0.getNorm() * r2); // estimate mass step, applying the same relative value as position final double hMass = mass.getValue() * hPos / FastMath.sqrt(r2); // compute nominal acceleration final AccelerationRetriever nominal = new AccelerationRetriever(); computeShiftedAcceleration(nominal, date, frame, p0, v0, rotation.toRotation(), mass.getValue()); final double[] a0 = nominal.getAcceleration().toArray(); // compute accelerations with shifted states final AccelerationRetriever shifted = new AccelerationRetriever(); // shift position by hPos alon x, y and z computeShiftedAcceleration(shifted, date, frame, shift(position, 0, hPos), v0, shift(rotation, 0, hPos), shift(mass, 0, hPos)); final double[] derPx = new Vector3D(1 / hPos, shifted.getAcceleration(), -1 / hPos, nominal.getAcceleration()).toArray(); computeShiftedAcceleration(shifted, date, frame, shift(position, 1, hPos), v0, shift(rotation, 1, hPos), shift(mass, 1, hPos)); final double[] derPy = new Vector3D(1 / hPos, shifted.getAcceleration(), -1 / hPos, nominal.getAcceleration()).toArray(); computeShiftedAcceleration(shifted, date, frame, shift(position, 2, hPos), v0, shift(rotation, 2, hPos), shift(mass, 2, hPos)); final double[] derPz = new Vector3D(1 / hPos, shifted.getAcceleration(), -1 / hPos, nominal.getAcceleration()).toArray(); // shift velocity by hVel alon x, y and z computeShiftedAcceleration(shifted, date, frame, p0, shift(velocity, 3, hVel), shift(rotation, 3, hVel), shift(mass, 3, hVel)); final double[] derVx = new Vector3D(1 / hVel, shifted.getAcceleration(), -1 / hVel, nominal.getAcceleration()).toArray(); computeShiftedAcceleration(shifted, date, frame, p0, shift(velocity, 4, hVel), shift(rotation, 4, hVel), shift(mass, 4, hVel)); final double[] derVy = new Vector3D(1 / hVel, shifted.getAcceleration(), -1 / hVel, nominal.getAcceleration()).toArray(); computeShiftedAcceleration(shifted, date, frame, p0, shift(velocity, 5, hVel), shift(rotation, 5, hVel), shift(mass, 5, hVel)); final double[] derVz = new Vector3D(1 / hVel, shifted.getAcceleration(), -1 / hVel, nominal.getAcceleration()).toArray(); final double[] derM; if (parameters < 7) { derM = null; } else { // shift mass by hMass computeShiftedAcceleration(shifted, date, frame, p0, v0, shift(rotation, 6, hMass), shift(mass, 6, hMass)); derM = new Vector3D(1 / hMass, shifted.getAcceleration(), -1 / hMass, nominal.getAcceleration()) .toArray(); } final double[] derivatives = new double[1 + parameters]; final DerivativeStructure[] accDer = new DerivativeStructure[3]; for (int i = 0; i < 3; ++i) { // first element is value of acceleration derivatives[0] = a0[i]; // next three elements are derivatives with respect to position derivatives[1] = derPx[i]; derivatives[2] = derPy[i]; derivatives[3] = derPz[i]; // next three elements are derivatives with respect to velocity derivatives[4] = derVx[i]; derivatives[5] = derVy[i]; derivatives[6] = derVz[i]; if (derM != null) { derivatives[7] = derM[i]; } accDer[i] = new DerivativeStructure(parameters, order, derivatives); } return new FieldVector3D<DerivativeStructure>(accDer); }