/*
 * plot program for VECR/iPlot
 */
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <math.h>
 
typedef double (*Func)( double);

typedef struct { double x, y; } Point;

typedef Point (*Param)( double);

#if P1 /* array of y=func(x) */

  extern Func f[];

#elif P2 /* array of parameterized (x,y) = func(t) */

  extern Param f[];

#elif P3 /* Newton's method, function f(x) and derivative fp(x) */

  extern double f(double);

  extern double fp(double);

#endif

#define Npts 101
 
static double dt, tmin, tmax, xmin, xmax, ymin, ymax;

#if P1 || P3
static void plot_func( Func f)
{
  int i;
  double x, y;

  for( i = 0; i < Npts; ++i)
  {
    x = tmin + i*dt;
 
    y = f(x);
 
    if( y == HUGE_VAL) // skip HUGE_VAL points
    {
      putchar( '\n');
    }
    else
    {
      printf( "%g %g\n", x, y);
 
      if( y < ymin) ymin = y;
      if( y > ymax) ymax = y;
    }
  }
 
  putchar( '\n');
}
#elif P2
static void plot_param( Param f)
{
  int i;
  double t;
  Point r;

  for( i = 0; i < Npts; ++i)
  {
    t = tmin + i*dt;
 
    r = f(t);
 
    if( r.x == HUGE_VAL || r.y == HUGE_VAL) // skip HUGE_VAL points
    {
      putchar( '\n');
    }
    else
    {
      printf( "%g %g\n", r.x, r.y);
 
      if( r.y < ymin) ymin = r.y;
      if( r.y > ymax) ymax = r.y;

      if( r.x < xmin) xmin = r.x;
      if( r.x > xmax) xmax = r.x;
    }
  }
 
  putchar( '\n');
}
#endif

int main( int argc, char *argv[])
{
  char *progname = argv[0];
  int xopt = 1, yopt = 1;

#if P3
  int nflag = 0;
  double nstart;
#endif
 
  /* check for -n, -x, and -y options */
 
#if P3
  if( argc > 1 && strcmp( argv[1], "-n") == 0)
  {
    nflag = 1;
    nstart = atof( argv[2]);
    argc -= 2, argv += 2;
  }
#endif
 
  if( argc > 1 && strcmp( argv[1], "-x") == 0)
  {
    xopt = 0;
    --argc, ++argv;
  }
 
  if( argc > 1 && strcmp( argv[1], "-y") == 0)
  {
    yopt = 0;
    --argc, ++argv;
  }
 
  /* now need tmin tmax args */
 
  if( argc != 3)
  {
    fprintf( stderr, "Usage: %s [-n nstart] [-x] [-y] tmin tmax\n", progname);
    return 1;
  }
 
  tmin = atof( argv[1]);
 
  tmax = atof( argv[2]);
 
  dt = (tmax - tmin)/(Npts - 1);
 
  ymin = HUGE_VAL;
  ymax = -HUGE_VAL;

#if P1 || P3
  xmin = tmin;
  xmax = tmax;
#elif P2
  xmin = HUGE_VAL;
  xmax = -HUGE_VAL;
#endif

  /* create the plot(s), track ymin and ymax */
 
#if P1
  {
    int i;

    for( i = 0; f[i]; ++i)
      plot_func( f[i]);
  }
#elif P2
  {
    int i;

    for( i = 0; f[i]; ++i)
      plot_param( f[i]);
  }
#elif P3
  {
    plot_func( f);
  }
#endif

  /* display x and y axes */
 
  if( xopt)
    printf( "\n%g %g\n%g %g\n", xmin, 0.0, xmax, 0.0);
 
  if( yopt)
    printf( "\n%g %g\n%g %g\n", 0.0, ymin, 0.0, ymax);
 
#if P3 /* Newton's method */

  if( nflag)
  {
    int i;
    double x, f0, f1;

    FILE *fout = fopen( "newton.out", "w");
    if( fout == NULL)
    {
      perror( "newton.out");
      return 1;
    }

    x = nstart;
    f0 = f(x);
    f1 = fp(x);

    printf( "\n%g %g\n%g %g\n", x, 0.0, x, f0);

    fprintf( fout, "%25s %25s %25s\n", "x", "f", "f'");
    fprintf( fout, "%25.18g %25.18g %25.18g\n", x, f0, f1);

    for( i = 0; i < 7; ++i)
    {
      x -= f0/f1;
      f0 = f(x);
      f1 = fp(x);
      printf( "%g %g\n%g %g\n", x, 0.0, x, f0);
      fprintf( fout, "%25.18g %25.18g %25.18g\n", x, f0, f1);
    }

    fclose( fout);
  }
#endif

  return 0;
}