/*
  fallingclimber : calculates how a climber falling from a space elevator
  will hit the earth.
  Copyright (C) 2003 Blaise Gassend <blaise at gassend dot com> 
  
  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.
*/

#include <stdio.h>
#include <math.h>

#define STEP (1) 
#define MAXITERS (1000000)
#define RK2

#define EARTH_PULSATION (2 * M_PI / 86164.1)
#define EARTH_R (40e6 / (2 * M_PI))
#define EARTH_GM (9.81 * EARTH_R * EARTH_R)

typedef struct
{
  double t;
  double x;
  double y;
  double xd;
  double yd;
} state;

int dispflag = 0;

void getder(state *s, state *ds);
void update(state *dst, state *src, state *ds, double step);

double deg(double r)
{
  return r * 180 / M_PI;
}

void simul(double r0)
{
  int i;
  state s, ts, ds;
  
  s.x = 0;
  s.y = r0;
  s.xd = r0 * EARTH_PULSATION;
  s.yd = 0;
  s.t = 0;
  
  for (i = 0; i < MAXITERS; i++)
  {
    double step = STEP;
    double r;
    double theta;
    double longitude;
    double vtheta;
    double v;
    double vx;
    double vy;
    
#ifdef EULER
    getder(&s, &ds);
    update(&s, &s, &ds, step);
#else // EULER
#ifdef RK2
    getder(&s, &ds);
    update(&ts, &s, &ds, step / 2);
    getder(&ts, &ds);
    update(&s, &s, &ds, step);
#else // RK2
#ifdef RK4
#error RK4 not defined yet
#else // RK4
#error No integration algorithm specified
#endif // RK4
#endif // RK2
#endif // EULER
  
    r = sqrt(s.x * s.x + s.y * s.y);
    theta = atan2(-s.x, -s.y) + M_PI;
    longitude = theta - s.t * EARTH_PULSATION;
    vx = s.xd - EARTH_PULSATION * s.y;
    vy = s.yd + EARTH_PULSATION * s.x;
    v = sqrt(vx * vx + vy * vy);
    vtheta = atan2(-vx, -vy) + M_PI / 2;
    
    if (dispflag == 1)
    {
      printf("%f %f %f %f\n", 
	  s.t, r, theta, longitude);
      fflush(stdout);
    }

    if (r < EARTH_R)
    {
      if (dispflag == 2)
      {
	printf("%f %f %f %f %f %f\n", r0 - EARTH_R, s.t, deg(theta), deg(longitude), v, deg(vtheta - theta));
	fflush(stdout);
      }
      printf("#too low %f < %f\n", r, EARTH_R);
      return;
    }

    if (theta > M_PI)
    {
      printf("#timed out\n\n");
      return;
    }
  }
}
    
void getder(state *s, state *ds)
{
  double r3 = pow(s->x * s->x + s->y * s->y, 1.5);
  
  ds->x = s->xd;
  ds->y = s->yd;
  ds->xd = -EARTH_GM / r3 * s->x; 
  ds->yd = -EARTH_GM / r3 * s->y; 
  ds->t = 1;
}

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;
}

int main()
{
  double r0;
  
  dispflag = 2;

  for (r0 = EARTH_R; r0 < 30e6; r0+=r0 < 2.9e7 ? 1e5 : 1e3)
    simul(r0);
  
  return 0;
}