/* -*- Mode:C; Coding:us-ascii-unix; fill-column:132 -*- */
/**********************************************************************************************************************************/
/**
   @file      lorenz2.c
   @author    Mitch Richling <http://www.mitchr.me/>
   @Copyright Copyright 1994,1997,2001,2014 by Mitch Richling.  All rights reserved.
   @Revision  $Revision: 1.2 $ 
   @SCMdate   $Date: 2014/10/10 19:40:05 $
   @brief     Generate a PovRay include file for the lorenz strange attractor. @EOL
   @Keywords  lorenz strange attractor povray adaptive Euler
   @Std       C89

   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.  Alternate implimentation lorenz2.pl.
              
***********************************************************************************************************************************/

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

#define DEBUG 0
#define CYL 1

int main () {  
  double distToGo = 2500.0;             /* How long should the curve be? */

  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. */

  /* 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 curMaxTdelta, curMinTdelta, tDelta, dist, Xdelta, Ydelta, Zdelta, movDelta, xOld, yOld, zOld;
  int    numBisect, doneBisecting, numBalls, numOver;
  char   moveCh;

 /*  Solve the equations..... */
  numBalls = 0;
  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 (fabs (tDelta - curMinTdelta) < tDeltaZeroThresh) {
			doneBisecting = 1;
		  } else {
			curMaxTdelta = tDelta;
			tDelta = (tDelta + curMinTdelta) / 2;
			if (tDelta < minTdelta) {
			  tDelta = minTdelta;
			}
		  }
		} else if (movDelta < minMovDelta) {
		  if (fabs (tDelta - curMaxTdelta) < tDeltaZeroThresh) {
			doneBisecting = 1;
		  } else {
			curMinTdelta = tDelta;
			tDelta = (curMaxTdelta + tDelta) / 2;
			if (tDelta > maxTdelta) {
			  tDelta = maxTdelta;
			}
		  }
		} else {
		  doneBisecting = 1;
		}
      }
    }
    dist += movDelta;
    xOld = x;
    yOld = y;
    zOld = z;
    x = x + Xdelta;
    y = y + Ydelta;
    z = z + Zdelta;
    if (movDelta > maxMovDelta) {
      moveCh = '*';
      numOver++;
    } else {
      moveCh = ' ';
    }
#if DEBUG
    fprintf (stderr, "%d 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",
			 numBisect, 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);
#if CYL
    printf("cylinder {<%f,%f,%f>, <%f,%f,%f>, %f texture { crvTex } }\n", x, y, z, xOld, yOld, zOld, sphereRad);
#endif

  }
  fprintf (stderr, "Balls: %7d Dist: %7.1f #Over: %7d\n", numBalls, dist, numOver);

  return 0;
}