#!/usr/local/bin/perl
#
#  Author:   Mitch Richling<http://www.mitchr.me>
#  IP:       Copyright 1994 by Mitch Richling.  All rights reserved.
#  Key word: lorenz strange attractor povray adaptive Euler
#  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.

# How long should the curve be?
$distToGo = 2500;

# Size of the spheres to use
$sphereRad = 0.2;

# The maximum delta for the function.
$maxMovDelta   = $sphereRad / 4;
$minMovDelta   = $sphereRad / 5;

# t Delta parameters.
$minTdelta    = 0.000001;
$maxTdelta    = 0.01;
$maxNumBisect = 10;
$tDeltaZeroThresh = 0.00000001;

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

# Solve the equations.....
$tDelta       = 0.01;
while($dist < $distToGo) {
  $numBalls++;
  # Take a step up so that we 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(abs($Xdelta**2+$Ydelta**2+$Zdelta**2));
    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; }
        }
      } elsif($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 = ' ';
  }
  
  #printf(STDERR "Balls: %7d tDel: %9.6f x: %9.2f y: %9.2f z: %9.2f Dist: %7.1f fDel: %9.3f #Over: %7d %1s-%1s\n", $numBalls, $tDelta, $x, $y, $z, $dist, $movDelta, $numOver, $moveCh, ($numBisect>$maxNumBisect?'*':' '));
  print "sphere {<$x,$y,$z>, $sphereRad texture { crvTex } }\n"; 
}
printf(STDERR "Balls: %7d Dist: %7.1f #Over: %7d\n", $numBalls, $dist, $numOver);