/**
 * Simple calculation of the trajectory of a climber falling off of a space elevator
 * based on the C program by Blaise Gassend
 * translated to Java by David Forslund (Los Alamos)
 *
  This program is free software; you can redistribute it and/or
  modify it under the terms of the GNU General Public License
  as published by the Free Software Foundation; either version 2
  of the License, or (at your option) any later version.

  This program is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with this program; if not, write to the Free Software
  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
 */

public class FallingClimber {
    double pi = Math.PI;
    double earthAngularRotationRate = 2 * pi / 86164.1;   // radians per second

    double earthRadius = 40.0E6 / (2 * pi);   // in meters

    double earthGM = 9.81 * earthRadius * earthRadius;
    // time step in seconds
    int STEP = 1;
    //   int MAXITERS = 1000000;
    int dispflag = 0;
    // configure the format for displaying the results
    java.text.DecimalFormat format = new java.text.DecimalFormat();

    // initializer flag for displaying labels for the output
    boolean init;

    public FallingClimber() {
        init = false;
    }

    /**
     * Convert radians to degreees
     */
    public double deg(double r) {
        return r * 180 / pi;
    }

    /**
     * simulator
     * @param r0 is the starting radious
     * @return  the ending radius
     */
    public double simulate(double r0) {

        double step = STEP;
        double r;
        double theta;
        double longitude;
        double vtheta;
        double v;
        double vx;
        double vy;
        State s = new State();

        s.x = 0;
        s.y = r0;
        s.xd = r0 * earthAngularRotationRate;
        s.yd = 0;
        s.t = 0;

        //for (i = 0; i < MAXITERS; i++) {
        while (true) {
           // euler(s, step);
            rk2(s, step);
            r = Math.sqrt(s.x * s.x + s.y * s.y);
            theta = Math.atan2(-s.x, -s.y) + pi;
            longitude = theta - s.t * earthAngularRotationRate;
            vx = s.xd - earthAngularRotationRate * s.y;
            vy = s.yd + earthAngularRotationRate * s.x;
            v = Math.sqrt(vx * vx + vy * vy);
            vtheta = Math.atan2(-vx, -vy) + pi / 2;

            if (dispflag == 1) {
                System.out.println(s.t + " " + r + " " + theta + " " + longitude);
            }
            if (dispflag == 2 && !init) {
                System.out.println("Height(m)\tTime(s)\tAngle\tLongitude\tVelocity(m/s)\tAngleOfImpact");
                init = true;
            }
            if (r < earthRadius) {
                // The climber has hit the earth service
                if (dispflag == 2) {
                    System.out.println(format.format(r0 - earthRadius) + "\t" + format.format(s.t) + "\t" +
                            format.format(deg(theta)) + "\t" + format.format(deg(longitude)) + "\t" + format.format(v) + "\t" +
                            format.format(deg(vtheta - theta)));
                }
                //  System.out.println("#too low "+r+" < "+earthRadius);
                return r;
            }

            if (theta > pi) {
                // perigee has been reached, so the climber is actually in orbit
                return r;
            }


        }
        //System.err.println("failed to complete!");
    }

    /**
     * Compute the derivative
     * @param s contains input position and velocity
     * @param ds returns velocity and acceleration
     */

    private void getder(State s, State ds) {
        double r3 = Math.pow(s.x * s.x + s.y * s.y, 1.5);

        ds.x = s.xd;
        ds.y = s.yd;
        ds.xd = -earthGM / r3 * s.x;
        ds.yd = -earthGM / r3 * s.y;
        ds.t = 1;
    }

    /**
     * udate the position and velocity with the velocity and acceleration
     *
     */

    private void update(State dst, State src, State ds, double step) {
        dst.x = src.x + ds.x * step;
        dst.y = src.y + ds.y * step;
        dst.xd = src.xd + ds.xd * step;
        dst.yd = src.yd + ds.yd * step;
        dst.t = src.t + ds.t * step;
    }

    /**
     * Euler integration step
     * @param s
     * @param step
     */
    private void euler(State s, double step) {
        State ds = new State();
        getder(s, ds);
        update(s, s, ds, step);
    }

    /**
     * runge kutte integration step
     * @param s
     * @param step timestep
     */
    private void rk2(State s, double step) {
        State ds = new State();
        getder(s, ds);
        State ts = new State();
        update(ts, s, ds, step / 2);
        getder(ts, ds);
        update(s, s, ds, step);
    }

    /**
     * main program to call the simulator
     * @param argv
     */
    public static void main(String[] argv) {
        double r0;

        FallingClimber c = new FallingClimber();
        c.dispflag = 2;
        for (r0 = c.earthRadius; r0 < 30e6; r0 += r0 < 2.9e7 ? 1e5 : 1e3)
            c.simulate(r0);

    }

    /**
     * State variables for climber can hold position and velocity as
     * as well as the velocity and acceleration
     */

    private class State {
        double t;
        double x;
        double y;
        double xd;
        double yd;
    }
}