#include <stdlib.h>
#include <math.h>
#include "threed.h"
#include "grafprt.h"
#include "arrays.h"

#define PI 3.141592653589793238
#define BIGNUM 1.7e38

/*--------------------------------------------------------------------*/
/*  module threed.c in turbo c version 2.0, large code model.         */
/*                                                                    */
/*  released to the public domain by Gregory K. Landheim on January   */
/*  1, 1989, (c) all rights reserved.  anybody is authorized to use   */
/*  this code for any purpose whatsoever on the condition that they   */
/*  realize I assume absolutely no liability and give no warrantee    */
/*  for its use or performance.  it is distributed "as is", and if    */
/*  you use it you assume full liability for its use.                 */
/*                                                                    */
/*--------------------------------------------------------------------*/

/*------------------------function prototypes-------------------------*/

  void degrees_to_rads(float *horangle,float *elangle);
  int get_quadrant(float horangle);
  void transform_angles(float horangle,float elangle,int quadrant);
  void get_max_min(float xmin,float xmax,float ymin,float ymax,
                   float zmin,float zmax,
                   float *yminp,float *ymaxp,float *zminp,float *zmaxp);
  void axonometric(float x,float *y,float *z);
  void dosurf(float xmin,float xmax,float ymin,float ymax,
              int nx,int ny,float **z,
              int i0,int i1,int inci,
              int j0,int j1,int incj,int incmode,
              int fillcolor,int edgecolor);
  void drawfillquadrangle(int p[2][2],int q[2][2],
                          int fillcolor,int edgecolor);
  void filltriangle(int x1,int y1,int x2,int y2,
                    int x3,int y3,int fillcolor);
  void swapcoords(int *x1,int* y1,int* x2,int* y2);

  void transformbox(void);
  void drawboxbottom(int boxedgecolor);
  void drawboxback(int quadrant,int boxedgecolor);
  void drawboxtop(int boxedgecolor);
  void drawboxfront(int quadrant,int boxedgecolor);

/*-------------------------global declarations------------------------*/

int **Xa;                                              /* storage for */
int **Ya;                               /* filltriangle() edge lines  */

float Cosphi,Sinphi,Cospsi,Sinpsi,Cosel; /* constants for axonometric */
                                         /* projection                */

float Yb[5],Zb[5],Yt[5],Zt[5]; /* arrays for box bottoms and tops in  */
                               /* projected real coordinates          */
int Pb[5],Qb[5],Pt[5],Qt[5]; /* arrays for box bottoms and tops in    */
                             /* transformed integer coordinates       */

void far (*Linef)(int x1,int y1,int x2,int y2);       /* address of   */
                                 /* plot device line drawing function */
void far (*Csetf)(int color);/* address of plot device set color func */

/*--------------------------------------------------------------------*/
/*                                                                    */
/*  int surface(float xmin,float xmax,                                */
/*              float ymin,float ymax,                                */
/*              float zmin,float zmax,                                */
/*              int xminw,int xmaxw,int yminw,int ymaxw,              */
/*              int hmax,float **z,                                   */
/*              float horangle,float elangle,                         */
/*              int nx,int ny,int box,                                */
/*              int fillcolor,int edgecolor,int boxedgecol,           */
/*              void far csetfunc(int color),                         */
/*              void far linefunc(int x1,int y1,int x2,int y2))       */
/*                                                                    */
/*  this is a general-purpose surface plotting routine to plot a      */
/*  surface described as a matrix of gridded z-values in the x-y      */
/*  domain. the painter's algorithm is used for hidden line removal,  */
/*  and the 3-D surface coordinate system is assumed to be right-     */
/*  handed and cartesian.                                             */
/*                                                                    */
/*  xmin,xmax   minimum and maximum x-values of surface               */
/*  ymin,ymax   minimum and maximum y-values of surface               */
/*  zmin,zmax   minimum and maximum z-values of surface               */
/*  xminw,yminw lower left hand corner of plot window                 */
/*              in integer "world" coordinates                        */
/*  xmaxw,ymaxw upper right hand corner of plot window                */
/*              in integer "world" coordinates                        */
/*  hmax        maximum number of pixels in the plot device vertical  */
/*              direction                                             */
/*  z           address of first element of two dimensional array,    */
/*              (or its pointer equivalent)                           */
/*  horangle    horizontal view angle, counter-clockwise from 0,      */
/*              through 360 (or whatever) in decimal degrees.         */
/*              negative angles are also acceptable.                  */
/*  elangle     elevation view angle, positive upward from the x-y    */
/*              plane, in decimal degrees.  range from -90 through 90 */
/*  nx          number of grid points in the x-y plane in x direction */
/*  ny          number of grid points in the x-y plane in y direction */
/*  box         flag to indicate drawing of box.  if 0, don't draw.   */
/*  fillcolor   color of surface grid interior                        */
/*  edgecolor   color of surface grid lines                           */
/*  boxedgecol  color of box edge lines                               */
/*  csetfunc    address of line color setting function.               */
/*  linefunc    address of integer device coordinate line function.   */
/*                                                                    */
/*  external functions and variables:                                 */
/*                                                                    */
/*  alloc_2d_array() dynamically allocates a two-dimensional array    */
/*              using pointer notation.  prototype in arrays.h.       */
/*  free_2d_array() frees memory allocated by alloc_2d_array          */
/*  defreg()    defines the correspondence between integer world      */
/*              coordinates and real world coordinates for later use  */
/*              by transf() and itransf(). prototype in grafprt.h.    */
/*  transf()    transforms real world coordinates to integer world    */
/*              coordinates.  prototype in grafprt.h.                 */
/*  itransf()   transforms integer world coordinates to integer plot  */
/*              device coordinates.  prototype in grafprt.h.          */
/*  Maxheight   maximum pixel height of plot device.  declared in     */
/*              module grafprt.c.                                     */
/*                                                                    */
/*  a note about worlds:                                              */
/*  the real world is right-handed and rectangular and is expressed in*/
/*  floating point.  it is the projection onto two dimensions of the  */
/*  (x,y,z) coordinates of the three-dimensional surface.             */
/*  the integer world is the right-handed and rectangular integer     */
/*  representation of the maximum pixel plotting domain of the        */
/*  plot device.                                                      */
/*  the plot device is typically left-handed and rectangular.         */
/*  the transformation from real world to plot device coordinates is  */
/*  performed by two functions because of the design of the string    */
/*  drawing functions in module grafstr.                              */
/*  on a CRT display, drawing will normally be in portrait mode, and  */
/*  you probably will not need module grafstr.  if a CRT is the only  */
/*  device that will be used by your program, you may want to replace */
/*  functions transf() and itransf() with only a single function.     */
/*                                                                    */
/*--------------------------------------------------------------------*/

int surface(float xmin,float xmax,
            float ymin,float ymax,
            float zmin,float zmax,
            int xminw,int xmaxw,int yminw,int ymaxw,
            int hmax,float **z,
            float horangle,float elangle,
            int nx,int ny,int box,
            int fillcolor,int edgecolor,int boxedgecol,
            void far csetfunc(int color),
            void far linefunc(int x1,int y1,int x2,int y2))

{
   extern int Maxheight;                /* declared in module grafprt */

   float temp;
   int quadrant,bufsiz,incmode,invert;
   float ymaxp,yminp,zmaxp,zminp;/* max and mins of projected surface */

   if (abs(elangle) > 90.0) return(1);
   (elangle < 0.0) ? (invert = 1) : (invert = 0);

/* allocate memory for the arrays:                                    */

    bufsiz = ymaxw - yminw;
    if ((xmaxw - xminw) > bufsiz) bufsiz = xmaxw - xminw;

    Xa = (int **) alloc_2d_array(2,bufsiz,sizeof(int *),sizeof(int));
    if (!Xa) return(2);

    Ya = (int **) alloc_2d_array(2,bufsiz,sizeof(int *),sizeof(int));
    if (!Ya) return(3);

/* change the angles from degrees to radians:                         */

   degrees_to_rads(&horangle,&elangle);

/* determine the viewing quadrant:                                    */

   quadrant = get_quadrant(horangle);

/* determine the global projection angles:                            */

   transform_angles(horangle,elangle,quadrant);

/* swap xmin with xmax and ymin with ymax in quadrants 2 and 3:       */

   switch (quadrant)
     {
       case 1 : break;
       case 2 :
       case 3 : temp = xmin; xmin = xmax; xmax = temp;
                temp = ymin; ymin = ymax; ymax = temp;
                break;
       case 4 : break;
     }

/* find the maxima and minima of the real coordinate domain:          */

   get_max_min(xmin,xmax,ymin,ymax,zmin,zmax,
               &yminp,&ymaxp,&zminp,&zmaxp);

/* set maximum pixel height of device and define the plot region:     */

   Maxheight = hmax;
   defreg(yminp,ymaxp,zminp,zmaxp,xminw,xmaxw,yminw,ymaxw);

/* assign the color setting and line drawing function addresses to    */
/* global variables:                                                  */

   Csetf = csetfunc;
   Linef = linefunc;

/* put an optional box about the surface:                             */

   if (box)
     {
       transformbox();
       invert ? drawboxtop(boxedgecol) : drawboxbottom(boxedgecol);
       drawboxback(quadrant,boxedgecol);
     }

/* plot the surface as a function of quadrant.  make certain it       */
/* draws from back to front:                                          */

   switch (quadrant)
     {
       case 1 : (horangle > PI/4.0) ? (incmode = 0) : (incmode = 1);
                dosurf(xmin,xmax,ymin,ymax,nx,ny,z,
                  0,nx-1,1, 0,ny-1,1,incmode,fillcolor,edgecolor);
                break;
       case 2 : (horangle > 3.0*PI/4.0) ? (incmode = 1) : (incmode = 0);
                dosurf(xmin,xmax,ymin,ymax,nx,ny,z,
                  nx-1,0,-1,0,ny-1,1,incmode,fillcolor,edgecolor);
                break;
       case 3 : (horangle > 5.0*PI/4.0) ? (incmode = 0) : (incmode = 1);
                dosurf(xmin,xmax,ymin,ymax,nx,ny,z,
                nx-1,0,-1, ny-1,0,-1,incmode,fillcolor,edgecolor);
                break;
       case 4 : (horangle > 7.0*PI/4.0) ? (incmode = 1) : (incmode = 0);
                dosurf(xmin,xmax,ymin,ymax,nx,ny,z,
                 0,nx-1,1, ny-1,0,-1,incmode,fillcolor,edgecolor);
                break;
     }

   if (box)
     {
       invert ? drawboxbottom(boxedgecol) : drawboxtop(boxedgecol);
       drawboxfront(quadrant,boxedgecol);
     }

  free_2d_array((void **) Xa,2);
  free_2d_array((void **) Ya,2);

  return(0);

} /* int surface() */

/*--------------------------------------------------------------------*/
/*                                                                    */
/*  void degrees_to_rads(float * horangle,float *elangle)             */
/*                                                                    */
/*  converts input angles from decimal degrees to radians.            */
/*                                                                    */
/*--------------------------------------------------------------------*/

void degrees_to_rads(float * horangle,float *elangle)

{

/* make certain the angles are in the range 0 to 360:                 */

   do
     if (*horangle < 0.0) *horangle += 360.0;
   while (*horangle < 0.0);

   do
     if (*horangle > 360.0) *horangle -= 360.0;
   while (*horangle > 360.0);

   *elangle  = *elangle  * PI / 180.0;
   *horangle = *horangle * PI / 180.0;

}  /* void degrees_to_rads() */

/*--------------------------------------------------------------------*/
/*                                                                    */
/*  int get_quadrant(float horangle)                                  */
/*                                                                    */
/*  determines the quadrant in which the horizontal angle places the  */
/*  viewer.  the quadrants correspond to the standard math quadrants. */
/*                                                                    */
/*--------------------------------------------------------------------*/

int get_quadrant(float horangle)

{
  if ((horangle >= 0.0) && (horangle < PI/2.0))       return(1);
  else if ((horangle >= PI/2.0) && (horangle < PI))   return(2);
  else if ((horangle >= PI) && (horangle < 3*PI/2.0)) return(3);
  else                                                return(4);

} /* int get_quadrant() */

/*--------------------------------------------------------------------*/
/*                                                                    */
/*  void transform_angles(float horangle,float elangle,int quadrant)  */
/*                                                                    */
/*  computes the global transformation angles and their cosines and   */
/*  sines for later use in calls to axonometric().                    */
/*                                                                    */
/*--------------------------------------------------------------------*/

void transform_angles(float horangle,float elangle,int quadrant)

{
   float phi,psi,sinel;

/* avoid dividing by zero.                                            */
/* all this messing around with quadrants and whatnot is because the  */
/* angles that generate the cosines of 0.0 are a function of the      */
/* floating point library of the compiler and of round-off error.     */

   if (cos(horangle) == 0.0)            /* horizontal angle 90 or 270 */
     switch (quadrant)
            {
              case 1 :
              case 3 : phi = PI/2.0;
                       psi = 0.0;
                       break;
              case 2 :
              case 4 : phi = -PI/2.0;
                       psi =  PI;
                       break;
            }
   else
     {
       phi = atan(sin(horangle)/cos(horangle));
       psi = PI/2 - phi;
     }

   sinel  = sin(elangle);
   Sinphi = sinel * sin(phi);
   Cosphi = cos(phi);
   Sinpsi = sinel * sin(psi);
   Cospsi = cos(psi);
   Cosel  = cos(elangle);

} /* void transform_angles() */

/*--------------------------------------------------------------------*/
/*                                                                    */
/*  void get_max_min(float xmin,float xmax,float ymin,float ymax,     */
/*                   float zmin,float zmax,                           */
/*                   float *yminp,float *ymaxp,                       */
/*                   float *zminp,float *zmaxp);                      */
/*                                                                    */
/*  determines maxima and minima of projected real coordinates.       */
/*  stores axonometric projected box corner values in global          */
/*  arrays Yb,Zb,Yt,Zt.                                               */
/*                                                                    */
/*--------------------------------------------------------------------*/

void get_max_min(float xmin,float xmax,float ymin,float ymax,
                 float zmin,float zmax,
                 float *yminp,float *ymaxp,float *zminp,float *zmaxp)

{
  float ytemp,ztemp;
  float xb[5],yb[5];
  int i;

  *ymaxp = *zmaxp = -BIGNUM;
  *yminp = *zminp =  BIGNUM;

/* the maxima and minima of the real domain and range will be defined */
/* by the extremes of the function values:                            */

  xb[0] = xmin;        yb[0] = ymin;
  xb[1] = xmin;        yb[1] = ymax;
  xb[2] = xmax;        yb[2] = ymax;
  xb[3] = xmax;        yb[3] = ymin;
  xb[4] = xmin;        yb[4] = ymin;

/* find the maxima and minima of the real domain and range:           */

  for (i = 0; i < 5; i++)
    {
/*    do the bottom:                                                  */

      ytemp = yb[i];
      ztemp = zmin;

      axonometric(xb[i],&ytemp,&ztemp);

      if (ytemp > *ymaxp) *ymaxp = ytemp;
      if (ytemp < *yminp) *yminp = ytemp;
      if (ztemp > *zmaxp) *zmaxp = ztemp;
      if (ztemp < *zminp) *zminp = ztemp;

/*    store the extremes from the bottom for later box drawing:       */

      Yb[i] = ytemp;
      Zb[i] = ztemp;

/*    now do the top:                                                 */

      ytemp = yb[i];
      ztemp = zmax;

      axonometric(xb[i],&ytemp,&ztemp);

      if (ytemp > *ymaxp) *ymaxp = ytemp;
      if (ytemp < *yminp) *yminp = ytemp;
      if (ztemp > *zmaxp) *zmaxp = ztemp;
      if (ztemp < *zminp) *zminp = ztemp;

/*    store the extremes from the top for later box drawing:          */

      Yt[i] = ytemp;
      Zt[i] = ztemp;
    }
} /* void get_max_min() */

/*--------------------------------------------------------------------*/
/*                                                                    */
/*  void axonometric(float x,float *y,float *z)                       */
/*                                                                    */
/*  performs an axonometric transformation from three-dimensional     */
/*  real coordinates to two-dimensional real coordinates. result is   */
/*  returned in y,z.  globals Cosphi,Cospsi,Sinphi,Sinpsi & Cosel     */
/*  must be defined elsewhere.                                        */
/*                                                                    */
/*--------------------------------------------------------------------*/

void axonometric(float x,float *y,float *z)

{
   float ytemp;

   ytemp = *y;

   *y = ytemp * Cosphi - x * Cospsi;
   *z = -ytemp * Sinphi - x * Sinpsi + *z * Cosel;

}  /* void axonometric() */

/*--------------------------------------------------------------------*/
/*                                                                    */
/*  void dosurf(float xmin,float xmax,float ymin,float ymax,          */
/*              int nx,int ny,float **z,                              */
/*              int i0,int i1,int inci,                               */
/*              int j0,int j1,int incj,int incmode,                   */
/*              int fillcolor,int edgecolor)                          */
/*                                                                    */
/*  draws a three dimensional surface using the painter's algorithm.  */
/*                                                                    */
/*  the arrays are indexed according to the original gridded surface. */
/*  dosurf() draws the grid from back to front.                       */
/*                                                                    */
/*--------------------------------------------------------------------*/

void dosurf(float xmin,float xmax,float ymin,float ymax,
            int nx,int ny,float **z,
            int i0,int i1,int inci,
            int j0,int j1,int incj,int incmode,
            int fillcolor,int edgecolor)

{
  register int i,j;
  int k,l;
  int p[2][2],q[2][2];
  float x[2][2],y[2][2];
  float ytemp,ztemp,xdif,ydif,xdifp,ydifp;

  xdif = (xmax - xmin)/(nx-1);
  ydif = (ymax - ymin)/(ny-1);

  if (incmode)      /* draw and fill the columns of quadrangles first */
    {
      ydifp = incj * ydif;
      for (i=i0; i!=i1; i+=inci)
        {
          {                              /* draw the first quadrangle */
            for (k=0; k<2; k++)
              for (l=0; l<2; l++)
                {
                  x[k][l] = xmin + (i + inci * k) * xdif;
                  y[k][l] = ymin + (j0 + incj * l) * ydif;
                  ytemp = y[k][l];
                  ztemp = z[i+inci * k][j0+incj*l];
                  axonometric(x[k][l],&ytemp,&ztemp);
                  transf(ytemp,ztemp,&p[k][l],&q[k][l]);
                  itransf(&p[k][l],&q[k][l]);
                }
            drawfillquadrangle(p,q,fillcolor,edgecolor);

            y[0][0] = y[0][1];
            y[1][0] = y[1][1];
          }

        for(j=j0+incj; j!=j1; j+=incj)    /* do the other quadrangles */
          {
             {
               for (k=0;k<2;k++)
                 {
                   p[k][0] = p[k][1];
                   q[k][0] = q[k][1];
                   y[k][1] = y[k][0] + ydifp;
                   y[k][0] = y[k][1];
                   ztemp = z[i+inci*k][j+incj];
                   axonometric(x[k][1],&y[k][1],&ztemp);
                   transf(y[k][1],ztemp,&p[k][1],&q[k][1]);
                   itransf(&p[k][1],&q[k][1]);
                 }
             }
           drawfillquadrangle(p,q,fillcolor,edgecolor);
          }
        }
    }
  else                 /* draw and fill the rows of quadrangles first */
    {
      xdifp = inci * xdif;
      for (j=j0; j!=j1; j+=incj)
        {
          {                              /* draw the first quadrangle */
            for (k=0; k<2; k++)
              for (l=0; l<2; l++)
                {
                  x[k][l] = xmin + (i0 + inci * k) * xdif;
                  y[k][l] = ymin + (j + incj * l) * ydif;
                  ytemp = y[k][l];
                  ztemp = z[i0+inci * k][j+incj*l];
                  axonometric(x[k][l],&ytemp,&ztemp);
                  transf(ytemp,ztemp,&p[k][l],&q[k][l]);
                  itransf(&p[k][l],&q[k][l]);
                }
            drawfillquadrangle(p,q,fillcolor,edgecolor);

            x[0][0] = x[1][0];
            x[0][1] = x[1][1];
          }

        for (i=i0+inci; i!=i1; i+=inci)   /* do the other quadrangles */
          {
             {
               for (l=0;l<2;l++)
                 {
                   p[0][l] = p[1][l];
                   q[0][l] = q[1][l];
                   x[1][l] = x[0][l] + xdifp;
                   x[0][l] = x[1][l];
                   ytemp =  y[1][l];
                   ztemp = z[i+inci][j+l*incj];
                   axonometric(x[1][l],&ytemp,&ztemp);
                   transf(ytemp,ztemp,&p[1][l],&q[1][l]);
                   itransf(&p[1][l],&q[1][l]);
                 }
             }
           drawfillquadrangle(p,q,fillcolor,edgecolor);
          }
        }
    }
} /* void dosurf() */

/*--------------------------------------------------------------------*/
/*                                                                    */
/*  void drawfillquadrangle(int p[2][2],int q[2][2],                  */
/*                          int fillcolor,int edgecolor)              */
/*                                                                    */
/*  divides a quadrangle into two triangles.  calls filltriangle() to */
/*  blank the triangle, then draws the quadrangle perimeter.          */
/*                                                                    */
/*  array p contains the integer device x-coords of the rect. corners */
/*  array q contains the integer device y-coords of the rect. corners */
/*                                                                    */
/*--------------------------------------------------------------------*/

void drawfillquadrangle(int p[2][2],int q[2][2],
                        int filcol,int edgecolor)

{
   filltriangle(p[0][0],q[0][0],p[1][0],q[1][0],p[1][1],q[1][1],filcol);
   filltriangle(p[0][0],q[0][0],p[0][1],q[0][1],p[1][1],q[1][1],filcol);

   Csetf(edgecolor);
   Linef(p[0][0],q[0][0],p[0][1],q[0][1]);
   Linef(p[0][1],q[0][1],p[1][1],q[1][1]);
   Linef(p[1][1],q[1][1],p[1][0],q[1][0]);
   Linef(p[1][0],q[1][0],p[0][0],q[0][0]);

}  /* void drawfillquadrangle() */

/*--------------------------------------------------------------------*/
/*                                                                    */
/*  void filltriangle(int x0,int y0,int x1,int y1,                    */
/*                    int x2,int y2,int fillcolor)                    */
/*                                                                    */
/*  draws a triangle described by the three points passed.  the       */
/*  coordinates of the perimeter are computed and stored in two       */
/*  global two-dimensional arrays, Xa and Ya.  Xa[0][j] and Ya[0][j]  */
/*  store the two short sides of the triangle in contiguous order,    */
/*  while Xa[1][j] and Ya[1][j] store the long side of the triangle.  */
/*                                                                    */
/*  the triangle is filled by drawing horizontal lines between the    */
/*  two arrays.                                                       */
/*                                                                    */
/*  based on a translation to C by Greg Landheim of a pascal          */
/*  implementation of Bresenham's line algorithm by Prof Richard      */
/*  Rasala of Northeastern University Computer Science Dept.          */
/*  used with Prof. Rasala's permission.                              */
/*--------------------------------------------------------------------*/

void filltriangle(int x0,int y0,int x1,int y1,
                  int x2,int y2,int fillcolor)

{
    int  z,a,b,dx,dy,d,deltap,deltaq,jstart;
    int xl[2][3],yl[2][3];
    register int x,y,i,j,k;

/* sort the points in ascending vertical order:                       */

   if (y1 < y0) swapcoords(&x1,&y1,&x0,&y0);
   if (y2 < y0) swapcoords(&x2,&y2,&x0,&y0);
   if (y2 < y1) swapcoords(&x2,&y2,&x1,&y1);

/* stick them in arrays for the triangle edge computation:            */

   xl[0][0] = x0;              yl[0][0] = y0;
   xl[0][1] = x1;              yl[0][1] = y1;
   xl[0][2] = x2;              yl[0][2] = y2;

/* for the sake of using loops instead of if statements, we're        */
/* pretending that the long side of the triangle is made up of two    */
/* lines.  one of them has zero length: (x2,y2) to (x2,y2).           */

   xl[1][0] = x0;              yl[1][0] = y0;
   xl[1][1] = x2;              yl[1][1] = y2;
   xl[1][2] = x2;              yl[1][2] = y2;

/* i loops over the two triangles:                                    */
/* k loops over the two short sides and the one long side:            */

   for (i = 0; i < 2; i++)
     {
       jstart = 0;
       for (k = 1; k < 3; k++)
         {
           dx = abs(xl[i][k] - xl[i][k-1]);
           dy = abs(yl[i][k] - yl[i][k-1]);

           x = xl[i][k-1];
           y = yl[i][k-1];

           if (dy <= dx)
             {
               z = xl[i][k];
               a = (xl[i][k-1] <= xl[i][k]) ? 1 : -1;
               b = (yl[i][k-1] <= yl[i][k]) ? 1 : -1;

               deltap = dy << 1;
               d      = deltap - dx;
               deltaq = d - dx;
               j = jstart;
               Xa[i][j] = x;
               Ya[i][j] = y;

               while (x != z)
                 {
                   x += a;
                   (d < 0) ? (d += deltap) :
                             (y += b,d += deltaq,j += 1);
                   Xa[i][j] = x;
                   Ya[i][j] = y;
                 }
             }
           else
             {
               z = yl[i][k];
               a = (yl[i][k-1] <= yl[i][k]) ? 1 : -1;
               b = (xl[i][k-1] <= xl[i][k]) ? 1 : -1;

               deltap = dx << 1;
               d      = deltap - dy;
               deltaq = d - dy;
               j = jstart;
               Xa[i][j] = x;
               Ya[i][j] = y;

               while (y != z)
                 {
                   y += a;
                   j += 1;
                   (d < 0) ? (d += deltap) :
                             (x += b,d += deltaq);
                   Xa[i][j] = x;
                   Ya[i][j] = y;
                 }
             }
           jstart = dy;
         }
     }

/* draw the blank horizontal lines:                                   */

   Csetf(fillcolor);
   for (j=0; j < (yl[0][2] - yl[0][0]); j++)
     Linef(Xa[0][j],Ya[0][j],Xa[1][j],Ya[0][j]);

}  /* void filltriangle() */

/*--------------------------------------------------------------------*/
/*                                                                    */
/*  void swapcoords(int *x1,int* y1,int *x2,int *y2)                  */
/*                                                                    */
/*  swaps the integer plot device coordinates of two points.          */
/*                                                                    */
/*--------------------------------------------------------------------*/

void swapcoords(int *x1,int *y1,int *x2,int *y2)

{
   int temp;

   temp = *x1;
   *x1  = *x2;
   *x2  = temp;

   temp = *y1;
   *y1  = *y2;
   *y2  = temp;

}  /* void swapcoords() */

/*--------------------------------------------------------------------*/
/*                                                                    */
/*  box drawing routines.....                                         */
/*                                                                    */
/*  draw a box around the plotted surface describing the region       */
/*  enclosed by the projected xmin,ymin,xmax,ymax,zmin,zmax.          */
/*                                                                    */
/*  the projected real coordinates of the box corners were saved      */
/*  in arrays Yb,Zb,Yt,Zt during the call to get_max_min().           */
/*                                                                    */
/*  void transformbox()                                               */
/*  transforms 8 box corners to integer plot device coordinates.      */
/*                                                                    */
/*  void drawboxbottom(int color)                                     */
/*  draws the box bottom.                                             */
/*                                                                    */
/*  void drawboxback(int quadrant,int color)                          */
/*  draws the box back.                                               */
/*                                                                    */
/*  void drawboxtop(int color)                                        */
/*  draws the box top.                                                */
/*                                                                    */
/*  void drawboxfront(int quadrant,int color)                         */
/*  draws the box front.                                              */
/*                                                                    */
/*--------------------------------------------------------------------*/

void transformbox()

{
   int i;

   for (i = 0; i < 5; i++)
     {
       transf(Yb[i],Zb[i],&Pb[i],&Qb[i]);
       itransf(&Pb[i],&Qb[i]);
       transf(Yt[i],Zt[i],&Pt[i],&Qt[i]);
       itransf(&Pt[i],&Qt[i]);
     }
}  /* void transformbox() */

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

void drawboxbottom(int color)

{
   int i;

/* first, draw the bottom of the box:                                 */

   Csetf(color);
   for (i = 0; i < 4; i++)
     Linef(Pb[i],Qb[i],Pb[i+1],Qb[i+1]);

}  /* void drawboxbottom() */

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

void drawboxback(int quadrant,int color)

{
   int i,j;
   static int back[4][2] = { { 0,1 },
                             { 0,3 },
                             { 2,3 },
                             { 1,2 } };

   Csetf(color);
   for (i=0; i<2; i++)
     {
       j = back[quadrant-1][i];
       Linef(Pb[j],Qb[j],Pt[j],Qt[j]);
     }
} /* void drawboxback() */

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

void drawboxtop(int color)

{
   int i;

   Csetf(color);
   for (i = 0; i < 4; i++)
     Linef(Pt[i],Qt[i],Pt[i+1],Qt[i+1]);

}  /* void drawboxtop() */

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

void drawboxfront(int quadrant,int color)

{
   int i,j;
   static int front[4][2] = { { 2,3 },
                              { 1,2 },
                              { 0,1 },
                              { 0,3 } };

   Csetf(color);
   for (i=0; i<2; i++)
     {
       j = front[quadrant-1][i];
       Linef(Pb[j],Qb[j],Pt[j],Qt[j]);
     }
}  /* void drawboxfront() */

/*-------------------end of surface module threed.c-------------------*/


