List of usage examples for org.apache.commons.math3.util FastMath cbrt
public static double cbrt(double x)
From source file:de.tuberlin.uebb.jbop.example.DSCompilerOnlyCompose.java
@Override public void rootN(final double[] operand, final int n, final double[] result) { // create the function value and derivatives // [x^(1/n), (1/n)x^((1/n)-1), (1-n)/n^2x^((1/n)-2), ... ] final double[] function = new double[1 + order]; double xk;// www . j a v a 2s.c o m if (n == 2) { function[0] = FastMath.sqrt(operand[0]); xk = 0.5 / function[0]; } else if (n == 3) { function[0] = FastMath.cbrt(operand[0]); xk = 1.0 / (3.0 * function[0] * function[0]); } else { function[0] = FastMath.pow(operand[0], 1.0 / n); xk = 1.0 / (n * FastMath.pow(function[0], n - 1)); } final double nReciprocal = 1.0 / n; final double xReciprocal = 1.0 / operand[0]; for (int i = 1; i <= order; ++i) { function[i] = xk; xk *= xReciprocal * (nReciprocal - i); } // apply function composition compose(operand, function, result); }
From source file:de.tuberlin.uebb.jbop.example.DSCompiler.java
@Override @Optimizable// www .ja va 2 s . c o m @StrictLoops public void rootN(final double[] operand, final int n, final double[] result) { // create the function value and derivatives // [x^(1/n), (1/n)x^((1/n)-1), (1-n)/n^2x^((1/n)-2), ... ] final double[] function = new double[1 + order]; double xk; if (n == 2) { function[0] = FastMath.sqrt(operand[0]); xk = 0.5 / function[0]; } else if (n == 3) { function[0] = FastMath.cbrt(operand[0]); xk = 1.0 / (3.0 * function[0] * function[0]); } else { function[0] = FastMath.pow(operand[0], 1.0 / n); xk = 1.0 / (n * FastMath.pow(function[0], n - 1)); } final double nReciprocal = 1.0 / n; final double xReciprocal = 1.0 / operand[0]; for (int i = 1; i <= order; ++i) { function[i] = xk; xk *= xReciprocal * (nReciprocal - i); } // apply function composition compose(operand, function, result); }
From source file:fr.cs.examples.bodies.Phasing.java
/** Guess an initial orbit from theoretical model. * @param date orbit date//from w w w . j a va 2s. c o m * @param frame frame to use for defining orbit * @param nbOrbits number of orbits in the phasing cycle * @param nbDays number of days in the phasing cycle * @param latitude reference latitude for Sun synchronous orbit * @param ascending if true, crossing latitude is from South to North * @param mst desired mean solar time at reference latitude crossing * @return an initial guess of Earth phased, Sun synchronous orbit * @exception OrekitException if mean solar time cannot be computed */ private CircularOrbit guessOrbit(AbsoluteDate date, Frame frame, int nbOrbits, int nbDays, double latitude, boolean ascending, double mst) throws OrekitException { double mu = gravityField.getMu(); NormalizedSphericalHarmonicsProvider.NormalizedSphericalHarmonics harmonics = gravityField.onDate(date); // initial semi major axis guess based on Keplerian period double period0 = (nbDays * Constants.JULIAN_DAY) / nbOrbits; double n0 = 2 * FastMath.PI / period0; double a0 = FastMath.cbrt(mu / (n0 * n0)); // initial inclination guess based on ascending node drift due to J2 double[][] unnormalization = GravityFieldFactory.getUnnormalizationFactors(3, 0); double j2 = -unnormalization[2][0] * harmonics.getNormalizedCnm(2, 0); double j3 = -unnormalization[3][0] * harmonics.getNormalizedCnm(3, 0); double raanRate = 2 * FastMath.PI / Constants.JULIAN_YEAR; double ae = gravityField.getAe(); double i0 = FastMath.acos(-raanRate * a0 * a0 / (1.5 * ae * ae * j2 * n0)); // initial eccentricity guess based on J2 and J3 double ex0 = 0; double ey0 = -j3 * ae * FastMath.sin(i0) / (2 * a0 * j2); // initial ascending node guess based on mean solar time double alpha0 = FastMath.asin(FastMath.sin(latitude) / FastMath.sin(i0)); if (!ascending) { alpha0 = FastMath.PI - alpha0; } double h = meanSolarTime( new CircularOrbit(a0, ex0, ey0, i0, 0.0, alpha0, PositionAngle.TRUE, frame, date, mu)); double raan0 = FastMath.PI * (mst - h) / 12.0; return new CircularOrbit(a0, ex0, ey0, i0, raan0, alpha0, PositionAngle.TRUE, frame, date, mu); }
From source file:org.esa.beam.util.math.FastMathPerformance.java
public void testCbrt() { System.gc();/*from ww w.jav a 2s . c om*/ double x = 0; long time = System.nanoTime(); for (int i = 0; i < RUNS; i++) x += StrictMath.cbrt(i * F1); long strictTime = System.nanoTime() - time; System.gc(); double y = 0; time = System.nanoTime(); for (int i = 0; i < RUNS; i++) y += FastMath.cbrt(i * F1); long fastTime = System.nanoTime() - time; System.gc(); double z = 0; time = System.nanoTime(); for (int i = 0; i < RUNS; i++) z += Math.cbrt(i * F1); long mathTime = System.nanoTime() - time; report("cbrt", x + y + z, strictTime, fastTime, mathTime); }
From source file:org.orekit.files.ccsds.OMMFile.java
/** Generate a {@link KeplerianOrbit} based on the OMM mean keplerian elements. * If the reference frame is not pseudo-inertial, an exception is raised. * @return the {@link KeplerianOrbit} generated from the OMM information * @exception OrekitException if the reference frame is not pseudo-inertial or if the central body * gravitational coefficient cannot be retrieved from the OMM *//* w ww . ja v a 2 s . c o m*/ public KeplerianOrbit generateKeplerianOrbit() throws OrekitException { setMuUsed(); final double a; if (Double.isNaN(getA())) { a = FastMath.cbrt(getMuUsed() / (meanMotion * meanMotion)); } else { a = getA(); } return new KeplerianOrbit(a, getE(), getI(), getPa(), getRaan(), getAnomaly(), PositionAngle.MEAN, metaData.getFrame(), getEpoch(), getMuUsed()); }
From source file:org.orekit.orbits.KeplerianOrbit.java
/** Computes the elliptic eccentric anomaly from the mean anomaly. * <p>// ww w . j av a 2 s. c o m * The algorithm used here for solving Kepler equation has been published * in: "Procedures for solving Kepler's Equation", A. W. Odell and * R. H. Gooding, Celestial Mechanics 38 (1986) 307-334 * </p> * @param M mean anomaly (rad) * @return v the true anomaly */ private double meanToEllipticEccentric(final double M) { // reduce M to [-PI PI) interval final double reducedM = MathUtils.normalizeAngle(M, 0.0); // compute start value according to A. W. Odell and R. H. Gooding S12 starter double E; if (FastMath.abs(reducedM) < 1.0 / 6.0) { E = reducedM + e * (FastMath.cbrt(6 * reducedM) - reducedM); } else { if (reducedM < 0) { final double w = FastMath.PI + reducedM; E = reducedM + e * (A * w / (B - w) - FastMath.PI - reducedM); } else { final double w = FastMath.PI - reducedM; E = reducedM + e * (FastMath.PI - A * w / (B - w) - reducedM); } } final double e1 = 1 - e; final boolean noCancellationRisk = (e1 + E * E / 6) >= 0.1; // perform two iterations, each consisting of one Halley step and one Newton-Raphson step for (int j = 0; j < 2; ++j) { double f; double fd; final double fdd = e * FastMath.sin(E); final double fddd = e * FastMath.cos(E); if (noCancellationRisk) { f = (E - fdd) - reducedM; fd = 1 - fddd; } else { f = eMeSinE(E) - reducedM; final double s = FastMath.sin(0.5 * E); fd = e1 + 2 * e * s * s; } final double dee = f * fd / (0.5 * f * fdd - fd * fd); // update eccentric anomaly, using expressions that limit underflow problems final double w = fd + 0.5 * dee * (fdd + dee * fddd / 3); fd += dee * (fdd + 0.5 * dee * fddd); E -= (f - dee * (fd - w)) / fd; } // expand the result back to original range E += M - reducedM; return E; }
From source file:org.orekit.propagation.semianalytical.dsst.forces.DSSTSolarRadiationPressure.java
/** * Compute the real roots of a cubic equation. * * <pre>/*from w w w . ja v a 2s . com*/ * a[0] * x + a[1] * x + a[2] * x + a[3] = 0. * </pre> * * @param a the 4 coefficients * @param y the real roots sorted in descending order * @return the number of real roots */ private int realCubicRoots(final double[] a, final double[] y) { if (Precision.equals(a[0], 0.)) { // Treat the degenerate cubic as quadratic final double[] aa = new double[a.length - 1]; System.arraycopy(a, 1, aa, 0, aa.length); return realQuadraticRoots(aa, y); } // Transform coefficients final double b = -a[1] / (3. * a[0]); final double c = a[2] / a[0]; final double d = a[3] / a[0]; final double b2 = b * b; final double p = b2 - c / 3.; final double q = b * (b2 - c * 0.5) - d * 0.5; // Compute discriminant final double disc = p * p * p - q * q; if (disc < 0.) { // One real root final double alpha = q + FastMath.copySign(FastMath.sqrt(-disc), q); final double cbrtAl = FastMath.cbrt(alpha); final double cbrtBe = p / cbrtAl; if (p < 0.) { y[0] = b + 2. * q / (cbrtAl * cbrtAl + cbrtBe * cbrtBe - p); } else if (p > 0.) { y[0] = b + cbrtAl + cbrtBe; } else { y[0] = b + cbrtAl; } return 1; } else if (disc > 0.) { // Three distinct real roots final double phi = FastMath.atan2(FastMath.sqrt(disc), q) / 3.; final double sqP = 2.0 * FastMath.sqrt(p); y[0] = b + sqP * FastMath.cos(phi); y[1] = b - sqP * FastMath.cos(FastMath.PI / 3. + phi); y[2] = b - sqP * FastMath.cos(FastMath.PI / 3. - phi); return 3; } else { // One distinct and two equals real roots final double cbrtQ = FastMath.cbrt(q); final double root1 = b + 2. * cbrtQ; final double root2 = b - cbrtQ; if (q < 0.) { y[0] = root2; y[1] = root2; y[2] = root1; } else { y[0] = root1; y[1] = root2; y[2] = root2; } return 3; } }