List of usage examples for org.apache.commons.math3.analysis.differentiation DerivativeStructure getPartialDerivative
public double getPartialDerivative(final int... orders) throws DimensionMismatchException, NumberIsTooLargeException
From source file:com.bwc.ora.models.Lrp.java
/** * Get the local maximums from a collection of Points. * * @param lrpSeries//from ww w . j av a 2 s . com * @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: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(); }/*w ww . j av a 2s. c o 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.propagation.events.ElevationExtremumDetector.java
/** Compute the value of the detection function. * <p>/*w w w.j a v a 2 s.co m*/ * The value is the spacecraft elevation first time derivative. * </p> * @param s the current state information: date, kinematics, attitude * @return spacecraft elevation first time derivative * @exception OrekitException if some specific error occurs */ public double g(final SpacecraftState s) throws OrekitException { // get position, velocity acceleration of spacecraft in topocentric frame final Transform inertToTopo = s.getFrame().getTransformTo(topo, s.getDate()); final TimeStampedPVCoordinates pvTopo = inertToTopo.transformPVCoordinates(s.getPVCoordinates()); // convert the coordinates to DerivativeStructure based vector // instead of having vector position, then vector velocity then vector acceleration // we get one vector and each coordinate is a DerivativeStructure containing // value, first time derivative (we don't need second time derivative here) final FieldVector3D<DerivativeStructure> pvDS = pvTopo.toDerivativeStructureVector(1); // compute elevation and its first time derivative final DerivativeStructure elevation = pvDS.getZ().divide(pvDS.getNorm()).asin(); // return elevation first time derivative return elevation.getPartialDerivative(1); }
From source file:org.orekit.propagation.semianalytical.dsst.forces.TesseralContribution.java
/** Compute the n-SUM for potential derivatives components. * @param date current date/*from w w w . j av a2 s . c o m*/ * @param j resonant index <i>j</i> * @param m resonant order <i>m</i> * @param s d'Alembert characteristic <i>s</i> * @param maxN maximum possible value for <i>n</i> index * @param roaPow powers of R/a up to degree <i>n</i> * @param ghMSJ G<sup>j</sup><sub>m,s</sub> and H<sup>j</sup><sub>m,s</sub> polynomials * @param gammaMNS Γ<sup>m</sup><sub>n,s</sub>() function * @return Components of U<sub>n</sub> derivatives for fixed j, m, s * @throws OrekitException if some error occurred */ private double[][] computeNSum(final AbsoluteDate date, final int j, final int m, final int s, final int maxN, final double[] roaPow, final GHmsjPolynomials ghMSJ, final GammaMnsFunction gammaMNS) throws OrekitException { //spherical harmonics final UnnormalizedSphericalHarmonics harmonics = provider.onDate(date); // Potential derivatives components double dUdaCos = 0.; double dUdaSin = 0.; double dUdhCos = 0.; double dUdhSin = 0.; double dUdkCos = 0.; double dUdkSin = 0.; double dUdlCos = 0.; double dUdlSin = 0.; double dUdAlCos = 0.; double dUdAlSin = 0.; double dUdBeCos = 0.; double dUdBeSin = 0.; double dUdGaCos = 0.; double dUdGaSin = 0.; // I^m final int Im = I > 0 ? 1 : (m % 2 == 0 ? 1 : -1); // jacobi v, w, indices from 2.7.1-(15) final int v = FastMath.abs(m - s); final int w = FastMath.abs(m + s); // Initialise lower degree nmin = (Max(2, m, |s|)) for summation over n final int nmin = FastMath.max(FastMath.max(2, m), FastMath.abs(s)); //Get the corresponding Hansen object final int sIndex = maxDegree + (j < 0 ? -s : s); final int jIndex = FastMath.abs(j); final HansenTesseralLinear hans = this.hansenObjects[sIndex][jIndex]; // n-SUM from nmin to N for (int n = nmin; n <= maxN; n++) { // If (n - s) is odd, the contribution is null because of Vmns if ((n - s) % 2 == 0) { // Vmns coefficient final double fns = fact[n + FastMath.abs(s)]; final double vMNS = CoefficientsFactory.getVmns(m, n, s, fns, fact[n - m]); // Inclination function Gamma and derivative final double gaMNS = gammaMNS.getValue(m, n, s); final double dGaMNS = gammaMNS.getDerivative(m, n, s); // Hansen kernel value and derivative final double kJNS = hans.getValue(-n - 1, chi); final double dkJNS = hans.getDerivative(-n - 1, chi); // Gjms, Hjms polynomials and derivatives final double gMSJ = ghMSJ.getGmsj(m, s, j); final double hMSJ = ghMSJ.getHmsj(m, s, j); final double dGdh = ghMSJ.getdGmsdh(m, s, j); final double dGdk = ghMSJ.getdGmsdk(m, s, j); final double dGdA = ghMSJ.getdGmsdAlpha(m, s, j); final double dGdB = ghMSJ.getdGmsdBeta(m, s, j); final double dHdh = ghMSJ.getdHmsdh(m, s, j); final double dHdk = ghMSJ.getdHmsdk(m, s, j); final double dHdA = ghMSJ.getdHmsdAlpha(m, s, j); final double dHdB = ghMSJ.getdHmsdBeta(m, s, j); // Jacobi l-index from 2.7.1-(15) final int l = FastMath.min(n - m, n - FastMath.abs(s)); // Jacobi polynomial and derivative final DerivativeStructure jacobi = JacobiPolynomials.getValue(l, v, w, new DerivativeStructure(1, 1, 0, gamma)); // Geopotential coefficients final double cnm = harmonics.getUnnormalizedCnm(n, m); final double snm = harmonics.getUnnormalizedSnm(n, m); // Common factors from expansion of equations 3.3-4 final double cf_0 = roaPow[n] * Im * vMNS; final double cf_1 = cf_0 * gaMNS * jacobi.getValue(); final double cf_2 = cf_1 * kJNS; final double gcPhs = gMSJ * cnm + hMSJ * snm; final double gsMhc = gMSJ * snm - hMSJ * cnm; final double dKgcPhsx2 = 2. * dkJNS * gcPhs; final double dKgsMhcx2 = 2. * dkJNS * gsMhc; final double dUdaCoef = (n + 1) * cf_2; final double dUdlCoef = j * cf_2; final double dUdGaCoef = cf_0 * kJNS * (jacobi.getValue() * dGaMNS + gaMNS * jacobi.getPartialDerivative(1)); // dU / da components dUdaCos += dUdaCoef * gcPhs; dUdaSin += dUdaCoef * gsMhc; // dU / dh components dUdhCos += cf_1 * (kJNS * (cnm * dGdh + snm * dHdh) + h * dKgcPhsx2); dUdhSin += cf_1 * (kJNS * (snm * dGdh - cnm * dHdh) + h * dKgsMhcx2); // dU / dk components dUdkCos += cf_1 * (kJNS * (cnm * dGdk + snm * dHdk) + k * dKgcPhsx2); dUdkSin += cf_1 * (kJNS * (snm * dGdk - cnm * dHdk) + k * dKgsMhcx2); // dU / dLambda components dUdlCos += dUdlCoef * gsMhc; dUdlSin += -dUdlCoef * gcPhs; // dU / alpha components dUdAlCos += cf_2 * (dGdA * cnm + dHdA * snm); dUdAlSin += cf_2 * (dGdA * snm - dHdA * cnm); // dU / dBeta components dUdBeCos += cf_2 * (dGdB * cnm + dHdB * snm); dUdBeSin += cf_2 * (dGdB * snm - dHdB * cnm); // dU / dGamma components dUdGaCos += dUdGaCoef * gcPhs; dUdGaSin += dUdGaCoef * gsMhc; } } return new double[][] { { dUdaCos, dUdaSin }, { dUdhCos, dUdhSin }, { dUdkCos, dUdkSin }, { dUdlCos, dUdlSin }, { dUdAlCos, dUdAlSin }, { dUdBeCos, dUdBeSin }, { dUdGaCos, dUdGaSin } }; }
From source file:org.orekit.propagation.semianalytical.dsst.utilities.hansen.HansenTesseralLinear.java
/** * Compute the values for the first four coefficients and their derivatives by means of series. * * @param e2 e/*from ww w . j a v a 2 s . c om*/ * @param chi Χ * @param chi2 Χ */ public void computeInitValues(final double e2, final double chi, final double chi2) { // compute the values for n, n+1, n+2 and n+3 by series // See Danielson 2.7.3-(10) //Ensure that only the needed terms are computed final int maxRoots = FastMath.min(4, N0 - Nmin + 4); for (int i = 0; i < maxRoots; i++) { final DerivativeStructure hansenKernel = hansenInit[i].getValue(e2, chi, chi2); this.hansenRoot[0][i] = hansenKernel.getValue(); this.hansenDerivRoot[0][i] = hansenKernel.getPartialDerivative(1); } for (int i = 1; i < numSlices; i++) { for (int k = 0; k < 4; k++) { final PolynomialFunction[] mv = mpvec[N0 - (i * SLICE) - k + 3 + offset]; final PolynomialFunction[] sv = mpvecDeriv[N0 - (i * SLICE) - k + 3 + offset]; hansenDerivRoot[i][k] = mv[3].value(chi) * hansenDerivRoot[i - 1][3] + mv[2].value(chi) * hansenDerivRoot[i - 1][2] + mv[1].value(chi) * hansenDerivRoot[i - 1][1] + mv[0].value(chi) * hansenDerivRoot[i - 1][0] + sv[3].value(chi) * hansenRoot[i - 1][3] + sv[2].value(chi) * hansenRoot[i - 1][2] + sv[1].value(chi) * hansenRoot[i - 1][1] + sv[0].value(chi) * hansenRoot[i - 1][0]; hansenRoot[i][k] = mv[3].value(chi) * hansenRoot[i - 1][3] + mv[2].value(chi) * hansenRoot[i - 1][2] + mv[1].value(chi) * hansenRoot[i - 1][1] + mv[0].value(chi) * hansenRoot[i - 1][0]; } } }
From source file:org.orekit.utils.FieldPVCoordinatesTest.java
@Test public void testNormalize() { RandomGenerator generator = new Well19937a(0x7ede9376e4e1ab5al); FiniteDifferencesDifferentiator differentiator = new FiniteDifferencesDifferentiator(5, 1.0e-3); for (int i = 0; i < 200; ++i) { final FieldPVCoordinates<DerivativeStructure> pv = randomPVCoordinates(generator, 1e6, 1e3, 1.0); DerivativeStructure x = differentiator.differentiate(new UnivariateFunction() { public double value(double t) { return pv.shiftedBy(t).getPosition().normalize().getX().getValue(); }/*w w w. j a v a 2 s. c o m*/ }).value(new DerivativeStructure(1, 2, 0, 0.0)); DerivativeStructure y = differentiator.differentiate(new UnivariateFunction() { public double value(double t) { return pv.shiftedBy(t).getPosition().normalize().getY().getValue(); } }).value(new DerivativeStructure(1, 2, 0, 0.0)); DerivativeStructure z = differentiator.differentiate(new UnivariateFunction() { public double value(double t) { return pv.shiftedBy(t).getPosition().normalize().getZ().getValue(); } }).value(new DerivativeStructure(1, 2, 0, 0.0)); FieldPVCoordinates<DerivativeStructure> normalized = pv.normalize(); Assert.assertEquals(x.getValue(), normalized.getPosition().getX().getValue(), 1.0e-16); Assert.assertEquals(y.getValue(), normalized.getPosition().getY().getValue(), 1.0e-16); Assert.assertEquals(z.getValue(), normalized.getPosition().getZ().getValue(), 1.0e-16); Assert.assertEquals(x.getPartialDerivative(1), normalized.getVelocity().getX().getValue(), 3.0e-13); Assert.assertEquals(y.getPartialDerivative(1), normalized.getVelocity().getY().getValue(), 3.0e-13); Assert.assertEquals(z.getPartialDerivative(1), normalized.getVelocity().getZ().getValue(), 3.0e-13); Assert.assertEquals(x.getPartialDerivative(2), normalized.getAcceleration().getX().getValue(), 6.0e-10); Assert.assertEquals(y.getPartialDerivative(2), normalized.getAcceleration().getY().getValue(), 6.0e-10); Assert.assertEquals(z.getPartialDerivative(2), normalized.getAcceleration().getZ().getValue(), 6.0e-10); } }
From source file:org.orekit.utils.IERSConventionsTest.java
@Test public void testGMST2000vs82() throws OrekitException { final TimeFunction<DerivativeStructure> gmst82 = IERSConventions.IERS_1996 .getGMSTFunction(TimeScalesFactory.getUT1(IERSConventions.IERS_1996, true)); final TimeFunction<DerivativeStructure> gmst00 = IERSConventions.IERS_2003 .getGMSTFunction(TimeScalesFactory.getUT1(IERSConventions.IERS_2003, true)); for (double dt = 0; dt < 10 * Constants.JULIAN_YEAR; dt += 10 * Constants.JULIAN_DAY) { AbsoluteDate date = new AbsoluteDate(AbsoluteDate.J2000_EPOCH, dt); DerivativeStructure delta = gmst00.value(date).subtract(gmst82.value(date)); // GMST82 and GMST2000 are similar to about 15 milli-arcseconds between 2000 and 2010 Assert.assertEquals(0.0, MathUtils.normalizeAngle(delta.getValue(), 0.0), 7.1e-8); Assert.assertEquals(0.0, delta.getPartialDerivative(1), 1.0e-15); }//from w w w .j ava2s . co m }
From source file:org.orekit.utils.IERSConventionsTest.java
@Test public void testGMST2000vs2006() throws OrekitException { final UT1Scale ut1 = TimeScalesFactory.getUT1(IERSConventions.IERS_2010, true); final TimeFunction<DerivativeStructure> gmst00 = IERSConventions.IERS_2003.getGMSTFunction(ut1); final TimeFunction<DerivativeStructure> gmst06 = IERSConventions.IERS_2010.getGMSTFunction(ut1); for (double dt = 0; dt < 10 * Constants.JULIAN_YEAR; dt += 10 * Constants.JULIAN_DAY) { AbsoluteDate date = new AbsoluteDate(AbsoluteDate.J2000_EPOCH, dt); DerivativeStructure delta = gmst06.value(date).subtract(gmst00.value(date)); // GMST2006 and GMST2000 are similar to about 0.15 milli-arcseconds between 2000 and 2010 Assert.assertEquals(0.0, MathUtils.normalizeAngle(delta.getValue(), 0.0), 7e-10); Assert.assertEquals(0.0, delta.getPartialDerivative(1), 3.0e-18); }//w ww . ja v a 2s. c o m }
From source file:org.orekit.utils.IERSConventionsTest.java
private void checkDerivative(final TimeFunction<DerivativeStructure> function, final AbsoluteDate date, final double span, final double sampleStep, final double h, final double tolerance) { UnivariateDifferentiableFunction differentiated = new FiniteDifferencesDifferentiator(4, h) .differentiate(new UnivariateFunction() { @Override//from w w w. ja v a 2 s . c om public double value(final double dt) { return function.value(date.shiftedBy(dt)).getValue(); } }); for (double dt = 0; dt < span; dt += sampleStep) { DerivativeStructure yRef = differentiated.value(new DerivativeStructure(1, 1, 0, dt)); DerivativeStructure y = function.value(date.shiftedBy(dt)); Assert.assertEquals(yRef.getPartialDerivative(1), y.getPartialDerivative(1), tolerance); } }
From source file:org.orekit.utils.PVCoordinatesTest.java
@Test public void testNormalize() { RandomGenerator generator = new Well19937a(0xb2011ffd25412067l); FiniteDifferencesDifferentiator differentiator = new FiniteDifferencesDifferentiator(5, 1.0e-3); for (int i = 0; i < 200; ++i) { final PVCoordinates pv = randomPVCoordinates(generator, 1e6, 1e3, 1.0); DerivativeStructure x = differentiator.differentiate(new UnivariateFunction() { public double value(double t) { return pv.shiftedBy(t).getPosition().normalize().getX(); }// w w w . j a va 2 s.co m }).value(new DerivativeStructure(1, 2, 0, 0.0)); DerivativeStructure y = differentiator.differentiate(new UnivariateFunction() { public double value(double t) { return pv.shiftedBy(t).getPosition().normalize().getY(); } }).value(new DerivativeStructure(1, 2, 0, 0.0)); DerivativeStructure z = differentiator.differentiate(new UnivariateFunction() { public double value(double t) { return pv.shiftedBy(t).getPosition().normalize().getZ(); } }).value(new DerivativeStructure(1, 2, 0, 0.0)); PVCoordinates normalized = pv.normalize(); Assert.assertEquals(x.getValue(), normalized.getPosition().getX(), 1.0e-16); Assert.assertEquals(y.getValue(), normalized.getPosition().getY(), 1.0e-16); Assert.assertEquals(z.getValue(), normalized.getPosition().getZ(), 1.0e-16); Assert.assertEquals(x.getPartialDerivative(1), normalized.getVelocity().getX(), 3.0e-13); Assert.assertEquals(y.getPartialDerivative(1), normalized.getVelocity().getY(), 3.0e-13); Assert.assertEquals(z.getPartialDerivative(1), normalized.getVelocity().getZ(), 3.0e-13); Assert.assertEquals(x.getPartialDerivative(2), normalized.getAcceleration().getX(), 6.0e-10); Assert.assertEquals(y.getPartialDerivative(2), normalized.getAcceleration().getY(), 6.0e-10); Assert.assertEquals(z.getPartialDerivative(2), normalized.getAcceleration().getZ(), 6.0e-10); } }