List of usage examples for org.apache.commons.math3.util FastMath acos
public static double acos(double x)
From source file:net.sf.dsp4j.octave.packages.signal_1_0_11.Cheb.java
public static double[] cheb(int n, double[] x) { if (n <= 0) { throw new IllegalArgumentException("cheb: n has to be a positive integer"); }// w w w .ja v a2s . c o m if (x.length == 0) { return new double[0]; } //# avoid resizing latencies double[] T = new double[x.length]; for (int i = 0; i < x.length; i++) { if (Math.abs(x[i]) > 1) { T[i] = FastMath.cos(n * FastMath.acos(x[i])); } else { T[i] = FastMath.cosh(n * FastMath.acosh(x[i])); } } return T; }
From source file:com.clust4j.kernel.CircularKernel.java
@Override public double getPartialSimilarity(double[] a, double[] b) { final double lp = toHilbertPSpace(a, b); // Per corner case condition if (lp >= getSigma()) return 0.0; final double twoOverPi = (2d / FastMath.PI); final double lpOverSig = lp / getSigma(); /* Front segment */ final double front = twoOverPi * FastMath.acos(-lpOverSig); /* Back segment */ final double first = twoOverPi * lpOverSig; final double second = FastMath.sqrt(1.0 - FastMath.pow(lpOverSig, 2)); final double back = first * second; final double answer = front - back; return Double.isNaN(answer) ? Double.NEGATIVE_INFINITY : answer; }
From source file:net.nicoulaj.benchmarks.math.DoubleAcos.java
@GenerateMicroBenchmark public void commonsmath(BlackHole hole) { for (int i = 0; i < data.length - 1; i++) hole.consume(FastMath.acos(data[i])); }
From source file:cc.redberry.core.number.ComplexUtils.java
/** * Numeric arccosine form complex number. * * @param complex argument//from w ww . j a va2 s . c o m * @return arccosine */ public static Complex arccos(Complex complex) { if (complex.isReal()) { double x = complex.getReal().doubleValue(); if (x <= 1.0 && x >= -1) return new Complex(FastMath.acos(complex.getReal().doubleValue())); } return new Complex(new org.apache.commons.math3.complex.Complex(complex.getReal().doubleValue(), complex.getImaginary().doubleValue()).acos()); }
From source file:lambertmrev.Lambert.java
/** Constructs and solves a Lambert problem. * * \param[in] R1 first cartesian position * \param[in] R2 second cartesian position * \param[in] tof time of flight/*from w w w .j ava 2 s . co m*/ * \param[in] mu gravity parameter * \param[in] cw when 1 a retrograde orbit is assumed * \param[in] multi_revs maximum number of multirevolutions to compute */ public void lambert_problem(Vector3D r1, Vector3D r2, double tof, double mu, Boolean cw, int multi_revs) { // sanity checks if (tof <= 0) { System.out.println("ToF is negative! \n"); } if (mu <= 0) { System.out.println("mu is below zero"); } // 1 - getting lambda and T double m_c = FastMath.sqrt((r2.getX() - r1.getX()) * (r2.getX() - r1.getX()) + (r2.getY() - r1.getY()) * (r2.getY() - r1.getY()) + (r2.getZ() - r1.getZ()) * (r2.getZ() - r1.getZ())); double R1 = r1.getNorm(); double R2 = r2.getNorm(); double m_s = (m_c + R1 + R2) / 2.0; Vector3D ir1 = r1.normalize(); Vector3D ir2 = r2.normalize(); Vector3D ih = Vector3D.crossProduct(ir1, ir2); ih = ih.normalize(); if (ih.getZ() == 0) { System.out.println("angular momentum vector has no z component \n"); } double lambda2 = 1.0 - m_c / m_s; double m_lambda = FastMath.sqrt(lambda2); Vector3D it1 = new Vector3D(0.0, 0.0, 0.0); Vector3D it2 = new Vector3D(0.0, 0.0, 0.0); if (ih.getZ() < 0.0) { // Transfer angle is larger than 180 degrees as seen from abive the z axis m_lambda = -m_lambda; it1 = Vector3D.crossProduct(ir1, ih); it2 = Vector3D.crossProduct(ir2, ih); } else { it1 = Vector3D.crossProduct(ih, ir1); it2 = Vector3D.crossProduct(ih, ir2); } it1.normalize(); it2.normalize(); if (cw) { // Retrograde motion m_lambda = -m_lambda; it1.negate(); it2.negate(); } double lambda3 = m_lambda * lambda2; double T = FastMath.sqrt(2.0 * mu / m_s / m_s / m_s) * tof; // 2 - We now hava lambda, T and we will find all x // 2.1 - let us first detect the maximum number of revolutions for which there exists a solution int m_Nmax = FastMath.toIntExact(FastMath.round(T / FastMath.PI)); double T00 = FastMath.acos(m_lambda) + m_lambda * FastMath.sqrt(1.0 - lambda2); double T0 = (T00 + m_Nmax * FastMath.PI); double T1 = 2.0 / 3.0 * (1.0 - lambda3); double DT = 0.0; double DDT = 0.0; double DDDT = 0.0; if (m_Nmax > 0) { if (T < T0) { // We use Halley iterations to find xM and TM int it = 0; double err = 1.0; double T_min = T0; double x_old = 0.0, x_new = 0.0; while (true) { ArrayRealVector deriv = dTdx(x_old, T_min, m_lambda); DT = deriv.getEntry(0); DDT = deriv.getEntry(1); DDDT = deriv.getEntry(2); if (DT != 0.0) { x_new = x_old - DT * DDT / (DDT * DDT - DT * DDDT / 2.0); } err = FastMath.abs(x_old - x_new); if ((err < 1e-13) || (it > 12)) { break; } tof = x2tof(x_new, m_Nmax, m_lambda); x_old = x_new; it++; } if (T_min > T) { m_Nmax -= 1; } } } // We exit this if clause with Mmax being the maximum number of revolutions // for which there exists a solution. We crop it to multi_revs m_Nmax = FastMath.min(multi_revs, m_Nmax); // 2.2 We now allocate the memory for the output variables m_v1 = MatrixUtils.createRealMatrix(m_Nmax * 2 + 1, 3); RealMatrix m_v2 = MatrixUtils.createRealMatrix(m_Nmax * 2 + 1, 3); RealMatrix m_iters = MatrixUtils.createRealMatrix(m_Nmax * 2 + 1, 3); //RealMatrix m_x = MatrixUtils.createRealMatrix(m_Nmax*2+1, 3); ArrayRealVector m_x = new ArrayRealVector(m_Nmax * 2 + 1); // 3 - We may now find all solution in x,y // 3.1 0 rev solution // 3.1.1 initial guess if (T >= T00) { m_x.setEntry(0, -(T - T00) / (T - T00 + 4)); } else if (T <= T1) { m_x.setEntry(0, T1 * (T1 - T) / (2.0 / 5.0 * (1 - lambda2 * lambda3) * T) + 1); } else { m_x.setEntry(0, FastMath.pow((T / T00), 0.69314718055994529 / FastMath.log(T1 / T00)) - 1.0); } // 3.1.2 Householder iterations //m_iters.setEntry(0, 0, housOutTmp.getEntry(0)); m_x.setEntry(0, householder(T, m_x.getEntry(0), 0, 1e-5, 15, m_lambda)); // 3.2 multi rev solutions double tmp; double x0; for (int i = 1; i < m_Nmax + 1; i++) { // 3.2.1 left householder iterations tmp = FastMath.pow((i * FastMath.PI + FastMath.PI) / (8.0 * T), 2.0 / 3.0); m_x.setEntry(2 * i - 1, (tmp - 1) / (tmp + 1)); x0 = householder(T, m_x.getEntry(2 * i - 1), i, 1e-8, 15, m_lambda); m_x.setEntry(2 * i - 1, x0); //m_iters.setEntry(2*i-1, 0, housOutTmp.getEntry(0)); //3.2.1 right Householder iterations tmp = FastMath.pow((8.0 * T) / (i * FastMath.PI), 2.0 / 3.0); m_x.setEntry(2 * i, (tmp - 1) / (tmp + 1)); x0 = householder(T, m_x.getEntry(2 * i), i, 1e-8, 15, m_lambda); m_x.setEntry(2 * i, x0); //m_iters.setEntry(2*i, 0, housOutTmp.getEntry(0)); } // 4 - For each found x value we recontruct the terminal velocities double gamma = FastMath.sqrt(mu * m_s / 2.0); double rho = (R1 - R2) / m_c; double sigma = FastMath.sqrt(1 - rho * rho); double vr1, vt1, vr2, vt2, y; ArrayRealVector ir1_vec = new ArrayRealVector(3); ArrayRealVector ir2_vec = new ArrayRealVector(3); ArrayRealVector it1_vec = new ArrayRealVector(3); ArrayRealVector it2_vec = new ArrayRealVector(3); // set Vector3D values to a mutable type ir1_vec.setEntry(0, ir1.getX()); ir1_vec.setEntry(1, ir1.getY()); ir1_vec.setEntry(2, ir1.getZ()); ir2_vec.setEntry(0, ir2.getX()); ir2_vec.setEntry(1, ir2.getY()); ir2_vec.setEntry(2, ir2.getZ()); it1_vec.setEntry(0, it1.getX()); it1_vec.setEntry(1, it1.getY()); it1_vec.setEntry(2, it1.getZ()); it2_vec.setEntry(0, it2.getX()); it2_vec.setEntry(1, it2.getY()); it2_vec.setEntry(2, it2.getZ()); for (int i = 0; i < m_x.getDimension(); i++) { y = FastMath.sqrt(1.0 - lambda2 + lambda2 * m_x.getEntry(i) * m_x.getEntry(i)); vr1 = gamma * ((m_lambda * y - m_x.getEntry(i)) - rho * (m_lambda * y + m_x.getEntry(i))) / R1; vr2 = -gamma * ((m_lambda * y - m_x.getEntry(i)) + rho * (m_lambda * y + m_x.getEntry(i))) / R2; double vt = gamma * sigma * (y + m_lambda * m_x.getEntry(i)); vt1 = vt / R1; vt2 = vt / R2; for (int j = 0; j < 3; ++j) m_v1.setEntry(i, j, vr1 * ir1_vec.getEntry(j) + vt1 * it1_vec.getEntry(j)); for (int j = 0; j < 3; ++j) m_v2.setEntry(i, j, vr2 * ir2_vec.getEntry(j) + vt2 * it2_vec.getEntry(j)); } }
From source file:com.bc.jexp.impl.DefaultNamespace.java
private void registerDefaultFunctions() { registerFunction(new AbstractFunction.D("sin", 1) { public double evalD(final EvalEnv env, final Term[] args) { return FastMath.sin(args[0].evalD(env)); }/*from www . j a v a 2 s .c om*/ }); registerFunction(new AbstractFunction.D("cos", 1) { public double evalD(final EvalEnv env, final Term[] args) { return FastMath.cos(args[0].evalD(env)); } }); registerFunction(new AbstractFunction.D("tan", 1) { public double evalD(final EvalEnv env, final Term[] args) { return FastMath.tan(args[0].evalD(env)); } }); registerFunction(new AbstractFunction.D("asin", 1) { public double evalD(final EvalEnv env, final Term[] args) { return FastMath.asin(args[0].evalD(env)); } }); registerFunction(new AbstractFunction.D("acos", 1) { public double evalD(final EvalEnv env, final Term[] args) { return FastMath.acos(args[0].evalD(env)); } }); registerFunction(new AbstractFunction.D("atan", 1) { public double evalD(final EvalEnv env, final Term[] args) { return FastMath.atan(args[0].evalD(env)); } }); registerFunction(new AbstractFunction.D("atan2", 2) { public double evalD(final EvalEnv env, final Term[] args) { return Math.atan2(args[0].evalD(env), args[1].evalD(env)); } }); registerFunction(new AbstractFunction.D("log", 1) { public double evalD(final EvalEnv env, final Term[] args) { return Math.log(args[0].evalD(env)); } }); registerFunction(new AbstractFunction.D("log10", 1) { public double evalD(EvalEnv env, Term[] args) throws EvalException { return Math.log10(args[0].evalD(env)); } }); registerFunction(new AbstractFunction.D("exp", 1) { public double evalD(final EvalEnv env, final Term[] args) { return FastMath.exp(args[0].evalD(env)); } }); registerFunction(new AbstractFunction.D("exp10", 1) { public double evalD(EvalEnv env, Term[] args) throws EvalException { return FastMath.pow(10.0, args[0].evalD(env)); } }); registerFunction(new AbstractFunction.D("sqr", 1) { public double evalD(final EvalEnv env, final Term[] args) { double v = args[0].evalD(env); return v * v; } }); registerFunction(new AbstractFunction.D("sqrt", 1) { public double evalD(final EvalEnv env, final Term[] args) { return Math.sqrt(args[0].evalD(env)); } }); registerFunction(new AbstractFunction.D("pow", 2) { public double evalD(final EvalEnv env, final Term[] args) { return FastMath.pow(args[0].evalD(env), args[1].evalD(env)); } }); registerFunction(new AbstractFunction.I("min", 2) { public int evalI(final EvalEnv env, final Term[] args) { return Math.min(args[0].evalI(env), args[1].evalI(env)); } }); registerFunction(new AbstractFunction.D("min", 2) { public double evalD(final EvalEnv env, final Term[] args) { return Math.min(args[0].evalD(env), args[1].evalD(env)); } }); registerFunction(new AbstractFunction.I("max", 2) { public int evalI(final EvalEnv env, final Term[] args) { return Math.max(args[0].evalI(env), args[1].evalI(env)); } }); registerFunction(new AbstractFunction.D("max", 2) { public double evalD(final EvalEnv env, final Term[] args) { return Math.max(args[0].evalD(env), args[1].evalD(env)); } }); registerFunction(new AbstractFunction.D("floor", 1) { public double evalD(EvalEnv env, Term[] args) throws EvalException { return Math.floor(args[0].evalD(env)); } }); registerFunction(new AbstractFunction.D("round", 1) { public double evalD(EvalEnv env, Term[] args) throws EvalException { return Math.round(args[0].evalD(env)); } }); registerFunction(new AbstractFunction.D("ceil", 1) { public double evalD(EvalEnv env, Term[] args) throws EvalException { return Math.ceil(args[0].evalD(env)); } }); registerFunction(new AbstractFunction.D("rint", 1) { public double evalD(EvalEnv env, Term[] args) throws EvalException { return Math.rint(args[0].evalD(env)); } }); registerFunction(new AbstractFunction.I("sign", 1) { public int evalI(final EvalEnv env, final Term[] args) { return ExtMath.sign(args[0].evalI(env)); } }); registerFunction(new AbstractFunction.D("sign", 1) { public double evalD(final EvalEnv env, final Term[] args) { return ExtMath.sign(args[0].evalD(env)); } }); registerFunction(new AbstractFunction.I("abs", 1) { public int evalI(final EvalEnv env, final Term[] args) { return Math.abs(args[0].evalI(env)); } }); registerFunction(new AbstractFunction.D("abs", 1) { public double evalD(final EvalEnv env, final Term[] args) { return Math.abs(args[0].evalD(env)); } }); registerFunction(new AbstractFunction.D("deg", 1) { public double evalD(final EvalEnv env, final Term[] args) { return Math.toDegrees(args[0].evalD(env)); } }); registerFunction(new AbstractFunction.D("rad", 1) { public double evalD(final EvalEnv env, final Term[] args) { return Math.toRadians(args[0].evalD(env)); } }); registerFunction(new AbstractFunction.D("ampl", 2) { public double evalD(final EvalEnv env, final Term[] args) { final double a = args[0].evalD(env); final double b = args[1].evalD(env); return Math.sqrt(a * a + b * b); } }); registerFunction(new AbstractFunction.D("phase", 2) { public double evalD(final EvalEnv env, final Term[] args) { final double a = args[0].evalD(env); final double b = args[1].evalD(env); return Math.atan2(b, a); } }); registerFunction(new AbstractFunction.B("feq", 2) { public boolean evalB(EvalEnv env, Term[] args) throws EvalException { final double x1 = args[0].evalD(env); final double x2 = args[1].evalD(env); return ExtMath.feq(x1, x2, EPS); } }); registerFunction(new AbstractFunction.B("feq", 3) { public boolean evalB(EvalEnv env, Term[] args) throws EvalException { final double x1 = args[0].evalD(env); final double x2 = args[1].evalD(env); final double eps = args[2].evalD(env); return ExtMath.feq(x1, x2, eps); } }); registerFunction(new AbstractFunction.B("fneq", 2) { public boolean evalB(EvalEnv env, Term[] args) throws EvalException { final double x1 = args[0].evalD(env); final double x2 = args[1].evalD(env); return ExtMath.fneq(x1, x2, EPS); } }); registerFunction(new AbstractFunction.B("fneq", 3) { public boolean evalB(EvalEnv env, Term[] args) throws EvalException { final double x1 = args[0].evalD(env); final double x2 = args[1].evalD(env); final double eps = args[2].evalD(env); return ExtMath.fneq(x1, x2, eps); } }); registerFunction(new AbstractFunction.B("inf", 1) { public boolean evalB(EvalEnv env, Term[] args) throws EvalException { return Double.isInfinite(args[0].evalD(env)); } }); registerFunction(new AbstractFunction.B("nan", 1) { public boolean evalB(EvalEnv env, Term[] args) throws EvalException { return Double.isNaN(args[0].evalD(env)); } }); registerFunction(new AbstractFunction.D("distance", -1) { public double evalD(final EvalEnv env, final Term[] args) { double sqrSum = 0.0; final int n = args.length / 2; for (int i = 0; i < n; i++) { final double v = args[i + n].evalD(env) - args[i].evalD(env); sqrSum += v * v; } return Math.sqrt(sqrSum); } }); registerFunction(new AbstractFunction.D("distance_deriv", -1) { public double evalD(final EvalEnv env, final Term[] args) { double sqrSum = 0.0; final int n = args.length / 2; for (int i = 0; i < n - 1; i++) { final double v1 = args[i + 1].evalD(env) - args[i].evalD(env); final double v2 = args[i + n + 1].evalD(env) - args[i + n].evalD(env); sqrSum += (v1 - v2) * (v1 - v2); } return Math.sqrt(sqrSum); } }); registerFunction(new AbstractFunction.D("distance_integ", -1) { public double evalD(final EvalEnv env, final Term[] args) { double sqrSum = 0.0; double v1Sum = 0.0; double v2Sum = 0.0; final int n = args.length / 2; for (int i = 0; i < n; i++) { v1Sum += args[i].evalD(env); v2Sum += args[i + n].evalD(env); sqrSum += (v2Sum - v1Sum) * (v2Sum - v1Sum); } return Math.sqrt(sqrSum); } }); registerFunction(new AbstractFunction.B("inrange", -1) { public boolean evalB(final EvalEnv env, final Term[] args) { final int n1 = args.length / 3; final int n2 = n1 + args.length / 3; for (int i = 0; i < n1; i++) { final double v = args[i].evalD(env); final double v1 = args[i + n1].evalD(env); final double v2 = args[i + n2].evalD(env); if (v < v1 || v > v2) { return false; } } return true; } }); registerFunction(new AbstractFunction.B("inrange_deriv", -1) { public boolean evalB(final EvalEnv env, final Term[] args) { final int n1 = args.length / 3; final int n2 = 2 * n1; for (int i = 0; i < n1 - 1; i++) { final double v = args[i + 1].evalD(env) - args[i].evalD(env); final double v1 = args[i + n1 + 1].evalD(env) - args[i + n1].evalD(env); final double v2 = args[i + n2 + 1].evalD(env) - args[i + n2].evalD(env); if (v < v1 || v > v2) { return false; } } return true; } }); registerFunction(new AbstractFunction.B("inrange_integ", -1) { public boolean evalB(final EvalEnv env, final Term[] args) { final int n1 = args.length / 3; final int n2 = 2 * n1; double vSum = 0.0; double v1Sum = 0.0; double v2Sum = 0.0; for (int i = 0; i < n1; i++) { vSum += args[i].evalD(env); v1Sum += args[i + n1].evalD(env); v2Sum += args[i + n2].evalD(env); if (vSum < v1Sum || vSum > v2Sum) { return false; } } return true; } }); }
From source file:fr.amap.lidar.amapvox.jeeb.archimed.raytracing.voxel.DirectionalTransmittance.java
public double directionalTransmittance(Point3d origin, Vector3d direction) { List<Double> distances = distToVoxelWallsV2(origin, direction); //we can optimize this by storing the angle value to avoid repeating this for each position double directionAngle = FastMath.toDegrees(FastMath.acos(direction.z)); double dMoy;//from w ww. jav a 2 s . com Point3d pMoy; double d1 = 0; double transmitted = 1; for (Double d2 : distances) { double pathLength = d2 - d1; dMoy = (d1 + d2) / 2.0; pMoy = new Point3d(direction); pMoy.scale(dMoy); pMoy.add(origin); pMoy.sub(min); int i = (int) Math.floor(pMoy.x / voxSize.x); int j = (int) Math.floor(pMoy.y / voxSize.y); int k = (int) Math.floor(pMoy.z / voxSize.z); if (i < 0 || j < 0 || k < 0 || i >= splitting.x || j >= splitting.y || k >= splitting.z) { if (toricity) { while (i < 0) { i += splitting.x; } while (j < 0) { j += splitting.y; } while (i >= splitting.x) { i -= splitting.x; } while (j >= splitting.y) { j -= splitting.y; } if (k < 0 || k >= splitting.z) { break; } } else { break; } } // Test if current voxel is below the ground level if (pMoy.z < mnt[i][j]) { transmitted = 0; } else { if (Float.isNaN(voxels[i][j][k].padBV)) { //test //voxels[i][j][k].padBV = 3.536958f; return Double.NaN; } float coefficientGTheta = (float) direcTrans.getGThetaFromAngle(directionAngle, true); //input transmittance //double transmittedBefore = transmitted; transmitted *= Math.exp(-coefficientGTheta * voxels[i][j][k].padBV * pathLength); //output transmittance //double transmittedAfter = transmitted; //intercepted transmittance //double interceptedTrans = transmittedAfter - transmittedBefore; //transmitted *= Math.exp(-0.5 * voxels[i][j][k].padBV * pathLength)/*(default coeff)*/; } if (transmitted <= EPSILON && toricity) { break; } d1 = d2; } return transmitted; }
From source file:lambertmrev.Lambert.java
public double x2tof2(double x, int N, double m_lambda) { double a = 1.0 / (1.0 - x * x); double tof;// ww w .j a v a 2 s .co m if (a > 0) //ellipse { double alfa = 2.0 * FastMath.acos(x); double beta = 2.0 * FastMath.asin(FastMath.sqrt(m_lambda * m_lambda / a)); if (m_lambda < 0.0) beta = -beta; tof = ((a * FastMath.sqrt(a) * ((alfa - FastMath.sin(alfa)) - (beta - FastMath.sin(beta)) + 2.0 * FastMath.PI * N)) / 2.0); } else { double alfa = 2.0 * FastMath.acosh(x); double beta = 2.0 * FastMath.asinh(FastMath.sqrt(-m_lambda * m_lambda / a)); if (m_lambda < 0.0) beta = -beta; tof = (-a * FastMath.sqrt(-a) * ((beta - FastMath.sinh(beta)) - (alfa - FastMath.sinh(alfa))) / 2.0); } return tof; }
From source file:lambertmrev.Lambert.java
public double x2tof(double x, int N, double m_lambda) { double battin = 0.01; double lagrange = 0.2; double dist = FastMath.abs(x - 1); double tof;// w ww . ja v a2s .c o m if (dist < lagrange && dist > battin) { // We use Lagrange tof expression tof = x2tof2(x, N, m_lambda); return tof; } double K = m_lambda * m_lambda; double E = x * x - 1.0; double rho = FastMath.abs(E); double z = FastMath.sqrt(1 + K * E); if (dist < battin) { // We use Battin series tof expression double eta = z - m_lambda * x; double S1 = 0.5 * (1.0 - m_lambda - x * eta); double Q = hypergeometricF(S1, 1e-11); Q = 4.0 / 3.0 * Q; tof = (eta * eta * eta * Q + 4.0 * m_lambda * eta) / 2.0 + N * FastMath.PI / FastMath.pow(rho, 1.5); return tof; } else { // We use Lancaster tof expresion double y = FastMath.sqrt(rho); double g = x * z - m_lambda * E; double d = 0.0; if (E < 0) { double l = FastMath.acos(g); d = N * FastMath.PI + l; } else { double f = y * (z - m_lambda * x); d = FastMath.log(f + g); } tof = (x - m_lambda * z - d / y) / E; return tof; } //return tof; }
From source file:fr.cs.examples.bodies.Phasing.java
/** Guess an initial orbit from theoretical model. * @param date orbit date/*from w ww .j a v a 2 s.co 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); }