//
// Programmer:    Craig Stuart Sapp <craig@ccrma.stanford.edu>
// Programmer:    Juan Reyes <juanig@ccrma.stanford.edu>
// Creation Date: Wed May 17 16:05:38 PDT 2000
// Last Modified: Fri May 19 00:19:23 PDT 2000
// Filename:      ...sig/doc/examples/all/henontune/henon.cpp
// Syntax:        C++ 
// $Smake:        g++ -O -o %b %f && strip %b
//
// Description:   This program generates a sequence of values
//                based on the 2-D invertible iterated map with
//                chaotic solutions proposed by the French astronomer 
//                Michel Henon ["A two-dimensional mapping with a strange 
//                attractor," Commun. Math. Phys. 50, 69-77
//                (1976)].
//
//                The equations for generating the sequence are:
//
//                y[n] = x[n-1]
//                and
//                x[n] = 1 + a x[n-1]^2 + b y[n-1]
// 
//                Where y[n] is the current output and x[n] is the
//                current input.
//

#include <iostream.h>
#include <stdlib.h>

void generateHenonValues(double alpha, double beta, int count, 
      double x = 0.0, double y = 0.0);

int main(int argc, char* argv[]) {

   if (argc < 4) {
      cout << "Usage: "<< argv[0] << " Avalue Bvalue repetitions x0 y0" << endl;
      cout << "   Suggested ranges/values: \n";
      cout << "      A = -2.5 to -1\n";
      cout << "      B = -0.5 to 0.5\n";
      cout << "      x0 = 0.0 or 0.63135448\n";
      cout << "      y0 = 0.0 or 0.18940634\n";
      cout << endl;
      exit(1);
   }
   
   double alpha = strtod(argv[1], NULL);
   double beta = strtod(argv[2], NULL);
   int count = atoi(argv[3]);

   double x = 0.0;
   double y = 0.0;

   if (argc >= 5) {
      x = strtod(argv[4], NULL);
   }
   if (argc >= 6) {
      y = strtod(argv[5], NULL);
   }

   generateHenonValues(alpha, beta, count, x, y);

   return 0;
}

///////////////////////////////////////////////////////////////////////////

void generateHenonValues(double alpha, double beta, int count, 
      double x, double y) {
   double newx;
   double newy;
   for (int i=0; i<count; i++) {
      cout << x << '\n';
      newx = 1 + alpha * x * x + beta * y;
      newy = x;
 
      x = newx; 
      y = newy;
   }
}