/*

  CCATSL Example program: Ordinary Differential Equations

  D.N.Harris@damtp.cam.ac.uk
  2002-04-04

  Illustrates use of Rk4CL, Rkf45CL to solve planetary orbits

  You can investigate the effects of how changing step length in Rk4CL
  affects discretization errors and long-term errors caused by round-off.

  See how changing RELERR affects Rkf45CL

  Planetary motion equations are:
  
  x''(t) = -G*x(t)/R(t)
  y''(t) = -G*y(t)/R(t)
  where R(t) = (x^2+y^2)^(3/2)
      
If G=1 and x=1, y=0, x'=0, y'=1 at t=0 we get a circular orbit
with period 2*pi.
  
*/

#include <catam.h>

#define YEARS           1 
#define ORDER           4 /* system of 4 ODEs */

/* number of integration steps per year for fixed-step length Rk4CL */
#define MAXSTEPS        180

/* required errors in Rkf45CL which adjusts step-length automatically */ 
#define ABSERR          0.0
#define RELERR          1e-9

#define X_INIT       1.0
#define Y_INIT       0.0
#define XDASH_INIT   0.0
#define YDASH_INIT   1.0

/* ranges for x and y */
#define SIZE  1.0
#define XMIN -SIZE
#define XMAX  SIZE
#define YMIN -SIZE
#define YMAX  SIZE

/* If initial parameters are changed so that the orbit is no longer
   a circle you will need to calculate appropriate scales for x and y,
   either by updating minimum and maximum values of x and y as points
   are calculated, or by storing the point coordinates in arrays
   and plotting with XYCurveCL using the AUTOAXES parameter.
*/

/* If G is to be varied its value must be passed to f as an external variable:
   
   double G = 1.0;
   
   but I'll just define it as a constant */

#define G           1.0

/* function f defines the system of equations */

void f(double t, double u[], double udash[])
{
  double R;
  /* To rewrite the planetary motion equations as a system of first order
     equations let
     
     u[0] = x
     u[1] = y
     u[2] = x'
     u[3] = y'
  */
  
  R = u[0]*u[0] + u[1]*u[1];
  R = R * sqrt(R);
  udash[0] =  u[2];
  udash[1] =  u[3];
  udash[2] = -G * u[0] / R;
  udash[3] = -G * u[1] / R;
}

int
MainCL(void)
{
  double u[ORDER], udash[ORDER];
  double t, deltat, yearend;
  int step, maxsteps = MAXSTEPS, year, years = YEARS;
  int rkf_control;

  int menu_option;
  int solver = 1; /* 1=Rk4CL, 2=Rkf45CL */
  int display_points = -1; /* all points; 1 = year end points only */

  XRangeCL(XMIN, XMAX);
  YRangeCL(YMIN, YMAX);

  do {
    menu_option = MenuCL("Options ",
			 "Select Rk4CL$"
			 "Select Rkf45CL$"
			 "Show all points$"
			 "Show only end of year points$"
			 "Calculate$"
			 "quit");

    if (menu_option >= 6) break;

    switch (menu_option) {
    case 1:
      MessageCL("Solving with Rk4CL");
      solver = 1;
      break;

    case 2:
      MessageCL("Solving with Rkf45CL");
      solver = 2;
      break;

    case 3:
      XRangeCL(-1, 1);
      YRangeCL(-1, 1);
      display_points = -1;
      break;

    case 4:
      display_points = 1;
      XRangeCL(.9, 1.1);
      YRangeCL(-.1, .1);
      break;

    case 5:
      do {
	years = ReadIntCL("Number of years (0 to finish)", years);
	if (years <= 0) break;
	/* Initial values */
	t = 0;
	u[0] = X_INIT;
	u[1] = Y_INIT;
	u[2] = XDASH_INIT;
	u[3] = YDASH_INIT;
	
	WClearCL();
	GAxesCL();
	GLineColourCL(BlueCC);
	
	if (solver == 1) { /* use Rk4CL */
	  maxsteps = ReadIntCL("Time steps per year", maxsteps);
	  deltat = 2 * M_PI / maxsteps;
	  for (year = 1; year <= years; ++year) {
	    for (step = 1; step <= maxsteps ; ++step) {
	      Rk4CL(ORDER, &t, deltat, f, u, udash);
	      if (display_points == -1 ) { /* plot all points */
		XYSymbolCL(u[0], u[1], XCROSS);
	      }
	      /* Press Esc to interrupt lengthy calculation */
	      if (EscapeCL()) HaltCL();
	    }
	    XYSymbolCL(u[0], u[1], XCROSS); /* plot the year end point */
	    printf("\n%7.3f%7.3f%7.3f\n", t, u[0], u[1]);
	  }
	} else { /* Use Rkf45CL */
	  t = 0;
	  for (year = 1; year <= years; ++year) {
	    yearend = year * 2 * M_PI; /* integrate to end of year */
	    /* set rkf_control to -1 to return after each intermediate
	       point is found, 1 to return after end of interval */ 
	    rkf_control = display_points;
	    do {
	      Rkf45CL(ORDER, ABSERR, RELERR, yearend, &t, f, u, udash, &rkf_control);
	      XYSymbolCL(u[0], u[1], 'X');
	      if (EscapeCL()) HaltCL();
	    } while (rkf_control == -2); /* repeat for all intermediate values */
	    if (rkf_control != 2) {
	      /* Something is wrong - see the CCATSL manual */
	      printf("Error: rkf_control = %i\n", rkf_control);
	      HaltCL();
	    }
	    printf("%7.3f%7.3f%7.3f\n", t, u[0], u[1]);
	  } /* for year */
	} /* solver */
      } while (years > 0);
    } /* switch */
  } while (menu_option < 6);
} /* end MainCL */