/******************************************************************************/
/*                                                                            */
/*  COMPUTE - All routines for computing the ray-traced vortex                */
/*                                                                            */
/******************************************************************************/

#define STRICT
#include <windows.h>
#include <commctrl.h>
#include <assert.h>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <ctype.h>
#include "const.h"
#include "classes.h"
#include "funcdefs.h"
#include "vortex.rh"
#include "extern.h"

// This table does not depend on the eye location, and it is not changed
// after being created, so it can safely be static, even in a multi-threaded case.

static int table_length ;
static double *arc_radius, *arc_length ;


/*
-------------------------------------------------------------------------------------

   curve() - This defines the profile curve of the side of the vortex.

   It must return a value of -1.0 when radius=0.0, the bottom (center) of the vortex.
   Then it must rise as radius increases, reaching a value of 0.0 at the top,
   which is where radius=1 (the outer perimeter of the vortex).
   For best appearance, the first derivative of the curve must be 0 when radius=1.
   This implies that the edge of the vertex is flat.
   This will ensure 0<=radius<=1, and it returns Z between -1 and 0.

-------------------------------------------------------------------------------------
*/

#if 0
// This is a simple cubic
double curve ( double radius )
{
   double temp ;

   if (radius > 1.0)
      return 0.0 ;
   temp = 1.0 - radius ;
   return -temp * temp * temp ;
}
#else

// This combines two cubics for faster, steeper plunge

// These define the transition between the upper and lower cubics
#define CURVE_A 0.5
#define CURVE_B 0.5

double curve ( double radius )
{
   double x ;

   if (radius > 1.0)
      return 0.0 ;

   if (radius >= CURVE_A) {   // Upper curve
      x = 1.0 - (radius - CURVE_A) / (1.0 - CURVE_A) ;
      return -CURVE_B * x * x * x ;
      }

   else {               // Lower curve
      x = radius / CURVE_A ;
      return (1.0 - CURVE_B) * x * x * x - 1.0 ;
      }
}
#endif


/*
-------------------------------------------------------------------------------------

   make_arc_table - Create the table mapping radius to cumulative arc length
   search_arc_table - Given a radius, find the arc length from this table
   delete_arc_table - Free its memory

-------------------------------------------------------------------------------------
*/

int make_arc_table (
   int length  // Length of table; larger makes more accuracy
   )
{
   int i ;
   double eps, value ;

   table_length = length ;
   arc_radius = (double *) malloc ( table_length * sizeof(double) ) ;
   arc_length = (double *) malloc ( table_length * sizeof(double) ) ;

   if (arc_radius == NULL  ||  arc_length == NULL) {
      if (arc_radius != NULL) {
         free ( arc_radius ) ;
         arc_radius = NULL ;
         }
      if (arc_length != NULL) {
         free ( arc_length ) ;
         arc_length = NULL ;
         }
      return ERROR_INSUFFICIENT_MEMORY ;
      }

   eps = 1.0 / (table_length - 1.0) ;

   for (i=0 ; i<table_length ; i++) {
      arc_radius[i] = i * eps ;
      arc_length[i] = curve ( arc_radius[i] ) ;  // Runs from -1 to 0
      }

   // Pythagorean base is eps and height is arc_length[i] - arc_length[i-1] 
   // We must generate the table backwards because we are changing it in place.
   // When this loop is done it contains the lengths of the individual segments.

   for (i=table_length-1 ; i>0 ; i--) {
      value = arc_length[i] - arc_length[i-1] ;
      arc_length[i] = sqrt ( eps * eps + value * value ) ; // Pythagorean theorem
      }

   arc_length[0] = 0.0 ;  // By definition

   for (i=1 ; i<table_length ; i++)      // Cumulate individual segments
      arc_length[i] += arc_length[i-1] ; // i=1 here is redundant but clear

   for (i=1 ; i<table_length ; i++)      // Normalize to unit total
      arc_length[i] /= arc_length[table_length-1] ;

   return 0 ;
}


int search_arc_table ( double radius )
{
   int klo, khi, k ;

/*
   Use binary search to bound radius between klo and khi.
   This binary search algorithm assumes that the value being bounded
   cannot lie below the lowest entry.
   Both of these arrays range from 0 to 1, inclusive.
*/

   if (radius >= arc_radius[table_length-1])  // Outside range means last case
      k = table_length - 1 ;
   else {                        // We need to do binary search
      klo = 0 ;                  // Will have arc_radius[klo] <= radius < arc_radius[khi]
      khi = table_length - 1 ;
   
      while (khi > klo+1) {      // While not bounded between adjacent table entries
         k = (khi + klo) / 2 ;
         if (radius < arc_radius[k])
            khi = k ;
         else
            klo = k ;
         }
      k = klo ;
      }
   return k ;
}


void delete_arc_table ()
{
  if (arc_radius != NULL) {
      free ( arc_radius ) ;
      arc_radius = NULL ;
      }
   if (arc_length != NULL) {
      free ( arc_length ) ;
      arc_length = NULL ;
      }
}

 
/*
-------------------------------------------------------------------------------------

   find_diff - This computes the difference between the depth (Z) of the vortex
     and the value of Z obtained by projecting along the ray a distance of t.
     Thus, this is the equation whose smallest root must be found in order to
     locate the point on the object that corresponds to the given window pixel.
     Note that t is the ratio of the projected ray's length to the eye's length.

     This is called with the eye location (0, y_eye, z_eye), the window pixel
     location (row, col), and t.  The eye location is in the same units as
     the vortex. Note that without loss of generality (because the vortext
     has circular symmetry) we force the eye X to be 0.

     The caller must have z_eye be greater than 1, and y_eye not positive,
     so we are looking down on the object below us.
     In other words, the eye is always above the XY plane, and the object is
     entirely below the XY plane, except that the perimeter of the vortex
     touches the XY plane.  If y_eye is zero we are looking directly down
     on the vortex, while if y_eye is negative the vortex is below and
     in front of us.

     We let the window center (0, 0) be the object-world origin (0, 0, 0).
     The coordinates within the window typically range from about -1
     (bottom row or leftmost column) to about +1 (top row or rightmost column).
     These limits are not mandatory but are often decent.  Adjust as needed
     to get a satisfactory view of the vortex.

     The relative length of the ray from the eye, through a window pixel,
     to the object, is given by t.  When t=0 we are at the eye,
     and when t = 1 we are on the window plane.
     If y_eye = 0 the entire vortex will be behind (below) the window
     (with the perimeter on the window).  In other words, t>=1
     for all points on the vortex, with t=1 around the perimeter.
     But if y_eye < 0 part of the object will be in front of the window (t < 1).

     Two values of Z are computed by this routine.  First, it computes the Z
     coordinate of the end of the ray having relative length t.
     Then it computes the X and Y coordinates of the end of the ray and uses
     these to compute the Z of the vortex corresponding to this (X,Y) point.
     This routine returns the former minus the latter.

     This difference will be zero when the projected ray intersects with the
     boundary of the object.
     If we are looking nearly straight down (Y near zero) there will be just
     one such intersection. But if we are looking from a low angle
     there can be two or three intersections, and we are interested in only the
     closest (smallest t), because this is the one visible to the eye.
     The other two are hidden.

     It is strongly suggested that z_eye>1, which is not a practical limitation.
     Mathematical weirdness can happen otherwise.

     About passing delta as pointers... The delta vector is just scratch,
     used locally only and not needed by the caller.  So why do I pass pointers
     to it in the parameter list?  This routine (find_diff) will be called
     many, many times with the same eye and window pixel location, just varying t.
     Delta does not depend on t, so it's wasteful to keep recomputing it.
     The first time this routine is called for this eye and pixel, the caller
     would set first_call true (nonzero), and delta will be computed and saved in
     the caller's local area.  Subsequently, first_call should be set false (0)
     so instead of recomputing delta it will be fetched from the caller.
     If I rather made delta static, this routine would not be threadsafe.

-------------------------------------------------------------------------------------
*/

double find_diff (
   double y_eye ,    // <= 0
   double z_eye ,    // > 1
   double row ,      // Row location in view window
   double col ,      // And column; often between -1 and 1
   double t ,        // Relative length of ray; 0 at eye and 1 at window
   int first_call ,  // Flag: First call with this eye and window pixel?
   double *x_delta , // Work areas saved by caller
   double *y_delta ,
   double *z_delta
   )
{
   double theta, x_window, y_window, z_window ;
   double x_terminus, y_terminus, z_terminus, radius ;

/*
   The view window is the plane containing the object-world origin (0,0,0)
   that is perpendicular to the ray connecting the eye to (0,0,0).
   We are forcing x_eye to be 0, so it is easy to compute the angle theta
   of the eye above the XY plane.
   Project the vector from the eye to the window's (row, col) pixel
   and compute the object-world 3D coordinates of the intersection of this
   ray with the view window plane.  This points to what the eye sees
   for this view pixel.
*/

   if (first_call) { // Delta does not depend on t, so compute it only once
      theta = atan2 ( z_eye , -y_eye ) ; // ArcTan (z / -y) in range -Pi to Pi
      x_window = col ;
      y_window = row * sin ( theta ) ;
      z_window = row * cos ( theta ) ;

/*
   Compute the ray from the eye to this pixel on the view window.
   Recall that the X coordinate of the eye is 0.
*/

      *x_delta = x_window ;
      *y_delta = y_window - y_eye ;
      *z_delta = z_window - z_eye ;
      }

   // If z_delta is not negative, we are looking at or above the horizon
   // so the ray cannot possibly intersect the vortex.

   if (*z_delta >= 0.0)  // Projection from eye to window pixel rising!
      return 1.e60 ;     // Flag that we will never intersect (above horizon)

/*
   Adjust this ray to have a length equal to t times the length of the delta ray.
   Again, recall that the X coordinate of the eye is 0.
*/

   x_terminus = t * *x_delta ;
   y_terminus = t * *y_delta + y_eye ;
   z_terminus = t * *z_delta + z_eye ;

/*
   Use these (X,Y) coordinates to compute the Z coordinate of the vortex at (X,Y).
   Subtract this from the Z of the t-projected ray.  This difference will be zero
   when the end of the ray from the eye through the view window pixel intersects
   the vortex.  One, two, or three intersections are possible.
*/


   radius = sqrt ( x_terminus * x_terminus + y_terminus * y_terminus ) ;
   return z_terminus - curve ( radius ) ;
}


/*
-------------------------------------------------------------------------------------

   find_t - This computes t that is the smallest root of find_diff.

   If we are looking straight or almost straight down on the vortex
   (y_eye is small relative to z_eye), there is just one root (intersection).
   In this case we could use a very efficient root finder.
   But at a low angle the ray could pass through the vortex near an edge,
   then pass through again into the central part of the vortex, and finally
   pass through once more nearer the base of the vortex, exiting it.
   In this case, only the first (smallest t) is visible, so we MUST find the
   smallest root.  That's not easy.  I use a primitive method here,
   just marching along from small t to larger until the equation (find_diff)
   changes sign.
   Note that it's theoretically possible, though extremely unlikely, that
   the ray could first touch the vortext tangent to it, in which case there
   would be two roots.  This is of no consequence.

-------------------------------------------------------------------------------------
*/

double find_t (
   double y_eye ,   // <= 0
   double z_eye ,   // > 1
   double row ,     // Row location in view window
   double col       // And column
   )
{
   int i;
   double eye_length, increment, lower, upper, lower_f, upper_f ;
   double x_delta, y_delta, z_delta ;

   eye_length = sqrt ( y_eye * y_eye + z_eye * z_eye ) ;
   increment = 0.001 / eye_length ;

   lower = 1.0 - (-y_eye / z_eye) / eye_length ; // Decent approximation
   lower *= 0.99 ;                               // Safety factor
   lower_f = find_diff ( y_eye , z_eye , row , col , lower , 1 , &x_delta , &y_delta , &z_delta ) ;

   if (lower_f > 1.e50)  // Looking above horizon?
      return 1.e60 ;

   assert ( lower_f > 0.0 ) ;

/*
   Loop until we go from positive to negative.
   The loop limit is for insurance only and should never be hit.
*/
 
   for (i=0 ; i<10000000 ; i++) {
      upper = lower + increment ;
      upper_f = find_diff ( y_eye , z_eye , row , col , upper , 0 , &x_delta , &y_delta , &z_delta ) ;
      if (upper_f <= 0.0)
         break ;
      lower = upper ;
      lower_f = upper_f ;
      }

   assert (i < 10000000) ;

   return  0.5 * (lower + upper) ; // It's not worth tangent interpolation
}


/*
-------------------------------------------------------------------------------------

   arc_length_from_t - Given the eye, the window pixel, and t, this computes
                       the total arc length from the bottom of the vortex to
                       the point on the vortex boundary defined by t.
                       It is assumed that the supplied t is such that the ray
                       projected from the eye through the window pixel does
                       truly lie on the vortex boundary. So this t must have
                       come from calling find_t().


-------------------------------------------------------------------------------------
*/

double arc_length_from_t (
   double y_eye ,       // <= 0
   double z_eye ,       // > 1
   double row ,         // Row location in view window
   double col ,         // And column
   double t ,           // Relative length of projected ray
   double *x_terminus , // Returned X coordinate of intersection
   double *y_terminus   // Returned Y coordinate of intersection
   )
{
   int k ;
   double theta, x_window, y_window, x_delta, y_delta, radius ;

/*
   The view window is the plane containing the object-world origin (0,0,0)
   that is perpendicular to the ray connecting the eye to (0,0,0).
   We are forcing x_eye to be 0, so it is easy to compute the angle theta
   of the eye above the XY plane.
   Project the vector from the eye to the window's (row, col) pixel
   and compute the object-world 3D coordinates of the intersection of this
   ray with the view window plane.  This points to what the eye sees
   for this view pixel.
*/

   theta = atan2 ( z_eye , -y_eye ) ; // ArcTan (z / -y) in range -Pi to Pi
   x_window = col ;
   y_window = row * sin ( theta ) ;
   
/*
   Compute the ray from the eye to this pixel on the view window.
   Recall that the X coordinate of the eye is 0.
*/

   x_delta = x_window ;
   y_delta = y_window - y_eye ;

/*
   Adjust this ray to have a length equal to t times the length of the delta ray.
   Again, recall that the X coordinate of the eye is 0.
*/

   *x_terminus = t * x_delta ;
   *y_terminus = t * y_delta + y_eye ;
   
   radius = sqrt ( *x_terminus * *x_terminus + *y_terminus * *y_terminus ) ;
   k = search_arc_table ( radius ) ;
   return arc_length[k] ;
}


/*
-------------------------------------------------------------------------------------

   compute_checker - Compute checkerboard raster image of vortex

-------------------------------------------------------------------------------------
*/

double compute_checker (
   double row ,     // Row location in view window
   double col       // And column
   )
{
   double t, arc_len ;
   double x_terminus, y_terminus, radius ;
   double tone_from_length, tone_from_angle, final_tone ;

   t = find_t ( vortex_params.eye_Y , vortex_params.eye_Z , row , col ) ;

   if (t > 1.e50)  // Above horizon?
      return 1.0 ;

   arc_len = arc_length_from_t ( vortex_params.eye_Y , vortex_params.eye_Z , row , col , t ,
                                 &x_terminus , &y_terminus ) ;

   tone_from_length = ((int) (vortex_params.length_divisions * arc_len)) % 2 ;


/*
   Compute the tone based on rotational position on the vortex
*/

   tone_from_angle = (atan2 ( x_terminus , y_terminus ) + PI) / (2.0 * PI) ;  // Scale it 0-1
   tone_from_angle = ((int) (vortex_params.angle_divisions * tone_from_angle + 0.5)) % 2 ;

   final_tone = ((int) tone_from_length + (int) (tone_from_angle)) % 2 ;  // XOR

   radius = x_terminus * x_terminus + y_terminus * y_terminus ;
   if (radius >= 1.0 + vortex_params.perim_width)
      final_tone = 1.0 ;
   else if (radius >= 1.0)
      final_tone = 0.0 ;

   return final_tone ;
}

/*
-------------------------------------------------------------------------------------

   compute_spiral - Compute spiral raster image of vortex

-------------------------------------------------------------------------------------
*/

double compute_spiral (
   double row ,     // Row location in view window
   double col       // And column
   )
{
   double t, arc_len ;
   double x_terminus, y_terminus, radius, final_tone ;
   double angle_part, length_part ;

   t = find_t ( vortex_params.eye_Y , vortex_params.eye_Z , row , col ) ;

   if (t > 1.e50)  // Above horizon?
      return 1.0 ;

   arc_len = arc_length_from_t ( vortex_params.eye_Y , vortex_params.eye_Z , row , col , t ,
                                 &x_terminus , &y_terminus ) ;

   length_part = vortex_params.length_divisions * arc_len ;


/*
   Compute the tone based on rotational position on the vortex
*/

   angle_part = (atan2 ( x_terminus , y_terminus ) + PI) / (2.0 * PI) ;  // Scale it 0-1
   angle_part = vortex_params.angle_divisions * angle_part ;

   final_tone = ((int) (length_part + angle_part)) % 2 ;

   radius = x_terminus * x_terminus + y_terminus * y_terminus ;
   if (radius >= 1.0 + vortex_params.perim_width)
      final_tone = 1.0 ;
   else if (radius >= 1.0)
      final_tone = 0.0 ;

   return final_tone ;
}