#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "twin_solve.h"

bool logging = false;
FILE *log_file;

#define X_HIGH_OFFSET 0.0000000001

static const char *do_what = "porthole";
static double v = 0.5;
static double x_low = 0;
static double x_high = (M_PI + X_HIGH_OFFSET);
static double x_step = (M_PI / 40);

inline bool streq (const char *a, const char *b) { return (!strcmp (a, b)); }

/****************************************************************/

static double ph_phi (double theta, double phi)
{
  double lhs = theta - phi;
  double rhs = 2 * v * sin ((theta + phi)/2);
  double diff = lhs - rhs;

  db_printf (log_file, "ph_phi(%f,%f) -> %f\n", theta, phi, diff);

  return (diff);
}  /* ph_phi */

static double theta;

static double ph_phi_func_1var (double phi)
{
  return (ph_phi (theta, phi));
}  /* ph_phi_func */

/* Print a series of lines.  Each one has theta, phi, and the derivative
   of phi with respect to theta.  Theta and phi are the angles made by
   ship A (observer) and ship B (source) as viewed from the Earth frame.
   Theta and phi are v*gamma times ship's clock time on A an B, respectively,
   if the radius of the orbit is 1, so the derivative is also the relative
   speed of B's clock as seen through a telescope on A.
 */

static void do_porthole ()
{
  theta = 0;
  double phi = 0;
  for (theta = x_low; theta < x_high; theta += x_step)
    {
      db_printf (log_file, "do_porthole(theta=%f)\n", theta);

      phi = solve_func (ph_phi_func_1var, 0, phi); /* Use the old value of phi for the guess */

      double vcos = v * cos((theta + phi) * 0.5);
      double dphi_dtheta = (1 - vcos) / (1 + vcos);

      printf ("%f	%f	%f\n",
	      theta, phi, dphi_dtheta);
    }  /* for */
}  /* do_porthole */

/****************************************************************/
/* Similar but for the "simultaneous" (Lorentz-transform) view */

static double sim_phi (double phi)
{
  double lhs = theta - phi;
  double rhs = v * v * sin (theta + phi);
  double diff = lhs - rhs;

  db_printf (log_file, "sim_phi(%f,%f) -> %f\n", theta, phi, diff);

  return (diff);
}  /* sim_phi */

static void do_simultaneous ()
{
  theta = 0;
  double phi = 0;
  for (theta = x_low; theta < x_high; theta += x_step)
    {
      db_printf (log_file, "do_simultaneous(theta=%f)\n", theta);

      phi = solve_func (sim_phi, 0, phi); /* Use the old value of phi for the guess */

      double cos_tp = v * v * cos(theta + phi);
      double dphi_dtheta = (1 - cos_tp) / (1 + cos_tp);

      printf ("%f	%f	%f\n",
	      theta, phi, dphi_dtheta);
    }  /* for */
}  /* do_simultaneous */

/****************************************************************/

main (int argc, char **argv)
{
  bool fix_step = false;
  int steps = 40;

  for (int i = 1; i < argc; i++)
    {
      const char *arg = argv[i];
      if (streq (arg, "-do"))
	do_what = argv[++i];
      else if (streq (arg, "-v"))
	v = atof (argv[++i]);
      else if (streq (arg, "-steps"))
	{
	  steps = atoi (argv[++i]);
	  fix_step = true;
	}
      else if (streq (arg, "-to"))
	{
	  x_high = atof (argv[++i]);
	  x_high += X_HIGH_OFFSET;
	  fix_step = true;
	}
      else if (streq (arg, "-from"))
	{
	  x_low = atof (argv[++i]);
	  fix_step = true;
	}
      else if (streq (arg, "-db"))
	{
	  logging = true;
	  log_file = fopen ("twins_debug.log", "w");
	  db_printf (log_file, "twins starting\n");
	}
      else
	{
	  fprintf (stderr, "Argument %s not understood\n", arg);
	  exit(1);
	}
    }  /* for */

  if (fix_step)
    {
      x_step = x_high - x_low;
      x_step -= X_HIGH_OFFSET;	/* Shade it so we seem to hit the end exactly */
      x_step /= steps;
    }

  if (streq (do_what, "porthole"))
    do_porthole();
  else if (streq (do_what, "simultaneous")
	   || streq (do_what, "sim"))
    do_simultaneous();
  else
    {
      fprintf (stderr, "Don't know how to do this: %s\n", do_what);
      fprintf (stderr, "Must be one of 'porthole'\n");
      exit(1);
    }

  db_printf (log_file, "Done.\n");
  exit (0);
}  /* main */