/*
   Author:   Mitch Richling<http://www.mitchr.me>
   IP:       Copyright 2001 by Mitch Richling.  All rights reserved.
   Key word: 
   Notes:    This is an adaptive Euler's method to solve Lorenz's
             differential equations.  The value of delta f is bounded
             above and below so that the curve looks smooth when we
             are done.  The bounding is adjusted by adjusting t delta
             in each leap-frog step.  The value of delta t is bounded
             absolutely above and below to keep things from going
             nuts.  The method of finding an acceptable t delta is to
             use bisection between the bounds on t delta.  This is a
             strange adaptation of Euler's method, but it works well
             for this application. 
*/

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

#define DEBUG 0

int main () {  
  double distToGo = 2500.0; /* How long should the curve be? */
  double dist;              /* Total length of curve. */

  double sphereRad = 0.2;   /* Size of the spheres to use */

  double maxMovDelta = sphereRad / 4.0; /* The maximum delta for the function. */
  double minMovDelta = sphereRad / 5.0; /* The minimum delta for the function. */

  double minTdelta = 0.000001;          /* Minimum value for tDelta */
  double maxTdelta = 0.01;              /* Maximum value for tDelta */
  double tDeltaZeroThresh = 0.00000001; /* Episolon ball for zero. */
  double maxNumBisect = 10;             /* Max bisections on tDelta. */
  double curMaxTdelta, curMinTdelta;    /* Current bisection window. */
  int    numBisect, doneBisecting;      /* Bisecting vars. */
  double tDelta;                        /* Current tDelta value. */

  /* Initial conditions, and constants. */
  double x = 0.1;
  double y = 0.0;
  double z = 0;
  double a = 10;
  double b = 28;
  double c = 8.0 / 3;

  double Xdelta, Ydelta, Zdelta, movDelta; /* Each Leap-Frog jump. */
  int numBalls, numOver, numBisect, doneBisecting;
  char moveCh;

 /*  Solve the equations..... */
  tDelta = maxTdelta;
  dist = 0.0;
  numOver = 0;
  while (dist < distToGo) {
    numBalls++;
    /*  Take a big step up to minimize the number of balls. */
    tDelta = (maxTdelta + tDelta) / 2;
    if (tDelta > maxTdelta) {
      tDelta = maxTdelta;
    }
    curMaxTdelta = maxTdelta;
    curMinTdelta = minTdelta;
    /*  Bisect t delta until we have done it too much or */
    /*  until we have a good value for t delta. */
    numBisect = 0;
    doneBisecting = 0;
    while (!(doneBisecting)) {
      numBisect++;
      Xdelta = a * (y - x) * tDelta;
      Ydelta = (x * (b - z) - y) * tDelta;
      Zdelta = (x * y - c * z) * tDelta;
      movDelta = sqrt (fabs (Xdelta * Xdelta + Ydelta * Ydelta + Zdelta * Zdelta));
      if (numBisect > maxNumBisect) {
        doneBisecting = 1;
      } else {
        if (movDelta > maxMovDelta) {
          if (abs (tDelta - curMinTdelta) < tDeltaZeroThresh) {
            doneBisecting = 1;
          } else {
            curMaxTdelta = tDelta;
            tDelta = (tDelta + curMinTdelta) / 2;
            if (tDelta < minTdelta) {
              tDelta = minTdelta;
            }
          }
        } else if (movDelta < minMovDelta) {
          if (abs (tDelta - curMaxTdelta) < tDeltaZeroThresh) {
            doneBisecting = 1;
          } else {
            curMinTdelta = tDelta;
            tDelta = (curMaxTdelta + tDelta) / 2;
            if (tDelta > maxTdelta) {
              tDelta = maxTdelta;
            }
          }
        } else {
          doneBisecting = 1;
        }
      }
    }
    dist += movDelta;
    x = x + Xdelta;
    y = y + Ydelta;
    z = z + Zdelta;
    if (movDelta > maxMovDelta) {
      moveCh = '*';
      numOver++;
    } else {
      moveCh = ' ';
    }
#if DEBUG
    fprintf (stderr, "Balls: %7d tDel: %9.6f x: %9.2f y: %9.2f z: %9.2f Dist: %7.1f fDel: %9.3f #Over: %7d %1c-%1c\n",
             numBalls, tDelta, x, y, z, dist, movDelta, numOver, moveCh, (numBisect > maxNumBisect ? '*' : ' '));
#endif
    printf ("sphere {<%f,%f,%f>, %f texture { crvTex } }\n", x, y, z, sphereRad);
  }
  fprintf (stderr, "Balls: %7d Dist: %7.1f #Over: %7d\n", numBalls, dist, numOver);
}