#include <stdio.h>
#include <string.h>
#include <math.h>
#include <stdlib.h>
#include <iostream.h>
#include "titillat.h"
#include "plot3d.h"

#ifndef TRUE
#define TRUE -1
#endif
#ifndef FALSE
#define FALSE 0
#endif

typedef struct
          {
            int x;
            int y;
          } box_rec;

plot3d::plot3d()
  {
    prime_array_allocated=FALSE;
    plot_prepared=FALSE;
  }

plot3d::~plot3d()
  {
    if (prime_array_allocated)
      delete prime_array;
  }

int plot3d::prepare_plot(
  double (*f)(double,double),
  double x_min,
  double x_max,
  double y_min,
  double y_max,
  int    (*external_to_plot)(double,double),
  int    (*red)(double,double),
  int    x_division_count,
  int    y_division_count,
  double rotation_in_degrees,
  double tilt_in_degrees,
  double light_x,
  double light_y,
  double light_z)
//      This function prepares a plot for generation.  If returns TRUE if
// and only if it is successful.  If it is successful, "plot" may be called
// to actually generate the plot.  Its parameters are as follow:
//
//           f -- z=f(x,y), the function to be plotted.  Before the plot is
//      tilted or rotated, the z-axis runs from the bottom to the top of the
//      display, the y-axis runs from the left to the right of the display,
//      and the x-axis runs out of the display.
//
//           x_min -- the minimum value of x to be plotted.
//
//           x_max -- the maximum value of x to be plotted.
//
//           y_min -- the minimum value of y to be plotted.
//
//           y_max -- the maximum value of y to be plotted.
//
//           external_to_plot -- a function that returns TRUE if and only if a
//      point should be omitted from the plot.
//
//           red -- a function that returns TRUE if and only if a point should
//      be flagged for highlighting.  A point should be so flagged only if it
//      can be seen in the final plot.
//
//           x_division_count -- the number of x divisions to be used in
//      constructing the plot.  At least two must be specified.
//
//           y_division_count -- the number of y divisions to be used in
//      constructing the plot.  At least two must be specified.
//
//           rotation_in_degrees -- rotation (degrees) about an axis parallel to
//      the z-axis and through the center of the surface.
//
//           tilt_in_degrees -- tilt (degrees) about an axis through the center
//      of the surface and parallel to a line from the lower left hand corner of
//      the display to the lower right hand corner of the display.  The plot is
//      tilted after it is rotated.
//
//           (light_x,light_y,light_z) -- a vector pointing to the light source
//      (at infinity).  The light source remains fixed while the plot is rotated
//      or tilted.
    {
      int       display_ready;
      prime_rec initialized_prime_rec;
      int       result;
      int       x_byte_num;

      plot_prepared=FALSE;
      if (prime_array_allocated)
        {
          delete prime_array;
          prime_array_allocated=FALSE;
        }
      num_x_divisions=x_division_count;
      num_y_divisions=y_division_count;
      num_primes=(long) num_x_divisions;
      num_primes*=((long) num_y_divisions);
        // Number of quadrilaterals composing the plot.
      initialized_prime_rec.base_z=(unsigned char) '\0';
      initialized_prime_rec.color=(unsigned char) '\0';
      initialized_prime_rec.x=(float) 0.0;
      initialized_prime_rec.y=(float) 0.0;
      initialized_prime_rec.z=(float) 0.0;
      initialized_prime_rec.x_division_index=0;
      initialized_prime_rec.y_division_index=0;
      prime_array=new varray<prime_rec>(initialized_prime_rec,num_primes,5);
        // Virtual array of quadrilaterals composing the plot.
      if (prime_array_allocated=(prime_array->allocated()))
        {
          titillator_ptr=new titillator;
          rotation=rotation_in_degrees;
          tilt=tilt_in_degrees;
          light.x=light_x;
          light.y=light_y;
          light.z=light_z;
          evaluate_and_transform(f,x_min,x_max,y_min,y_max,
           num_x_divisions,num_y_divisions,rotation,tilt,
           external_to_plot,red);
            // Compute the vertices, etc. of the quadrilaterals composing the
            // plot.
          shade();
            // Compute the shade of gray for each quadrilateral composing the
            // plot.
          adjust_perspective();
            // Force parallel lines running away from the viewer to converge at
            // the horizon.
          sort_back_to_front();
          plot_prepared=TRUE;
          delete titillator_ptr;
          display_ready=display_initialized();
            // Do whatever is necessary to prepare the display.
        }
      return (prime_array_allocated & display_ready);
    }

void plot3d::evaluate_and_transform(
  double (*f)(double,double),
  double x_min,
  double x_max,
  double y_min,
  double y_max,
  int    num_x_divisions,
  int    num_y_divisions,
  double rotation,
  double tilt,
  int    (*external_to_plot)(double,double),
  int    (*red)(double,double))
// Compute the vertices, etc. for each quadrilateral composing the plot.
    {
      double    cos_rotation;
      double    cos_tilt;
      double    degrees_per_radian;
      double    magnitude;
      prime_rec *prime;
      long      prime_num;
      double    radians;
      double    sin_rotation;
      double    sin_tilt;
      double    tem_x;
      double    tem_y;
      double    tem_z;
      double    x;
      double    x_delta;
      int       x_division_num;
      double    x_rotated;
      double    y;
      double    y_delta;
      int       y_division_num;
      double    z;

      degrees_per_radian=45.0/atan(1.0);
      radians=tilt/degrees_per_radian;
      cos_tilt=cos(radians);
      sin_tilt=sin(radians);
      radians=rotation/degrees_per_radian;
      cos_rotation=cos(radians);
      sin_rotation=sin(radians);
      z=f(x_min,y_min);
      x_rotated=x_min*cos_rotation+y_min*sin_rotation;
      y_prime_min=-x_min*sin_rotation+y_min*cos_rotation;
      z_prime_min=-x_rotated*sin_tilt+z*cos_tilt;
      y_prime_max=y_prime_min;
      z_prime_max=z_prime_min;
      x_prime_max=x_rotated*cos_tilt+z*sin_tilt;
      x_delta=(double) (num_x_divisions-1);
      x_delta=(x_max-x_min)/x_delta;
      y_delta=(double) (num_y_divisions-1);
      y_delta=(y_max-y_min)/y_delta;
      x=x_min;
      prime_num=(long) 0;
      for (x_division_num=0; x_division_num < num_x_divisions; x_division_num++)
        {
          titillator_ptr->titillate();
          y=y_min;
          for (y_division_num=0; y_division_num < num_y_divisions;
           y_division_num++)
            {
              z=f(x,y);
              prime=prime_array->vm_ptr(prime_num);
              if (external_to_plot(x,y))
                prime->base_z=(unsigned char) 3;
              else
                if (red(x,y))
                  prime->base_z=(unsigned char) 2;
                else
                  prime->base_z=(unsigned char) 1;
              prime->x_division_index=x_division_num;
              prime->y_division_index=y_division_num;
              x_rotated=x*cos_rotation+y*sin_rotation;
              tem_y=(-x*sin_rotation+y*cos_rotation);
              prime->y=(float) tem_y;
              tem_x=(x_rotated*cos_tilt+z*sin_tilt);
              prime->x=(float) tem_x;
              tem_z=(-x_rotated*sin_tilt+z*cos_tilt);
              prime->z=(float) tem_z;
              if (tem_x > x_prime_max)
                x_prime_max=tem_x;
              if (tem_y < y_prime_min)
                y_prime_min=tem_y;
              if (tem_y > y_prime_max)
                y_prime_max=tem_y;
              if (tem_z < z_prime_min)
                z_prime_min=tem_z;
              if (tem_z > z_prime_max)
                z_prime_max=tem_z;
              y+=y_delta;
              prime_num++;
            }
          x+=x_delta;
        }
      magnitude=light.x*light.x+light.y*light.y+light.z*light.z;
      magnitude=sqrt(magnitude);
      light.x=light.x/magnitude;
      light.y=light.y/magnitude;
      light.z=light.z/magnitude;
      return;
    }

void plot3d::shade()
// Compute the shade of gray for each quadrilateral composing the plot.
    {
      double     magnitude;
      vertex_rec normal;
      prime_rec  *prime;
      long       prime_num;
      vertex_rec vertex [4];
      int        x_division_num;
      int        y_division_num;

      color_min=(unsigned char) (NUM_COLORS-1);
      color_max=(unsigned char) '\0';
      prime_num=num_primes;
      for (x_division_num=num_x_divisions-1; x_division_num >= 0;
       --x_division_num)
        {
          titillator_ptr->titillate();
          for (y_division_num=num_y_divisions-1; y_division_num >= 0;
           --y_division_num)
            {
              prime_num--;
              prime=prime_array->vm_ptr(prime_num);
              vertex[0].x=(double) (prime->x);
              vertex[0].y=(double) (prime->y);
              vertex[0].z=(double) (prime->z);
              if (x_division_num < (num_x_divisions-1))
                if (y_division_num < (num_y_divisions-1))
                  {
                    prime_num+=((long) num_y_divisions);
                    prime=prime_array->vm_ptr(prime_num);
                    vertex[1].x=(double) (prime->x);
                    vertex[1].y=(double) (prime->y);
                    vertex[1].z=(double) (prime->z);
                    prime_num++;
                    prime=prime_array->vm_ptr(prime_num);
                    vertex[2].x=(double) (prime->x);
                    vertex[2].y=(double) (prime->y);
                    vertex[2].z=(double) (prime->z);
                    prime_num-=((long) num_y_divisions);
                    prime=prime_array->vm_ptr(prime_num);
                    vertex[3].x=(double) (prime->x);
                    vertex[3].y=(double) (prime->y);
                    vertex[3].z=(double) (prime->z);
                    prime_num--;
                  }
                else
                  {
                    prime_num--;
                    prime=prime_array->vm_ptr(prime_num);
                    vertex[1].x=(double) (prime->x);
                    vertex[1].y=(double) (prime->y);
                    vertex[1].z=(double) (prime->z);
                    prime_num+=((long) num_y_divisions);
                    prime=prime_array->vm_ptr(prime_num);
                    vertex[2].x=(double) (prime->x);
                    vertex[2].y=(double) (prime->y);
                    vertex[2].z=(double) (prime->z);
                    prime_num++;
                    prime=prime_array->vm_ptr(prime_num);
                    vertex[3].x=(double) (prime->x);
                    vertex[3].y=(double) (prime->y);
                    vertex[3].z=(double) (prime->z);
                    prime_num-=((long) num_y_divisions);
                  }
              else
                if (y_division_num < (num_y_divisions-1))
                  {
                    prime_num++;
                    prime=prime_array->vm_ptr(prime_num);
                    vertex[1].x=(double) (prime->x);
                    vertex[1].y=(double) (prime->y);
                    vertex[1].z=(double) (prime->z);
                    prime_num-=((long) num_y_divisions);
                    prime=prime_array->vm_ptr(prime_num);
                    vertex[2].x=(double) (prime->x);
                    vertex[2].y=(double) (prime->y);
                    vertex[2].z=(double) (prime->z);
                    prime_num--;
                    prime=prime_array->vm_ptr(prime_num);
                    vertex[3].x=(double) (prime->x);
                    vertex[3].y=(double) (prime->y);
                    vertex[3].z=(double) (prime->z);
                    prime_num+=((long) num_y_divisions);
                  }
                else
                  {
                    prime_num-=((long) num_y_divisions);
                    prime=prime_array->vm_ptr(prime_num);
                    vertex[1].x=(double) (prime->x);
                    vertex[1].y=(double) (prime->y);
                    vertex[1].z=(double) (prime->z);
                    prime_num--;
                    prime=prime_array->vm_ptr(prime_num);
                    vertex[2].x=(double) (prime->x);
                    vertex[2].y=(double) (prime->y);
                    vertex[2].z=(double) (prime->z);
                    prime_num+=((long) num_y_divisions);
                    prime=prime_array->vm_ptr(prime_num);
                    vertex[3].x=(double) (prime->x);
                    vertex[3].y=(double) (prime->y);
                    vertex[3].z=(double) (prime->z);
                    prime_num++;
                  }
              // Compute the normal to a quadrilateral by averaging the
              // normals to each vertex of the quadrilateral.
              normal.x
               =(vertex[1].y-vertex[0].y)*(vertex[3].z-vertex[0].z)
               -(vertex[3].y-vertex[0].y)*(vertex[1].z-vertex[0].z)
               +(vertex[2].y-vertex[1].y)*(vertex[0].z-vertex[1].z)
               -(vertex[0].y-vertex[1].y)*(vertex[2].z-vertex[1].z)
               +(vertex[3].y-vertex[2].y)*(vertex[1].z-vertex[2].z)
               -(vertex[1].y-vertex[2].y)*(vertex[3].z-vertex[2].z)
               +(vertex[0].y-vertex[3].y)*(vertex[2].z-vertex[3].z)
               -(vertex[2].y-vertex[3].y)*(vertex[0].z-vertex[3].z);
              normal.y
               =(vertex[3].x-vertex[0].x)*(vertex[1].z-vertex[0].z)
               -(vertex[1].x-vertex[0].x)*(vertex[3].z-vertex[0].z)
               +(vertex[0].x-vertex[1].x)*(vertex[2].z-vertex[1].z)
               -(vertex[2].x-vertex[1].x)*(vertex[0].z-vertex[1].z)
               +(vertex[1].x-vertex[2].x)*(vertex[3].z-vertex[2].z)
               -(vertex[3].x-vertex[2].x)*(vertex[1].z-vertex[2].z)
               +(vertex[2].x-vertex[3].x)*(vertex[0].z-vertex[3].z)
               -(vertex[0].x-vertex[3].x)*(vertex[2].z-vertex[3].z);
              normal.z
               =(vertex[1].x-vertex[0].x)*(vertex[3].y-vertex[0].y)
               -(vertex[3].x-vertex[0].x)*(vertex[1].y-vertex[0].y)
               +(vertex[2].x-vertex[1].x)*(vertex[0].y-vertex[1].y)
               -(vertex[0].x-vertex[1].x)*(vertex[2].y-vertex[1].y)
               +(vertex[3].x-vertex[2].x)*(vertex[1].y-vertex[2].y)
               -(vertex[1].x-vertex[2].x)*(vertex[3].y-vertex[2].y)
               +(vertex[0].x-vertex[3].x)*(vertex[2].y-vertex[3].y)
               -(vertex[2].x-vertex[3].x)*(vertex[0].y-vertex[3].y);
              prime=prime_array->vm_ptr(prime_num);
              magnitude
               =sqrt(normal.x*normal.x+normal.y*normal.y+normal.z*normal.z);
              if (magnitude == 0.0)
                {
                  color_min=(unsigned char) '\0';
                  prime->color=(unsigned char) '\0';
                }
              else
                {
                  prime->color
                   =(unsigned char) int((double(NUM_COLORS)/2.0)
                   *(1.0+(light.x*normal.x+light.y*normal.y
                   +light.z*normal.z)/magnitude)); // shadows not absolute
                  if (prime->color
                   >= (unsigned char) (NUM_COLORS))
                    prime->color=(unsigned char) (NUM_COLORS-1);
                  if ((prime->color) < color_min)
                    color_min=(prime->color);
                  if ((prime->color) > color_max)
                    color_max=(prime->color);
                }
            }
        }
      return;
    }

void plot3d::adjust_perspective()
// Make parallel lines running away from the viewer appear to converge at the
// horizon. 
    {
      prime_rec  *prime;
      long       prime_num;
      double     tem;
      vertex_rec vertex [4];
      int        x_division_num;
      double     x_eye;
      double     y_center;
      int        y_division_num;
      double     z_center;

      if ((y_prime_max-y_prime_min) > (z_prime_max-z_prime_min))
        x_eye=1.1*(y_prime_max-y_prime_min)+x_prime_max;
      else
        x_eye=1.1*(z_prime_max-z_prime_min)+x_prime_max;
      if (((y_prime_max-y_prime_min) > (z_prime_max-z_prime_min))
      ||  (z_prime_max != z_prime_min))
        {
          y_center=(y_prime_max+y_prime_min)/2.0;
          z_center=(z_prime_max+z_prime_min)/2.0;
          prime_num=(long) 0;
          for (x_division_num=0; x_division_num < num_x_divisions;
           x_division_num++)
            {
              titillator_ptr->titillate();
              for (y_division_num=0; y_division_num < num_y_divisions;
               y_division_num++)
                {
                  prime=prime_array->vm_ptr(prime_num);
                  vertex[0].x=(double) (prime->x);
                  vertex[0].y=(double) (prime->y);
                  vertex[0].z=(double) (prime->z);
                  if (x_division_num < (num_x_divisions-1))
                    if (y_division_num < (num_y_divisions-1))
                      {
                        prime_num+=((long) num_y_divisions);
                        prime=prime_array->vm_ptr(prime_num);
                        vertex[1].x=(double) (prime->x);
                        vertex[1].y=(double) (prime->y);
                        vertex[1].z=(double) (prime->z);
                        prime_num++;
                        prime=prime_array->vm_ptr(prime_num);
                        vertex[2].x=(double) (prime->x);
                        vertex[2].y=(double) (prime->y);
                        vertex[2].z=(double) (prime->z);
                        prime_num-=((long) num_y_divisions);
                        prime=prime_array->vm_ptr(prime_num);
                        vertex[3].x=(double) (prime->x);
                        vertex[3].y=(double) (prime->y);
                        vertex[3].z=(double) (prime->z);
                        prime_num--;
                      }
                    else
                      {
                        prime_num--;
                        prime=prime_array->vm_ptr(prime_num);
                        vertex[1].x=(double) (prime->x);
                        vertex[1].y=(double) (prime->y);
                        vertex[1].z=(double) (prime->z);
                        prime_num+=((long) num_y_divisions);
                        prime=prime_array->vm_ptr(prime_num);
                        vertex[2].x=(double) (prime->x);
                        vertex[2].y=(double) (prime->y);
                        vertex[2].z=(double) (prime->z);
                        prime_num++;
                        prime=prime_array->vm_ptr(prime_num);
                        vertex[3].x=(double) (prime->x);
                        vertex[3].y=(double) (prime->y);
                        vertex[3].z=(double) (prime->z);
                        prime_num-=((long) num_y_divisions);
                      }
                  else
                    if (y_division_num < (num_y_divisions-1))
                      {
                        prime_num++;
                        prime=prime_array->vm_ptr(prime_num);
                        vertex[1].x=(double) (prime->x);
                        vertex[1].y=(double) (prime->y);
                        vertex[1].z=(double) (prime->z);
                        prime_num-=((long) num_y_divisions);
                        prime=prime_array->vm_ptr(prime_num);
                        vertex[2].x=(double) (prime->x);
                        vertex[2].y=(double) (prime->y);
                        vertex[2].z=(double) (prime->z);
                        prime_num--;
                        prime=prime_array->vm_ptr(prime_num);
                        vertex[3].x=(double) (prime->x);
                        vertex[3].y=(double) (prime->y);
                        vertex[3].z=(double) (prime->z);
                        prime_num+=((long) num_y_divisions);
                      }
                    else
                      {
                        prime_num-=((long) num_y_divisions);
                        prime=prime_array->vm_ptr(prime_num);
                        vertex[1].x=(double) (prime->x);
                        vertex[1].y=(double) (prime->y);
                        vertex[1].z=(double) (prime->z);
                        prime_num--;
                        prime=prime_array->vm_ptr(prime_num&