EM Drive acceleration estimator

package emcalc;

import java.util.function.Function;

/**
 * Derive EMDrive acceleration from NASA figures.
 *
 * @author ps200306
 * @email [email protected]
 * @since 17-Jan-2017
 *
 */
public class Emcalc {

    // Constants

    final double G = 6.67408e-11; // Gravitational constant, units m3 kg-1 s-2.
    final double M = 1.989e30; // Mass of Sun in kg.
    final double GM = G * M; // convenience constant.
    final double Re = 1.496e11; // Sun-Earth distance in metres.
    final double Rp = 5.906380e12; // Sun-Pluto distance in metres.
    final double Rx = 4.014e16; // Sun-Proxima distance in metres.
    final double yr = 365 * 24 * 3600; // year in seconds
    final double ft = yr * 1.5; // target flight time = 18 months
    final double v0 = 0; // initial radial velocity is zero (just added for completeness)

    double a0 = 0; // The EMDrive acceleration which we seek.


    public static void main(String args[]) {
        (new Emcalc()).calc();
    }


    /**
     * Main calculation. Take initial estimate of EMdrive acceleration
     * (in the range zero to something big) then binary chop until we
     * get within desired accuracy of the target flight time.
     */
    void calc() {

        final double precision = 0.01; // target precision 1%
        double error = 1; // try to reduce this error to within target precision
        final int numSteps = 100 * 1000; // divisions for numerical integration

        // initial range of EMDrive acceleration in m/s^2
        double aLo = 0;
        double aHi = 1000; // something big

        System.out.printf("target flight time=%s, in years = %s\n\n", ft, ft/yr);

        while (error > precision) {

            a0 = (aLo + aHi) / 2;

            final double aconst = a0; // final needed for Java 8 lambda

            // Perform numerical integration on 1/v
            double t = integrate(Re, Rp, numSteps, r ->
                1/Math.sqrt(2 * (v0 + aconst * r + GM/r) - 2 * (aconst * Re + GM/Re))
            );

            // Compare to required target
            error = Math.abs(1 - t / ft);

            System.out.printf("a0=%3.3e, t=%3.3e, error=%3.3e \n", a0, t, error);

            // We can get some singularities e.g. divide by zero, in which case we
            // just pretend the result was too high and let the binary chop do the rest.
            if (Double.isNaN(t)) {
                t = 2* ft;
                error = precision + 1;
            }

            // Binary chop depending on whether result was high or low.
            aLo = (t > ft)? a0 : aLo;
            aHi = (t > ft)? aHi : a0;

        }

        // Final acceleration estimate
        System.out.printf("\nFinal accelelation a0 = %s\n", a0);

        // Use this acceleration to fly to Proxima Centauri
        final double aconst = a0;
        double u = integrate(Re, Rx, numSteps * 10, r ->
            1/Math.sqrt(2 * (v0 + aconst * r + GM/r) - 2 * (aconst * Re + GM/Re))
        );
        System.out.printf("\nFlight time to Proxima = %s, in years = %s\n\n", u, u/yr);

    }

    /**
     *
     * Numerical integrator.
     *
     * @param lower - integration lower bound
     * @param upper - integration upper bound
     * @param numSteps - number of integration steps (affects precision)
     * @param integrand - function to integrate
     * @return the integral over the bounds
     */
    double integrate(double lower, double upper, int numSteps, Function<Double, Double> integrand) {

        double step = (upper - lower)/numSteps; // step size
        double integral = 0;

        for (int i = 0; i < numSteps; i++) {
            double L = lower + step * i + step/2; // use midpoint of each step
            integral += (step * integrand.apply(L));
        }

        return integral;

    }

}