xref: /dports/math/plplot-ada/plplot-5.15.0/src/plot3d.c

//! @file
//!
//! 3d plot routines.
//!
// Copyright (C) 2004-2014 Alan W. Irwin
// Copyright (C) 2004  Joao Cardoso
// Copyright (C) 2004  Andrew Ross
//
// This file is part of PLplot.
//
// PLplot is free software; you can redistribute it and/or modify
// it under the terms of the GNU Library General Public License as published
// by the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// PLplot is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU Library General Public License for more details.
//
// You should have received a copy of the GNU Library General Public License
// along with PLplot; if not, write to the Free Software
// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
//

#include "plplotP.h"

// Internal constants

#define  BINC    50             // Block size for memory allocation

static PLINT pl3mode = 0;       // 0 3d solid; 1 mesh plot
static PLINT pl3upv  = 1;       // 1 update view; 0 no update

static PLINT zbflg = 0, zbcol;
static PLFLT zbtck, zbwidth;

static PLINT *oldhiview = NULL;
static PLINT *oldloview = NULL;
static PLINT *newhiview = NULL;
static PLINT *newloview = NULL;
static PLINT *utmp      = NULL;
static PLINT *vtmp      = NULL;
static PLFLT *ctmp      = NULL;

static PLINT mhi, xxhi, newhisize;
static PLINT mlo, xxlo, newlosize;

// Light source for shading
static PLFLT xlight, ylight, zlight;
static PLINT falsecolor = 0;
static PLFLT fc_minz, fc_maxz;

// Prototypes for static functions

static void plgrid3( PLFLT );
static void plnxtv( PLINT *, PLINT *, PLFLT*, PLINT, PLINT );
static void
plside3( PLFLT_VECTOR x, PLFLT_VECTOR y, PLF2OPS zops, PLPointer zp, PLINT nx, PLINT ny, PLINT opt );
static void
plt3zz( PLINT x0, PLINT y0, PLINT dx, PLINT dy, PLINT flag, PLINT *init,
        PLFLT_VECTOR x, PLFLT_VECTOR y, PLF2OPS zops, PLPointer zp, PLINT nx, PLINT ny,
        PLINT *u, PLINT *v, PLFLT* c );
static void plnxtvhi( PLINT *, PLINT *, PLFLT*, PLINT, PLINT );
static void plnxtvlo( PLINT *, PLINT *, PLFLT*, PLINT, PLINT );
static void plnxtvhi_draw( PLINT *u, PLINT *v, PLFLT* c, PLINT n );

static void savehipoint( PLINT, PLINT );
static void savelopoint( PLINT, PLINT );
static void swaphiview( void );
static void swaploview( void );
static void myexit( PLCHAR_VECTOR );
static void myabort( PLCHAR_VECTOR );
static void freework( void );
static int  plabv( PLINT, PLINT, PLINT, PLINT, PLINT, PLINT );
static void pl3cut( PLINT, PLINT, PLINT, PLINT, PLINT,
                    PLINT, PLINT, PLINT, PLINT *, PLINT * );
static PLFLT plGetAngleToLight( PLFLT* x, PLFLT* y, PLFLT* z );
static void plP_draw3d( PLINT x, PLINT y, PLFLT *c, PLINT j, PLINT move );
//static void plxyindexlimits( PLINT instart, PLINT inn,
//                             PLINT *inarray_min, PLINT *inarray_max,
//                             PLINT *outstart, PLINT *outn, PLINT outnmax,
//                             PLINT *outarray_min, PLINT *outarray_max );


// #define MJL_HACK 1
#if MJL_HACK
static void plP_fill3( PLINT x0, PLINT y0, PLINT x1, PLINT y1,
                       PLINT x2, PLINT y2, PLINT j );
static void plP_fill4( PLINT x0, PLINT y0, PLINT x1, PLINT y1,
                       PLINT x2, PLINT y2, PLINT x3, PLINT y3, PLINT j );
#endif

//--------------------------------------------------------------------------
// void plsetlightsource(x, y, z)
//
// Sets the position of the light source.
//--------------------------------------------------------------------------

void
c_pllightsource( PLFLT x, PLFLT y, PLFLT z )
{
    xlight = x;
    ylight = y;
    zlight = z;
}

//--------------------------------------------------------------------------
// void plmesh(x, y, z, nx, ny, opt)
//
// Plots a mesh representation of the function z[x][y]. The x values
// are stored as x[0..nx-1], the y values as y[0..ny-1], and the
// z values are in the 2-d array z[][]. The integer "opt" specifies:
// see plmeshc() below.
//--------------------------------------------------------------------------

void
c_plmesh( PLFLT_VECTOR x, PLFLT_VECTOR y, PLFLT_MATRIX z, PLINT nx, PLINT ny, PLINT opt )
{
    plfplot3dc( x, y, plf2ops_c(), (PLPointer) z, nx, ny, opt | MESH, NULL, 0 );
}

void
plfmesh( PLFLT_VECTOR x, PLFLT_VECTOR y, PLF2OPS zops, PLPointer zp,
         PLINT nx, PLINT ny, PLINT opt )
{
    plfplot3dc( x, y, zops, zp, nx, ny, opt | MESH, NULL, 0 );
}

//--------------------------------------------------------------------------
// void plmeshc(x, y, z, nx, ny, opt, clevel, nlevel)
//
// Plots a mesh representation of the function z[x][y]. The x values
// are stored as x[0..nx-1], the y values as y[0..ny-1], and the
// z values are in the 2-d array z[][]. The integer "opt" specifies:
//
// DRAW_LINEX   draw lines parallel to the X axis
// DRAW_LINEY   draw lines parallel to the Y axis
// DRAW_LINEXY  draw lines parallel to both the X and Y axis
// MAG_COLOR    draw the mesh with a color dependent of the magnitude
// BASE_CONT    draw contour plot at bottom xy plane
// TOP_CONT     draw contour plot at top xy plane (not yet)
// DRAW_SIDES   draw sides
//
// or any bitwise combination, e.g. "MAG_COLOR | DRAW_LINEX"
//
//--------------------------------------------------------------------------

void
c_plmeshc( PLFLT_VECTOR x, PLFLT_VECTOR y, PLFLT_MATRIX z, PLINT nx, PLINT ny, PLINT opt,
           PLFLT_VECTOR clevel, PLINT nlevel )
{
    plfplot3dc( x, y, plf2ops_c(), (PLPointer) z, nx, ny, opt | MESH, clevel, nlevel );
}

void
plfmeshc( PLFLT_VECTOR x, PLFLT_VECTOR y, PLF2OPS zops, PLPointer zp,
          PLINT nx, PLINT ny, PLINT opt, PLFLT_VECTOR clevel, PLINT nlevel )
{
    plfplot3dc( x, y, zops, zp, nx, ny, opt | MESH, clevel, nlevel );
}

// clipping helper for 3d polygons

int
plP_clip_poly( int Ni, PLFLT *Vi[3], int axis, PLFLT dir, PLFLT offset )
{
    int   anyout = 0;
    PLFLT _in[PL_MAXPOLY], _T[3][PL_MAXPOLY];
    PLFLT *in, *T[3], *TT = NULL;
    int   No = 0;
    int   i, j, k;

    if ( Ni > PL_MAXPOLY )
    {
        in = (PLFLT *) malloc( sizeof ( PLFLT ) * (size_t) Ni );
        TT = (PLFLT *) malloc( 3 * sizeof ( PLFLT ) * (size_t) Ni );

        if ( in == NULL || TT == NULL )
        {
            plexit( "plP_clip_poly: insufficient memory for large polygon" );
        }

        T[0] = &TT[0];
        T[1] = &TT[Ni];
        T[2] = &TT[2 * Ni];
    }
    else
    {
        in   = _in;
        T[0] = &_T[0][0];
        T[1] = &_T[1][0];
        T[2] = &_T[2][0];
    }

    for ( i = 0; i < Ni; i++ )
    {
        in[i]   = Vi[axis][i] * dir + offset;
        anyout += in[i] < 0;
    }

    // none out
    if ( anyout == 0 )
        return Ni;

    // all out
    if ( anyout == Ni )
    {
        return 0;
    }

    // some out
    // copy over to a temp vector
    for ( i = 0; i < 3; i++ )
    {
        for ( j = 0; j < Ni; j++ )
        {
            T[i][j] = Vi[i][j];
        }
    }
    // copy back selectively
    for ( i = 0; i < Ni; i++ )
    {
        j = ( i + 1 ) % Ni;

        if ( in[i] >= 0 && in[j] >= 0 )
        {
            for ( k = 0; k < 3; k++ )
                Vi[k][No] = T[k][j];
            No++;
        }
        else if ( in[i] >= 0 && in[j] < 0 )
        {
            PLFLT u = in[i] / ( in[i] - in[j] );
            for ( k = 0; k < 3; k++ )
                Vi[k][No] = T[k][i] * ( 1 - u ) + T[k][j] * u;
            No++;
        }
        else if ( in[i] < 0 && in[j] >= 0 )
        {
            PLFLT u = in[i] / ( in[i] - in[j] );
            for ( k = 0; k < 3; k++ )
                Vi[k][No] = T[k][i] * ( 1 - u ) + T[k][j] * u;
            No++;
            for ( k = 0; k < 3; k++ )
                Vi[k][No] = T[k][j];
            No++;
        }
    }

    if ( Ni > PL_MAXPOLY )
    {
        free( in );
        free( TT );
    }

    return No;
}

// helper for plsurf3d, similar to c_plfill3()
static void
shade_triangle( PLFLT x0, PLFLT y0, PLFLT z0,
                PLFLT x1, PLFLT y1, PLFLT z1,
                PLFLT x2, PLFLT y2, PLFLT z2 )
{
    int   i;
    // arrays for interface to core functions
    short u[6], v[6];
    PLFLT x[6], y[6], z[6];
    int   n;
    PLFLT xmin, xmax, ymin, ymax, zmin, zmax, zscale;
    PLFLT *V[3];

    plP_gdom( &xmin, &xmax, &ymin, &ymax );
    plP_grange( &zscale, &zmin, &zmax );

    x[0] = x0; x[1] = x1; x[2] = x2;
    y[0] = y0; y[1] = y1; y[2] = y2;
    z[0] = z0; z[1] = z1; z[2] = z2;
    n    = 3;

    V[0] = x; V[1] = y; V[2] = z;

    n = plP_clip_poly( n, V, 0, 1, -xmin );
    n = plP_clip_poly( n, V, 0, -1, xmax );
    n = plP_clip_poly( n, V, 1, 1, -ymin );
    n = plP_clip_poly( n, V, 1, -1, ymax );
    n = plP_clip_poly( n, V, 2, 1, -zmin );
    n = plP_clip_poly( n, V, 2, -1, zmax );

    if ( n > 0 )
    {
        if ( falsecolor )
            plcol1( ( ( z[0] + z[1] + z[2] ) / 3. - fc_minz ) / ( fc_maxz - fc_minz ) );
        else
            plcol1( plGetAngleToLight( x, y, z ) );

        for ( i = 0; i < n; i++ )
        {
            u[i] = (short) plP_wcpcx( plP_w3wcx( x[i], y[i], z[i] ) );
            v[i] = (short) plP_wcpcy( plP_w3wcy( x[i], y[i], z[i] ) );
        }
        u[n] = u[0];
        v[n] = v[0];

#ifdef SHADE_DEBUG // show triangles
        plcol0( 1 );
        x[3] = x[0]; y[3] = y[0]; z[3] = z[0];
        plline3( 4, x, y, z );
#else   // fill triangles
        plP_fill( u, v, n + 1 );
#endif
    }
}

//--------------------------------------------------------------------------
// void plsurf3d(x, y, z, nx, ny, opt, clevel, nlevel)
//
// Plots the 3-d surface representation of the function z[x][y].
// The x values are stored as x[0..nx-1], the y values as y[0..ny-1],
//  and the z values are in the 2-d array z[][].  The integer "opt" specifies:
// see plsurf3dl() below.
//--------------------------------------------------------------------------

void
c_plsurf3d( PLFLT_VECTOR x, PLFLT_VECTOR y, PLFLT_MATRIX z, PLINT nx, PLINT ny,
            PLINT opt, PLFLT_VECTOR clevel, PLINT nlevel )
{
    plfsurf3d( x, y, plf2ops_c(), (PLPointer) z, nx, ny,
        opt, clevel, nlevel );
}

void
plfsurf3d( PLFLT_VECTOR x, PLFLT_VECTOR y, PLF2OPS zops, PLPointer zp,
           PLINT nx, PLINT ny, PLINT opt, PLFLT_VECTOR clevel, PLINT nlevel )
{
    PLINT i;
    PLINT *indexymin = (PLINT *) malloc( (size_t) nx * sizeof ( PLINT ) );
    PLINT *indexymax = (PLINT *) malloc( (size_t) nx * sizeof ( PLINT ) );

    if ( !indexymin || !indexymax )
        plexit( "plsurf3d: Out of memory." );
    for ( i = 0; i < nx; i++ )
    {
        indexymin[i] = 0;
        indexymax[i] = ny;
    }
    plfsurf3dl( x, y, zops, zp, nx, ny, opt, clevel, nlevel,
        0, nx, indexymin, indexymax );
    free_mem( indexymin );
    free_mem( indexymax );
}

//--------------------------------------------------------------------------
// void plsurf3dl(x, y, z, nx, ny, opt, clevel, nlevel,
// indexxmin, indexxmax, indexymin, indexymax)
//
// Plots the 3-d surface representation of the function z[x][y].
// The x values are stored as x[0..nx-1], the y values as y[0..ny-1],
//  and the z values are in the 2-d array z[][].
//
//
// BASE_CONT    draw contour plot at bottom xy plane
// TOP_CONT     draw contour plot at top xy plane (not implemented)
// SURF_CONT    draw contour at surface
// FACETED      each square that makes up the surface is faceted
// DRAW_SIDES   draw sides
// MAG_COLOR    the surface is colored according to the value of z;
//               if not defined, the surface is colored according to the
//               intensity of the reflected light in the surface from a light
//               source whose position is set using c_pllightsource()
//
// or any bitwise combination, e.g. "MAG_COLOR | SURF_CONT | BASE_CONT"
//
// indexymin and indexymax are arrays which specify the y index range
// (following the convention that the upper range limit is one more than
// actual index limit) for an x index range of indexxmin, indexxmax.
// This code is a complete departure from the approach taken in the old version
// of this routine. Formerly to code attempted to use the logic for the hidden
// line algorithm to draw the hidden surface. This was really hard. This code
// below uses a simple back to front (painters) algorithm. All the
// triangles are drawn.
//
// There are multitude of ways this code could be optimized. Given the
// problems with the old code, I tried to focus on clarity here.
//--------------------------------------------------------------------------

void
c_plsurf3dl( PLFLT_VECTOR x, PLFLT_VECTOR y, PLFLT_MATRIX z, PLINT nx, PLINT ny,
             PLINT opt, PLFLT_VECTOR clevel, PLINT nlevel,
             PLINT indexxmin, PLINT indexxmax, PLINT_VECTOR indexymin, PLINT_VECTOR indexymax )
{
    plfsurf3dl( x, y, plf2ops_c(), (PLPointer) z, nx, ny,
        opt, clevel, nlevel, indexxmin, indexxmax, indexymin, indexymax );
}

void
plfsurf3dl( PLFLT_VECTOR x, PLFLT_VECTOR y, PLF2OPS zops, PLPointer zp, PLINT nx, PLINT ny,
            PLINT opt, PLFLT_VECTOR clevel, PLINT nlevel,
            PLINT indexxmin, PLINT indexxmax, PLINT_VECTOR indexymin, PLINT_VECTOR indexymax )
{
    PLFLT      cxx, cxy, cyx, cyy, cyz;
    PLINT      i, j, k;
    PLINT      ixDir, ixOrigin, iyDir, iyOrigin, nFast, nSlow;
    PLINT      ixFast, ixSlow, iyFast, iySlow;
    PLINT      iFast, iSlow;
    PLFLT      xmin, xmax, ymin, ymax, zmin, zmax, zscale;
    PLFLT      xm, ym, zm;
    PLINT      ixmin = 0, ixmax = nx, iymin = 0, iymax = ny;
    PLFLT      xx[3], yy[3], zz[3];
    PLFLT      px[4], py[4], pz[4];
    CONT_LEVEL *cont, *clev;
    CONT_LINE  *cline;
    int        ct, ix, iy, iftriangle;
    PLINT      color = plsc->icol0;
    PLFLT      width = plsc->width;
    PLFLT      ( *getz )( PLPointer, PLINT, PLINT ) = zops->get;

    if ( plsc->level < 3 )
    {
        myabort( "plsurf3dl: Please set up window first" );
        return;
    }

    if ( nx <= 0 || ny <= 0 )
    {
        myabort( "plsurf3dl: Bad array dimensions." );
        return;
    }

    //
    // Don't use the data z value to scale the color, use the z axis
    // values set by plw3d()
    //
    // plMinMax2dGrid(z, nx, ny, &fc_maxz, &fc_minz);
    //

    fc_minz = plsc->ranmi;
    fc_maxz = plsc->ranma;
    if ( fc_maxz == fc_minz )
    {
        plwarn( "plsurf3dl: Maximum and minimum Z values are equal! \"fixing\"..." );
        fc_maxz = fc_minz + 1e-6;
    }

    if ( opt & MAG_COLOR )
        falsecolor = 1;
    else
        falsecolor = 0;

    plP_gdom( &xmin, &xmax, &ymin, &ymax );
    plP_grange( &zscale, &zmin, &zmax );
    if ( zmin > zmax )
    {
        PLFLT t = zmin;
        zmin = zmax;
        zmax = t;
    }

    // Check that points in x and in y are strictly increasing  and in range

    for ( i = 0; i < nx - 1; i++ )
    {
        if ( x[i] >= x[i + 1] )
        {
            myabort( "plsurf3dl: X array must be strictly increasing" );
            return;
        }
        if ( x[i] < xmin && x[i + 1] >= xmin )
            ixmin = i;
        if ( x[i + 1] > xmax && x[i] <= xmax )
            ixmax = i + 2;
    }
    for ( i = 0; i < ny - 1; i++ )
    {
        if ( y[i] >= y[i + 1] )
        {
            myabort( "plsurf3dl: Y array must be strictly increasing" );
            return;
        }
        if ( y[i] < ymin && y[i + 1] >= ymin )
            iymin = i;
        if ( y[i + 1] > ymax && y[i] <= ymax )
            iymax = i + 2;
    }

    // get the viewing parameters
    plP_gw3wc( &cxx, &cxy, &cyx, &cyy, &cyz );

    // we're going to draw from back to front

    // iFast will index the dominant (fastest changing) dimension
    // iSlow will index the slower changing dimension
    //
    // iX indexes the X dimension
    // iY indexes the Y dimension

    // get direction for X
    if ( cxy >= 0 )
    {
        ixDir    = 1;     // direction in X
        ixOrigin = ixmin; // starting point
    }
    else
    {
        ixDir    = -1;
        ixOrigin = ixmax - 1;
    }
    // get direction for Y
    if ( cxx >= 0 )
    {
        iyDir    = -1;
        iyOrigin = iymax - 1;
    }
    else
    {
        iyDir    = 1;
        iyOrigin = iymin;
    }
    // figure out which dimension is dominant
    if ( fabs( cxx ) > fabs( cxy ) )
    {
        // X is dominant
        nFast = ixmax - ixmin;  // samples in the Fast direction
        nSlow = iymax - iymin;  // samples in the Slow direction

        ixFast = ixDir; ixSlow = 0;
        iyFast = 0;     iySlow = iyDir;
    }
    else
    {
        nFast = iymax - iymin;
        nSlow = ixmax - ixmin;

        ixFast = 0;     ixSlow = ixDir;
        iyFast = iyDir; iySlow = 0;
    }

    // we've got to draw the background grid first, hidden line code has to draw it last
    if ( zbflg )
    {
        PLFLT bx[3], by[3], bz[3];
        PLFLT tick = zbtck, tp;
        PLINT nsub = 0;

        // get the tick spacing
        pldtik( zmin, zmax, &tick, &nsub, FALSE );

        // determine the vertices for the background grid line
        bx[0] = ( ixOrigin != ixmin && ixFast == 0 ) || ixFast > 0 ? xmax : xmin;
        by[0] = ( iyOrigin != iymin && iyFast == 0 ) || iyFast > 0 ? ymax : ymin;
        bx[1] = ixOrigin != ixmin ? xmax : xmin;
        by[1] = iyOrigin != iymin ? ymax : ymin;
        bx[2] = ( ixOrigin != ixmin && ixSlow == 0 ) || ixSlow > 0 ? xmax : xmin;
        by[2] = ( iyOrigin != iymin && iySlow == 0 ) || iySlow > 0 ? ymax : ymin;

        plwidth( zbwidth );
        plcol0( zbcol );
        for ( tp = tick * floor( zmin / tick ) + tick; tp <= zmax; tp += tick )
        {
            bz[0] = bz[1] = bz[2] = tp;
            plline3( 3, bx, by, bz );
        }
        // draw the vertical line at the back corner
        bx[0] = bx[1];
        by[0] = by[1];
        bz[0] = zmin;
        plline3( 2, bx, by, bz );
        plwidth( width );
        plcol0( color );
    }

    // If enabled, draw the contour at the base

    // The contour plotted at the base will be identical to the one obtained
    // with c_plcont(). The contour plotted at the surface is simple minded, but
    // can be improved by using the contour data available.
    //

    if ( clevel != NULL && opt & BASE_CONT )
    {
#define NPTS    100
        int      np = NPTS;
        PLFLT    **zstore;
        PLcGrid2 cgrid2;
        PLFLT    *zzloc = (PLFLT *) malloc( (size_t) NPTS * sizeof ( PLFLT ) );
        if ( zzloc == NULL )
            plexit( "plsurf3dl: Insufficient memory" );

        // get the contour lines

        // prepare cont_store input
        cgrid2.nx = nx;
        cgrid2.ny = ny;
        plAlloc2dGrid( &cgrid2.xg, nx, ny );
        plAlloc2dGrid( &cgrid2.yg, nx, ny );
        plAlloc2dGrid( &zstore, nx, ny );

        for ( i = indexxmin; i < indexxmax; i++ )
        {
            for ( j = 0; j < indexymin[i]; j++ )
            {
                cgrid2.xg[i][j] = x[i];
                cgrid2.yg[i][j] = y[indexymin[i]];
                zstore[i][j]    = getz( zp, i, indexymin[i] );
            }
            for ( j = indexymin[i]; j < indexymax[i]; j++ )
            {
                cgrid2.xg[i][j] = x[i];
                cgrid2.yg[i][j] = y[j];
                zstore[i][j]    = getz( zp, i, j );
            }
            for ( j = indexymax[i]; j < ny; j++ )
            {
                cgrid2.xg[i][j] = x[i];
                cgrid2.yg[i][j] = y[indexymax[i] - 1];
                zstore[i][j]    = getz( zp, i, indexymax[i] - 1 );
            }
        }
        // Fill cont structure with contours.
        cont_store( (PLFLT_MATRIX) zstore, nx, ny, indexxmin + 1, indexxmax, 1, ny,
            clevel, nlevel, pltr2, (void *) &cgrid2, &cont );

        // Free the 2D input arrays to cont_store since not needed any more.
        plFree2dGrid( zstore, nx, ny );
        plFree2dGrid( cgrid2.xg, nx, ny );
        plFree2dGrid( cgrid2.yg, nx, ny );

        // follow the contour levels and lines
        clev = cont;
        do  // for each contour level
        {
            cline = clev->line;
            do  // there are several lines that make up the contour
            {
                if ( cline->npts > np )
                {
                    np = cline->npts;
                    if ( ( zzloc = (PLFLT *) realloc( zzloc, (size_t) np * sizeof ( PLFLT ) ) ) == NULL )
                    {
                        plexit( "plsurf3dl: Insufficient memory" );
                    }
                }
                for ( j = 0; j < cline->npts; j++ )
                    zzloc[j] = plsc->ranmi;
                if ( cline->npts > 0 )
                {
                    plcol1( ( clev->level - fc_minz ) / ( fc_maxz - fc_minz ) );
                    plline3( cline->npts, cline->x, cline->y, zzloc );
                }
                cline = cline->next;
            }
            while ( cline != NULL );
            clev = clev->next;
        }
        while ( clev != NULL );

        cont_clean_store( cont ); // now release the memory
        free( zzloc );
    }

    // Now we can iterate over the grid drawing the quads
    for ( iSlow = 0; iSlow < nSlow - 1; iSlow++ )
    {
        for ( iFast = 0; iFast < nFast - 1; iFast++ )
        {
            // get the 4 corners of the Quad, which are
            //
            //       0--2
            //       |  |
            //       1--3
            //

            xm = ym = zm = 0.;

            iftriangle = 1;
            for ( i = 0; i < 2; i++ )
            {
                for ( j = 0; j < 2; j++ )
                {
                    // we're transforming from Fast/Slow coordinates to x/y coordinates
                    // note, these are the indices, not the values
                    ix = ixFast * ( iFast + i ) + ixSlow * ( iSlow + j ) + ixOrigin;
                    iy = iyFast * ( iFast + i ) + iySlow * ( iSlow + j ) + iyOrigin;

                    if ( indexxmin <= ix && ix < indexxmax &&
                         indexymin[ix] <= iy && iy < indexymax[ix] )
                    {
                        xm += px[2 * i + j] = x[ix];
                        ym += py[2 * i + j] = y[iy];
                        zm += pz[2 * i + j] = getz( zp, ix, iy );
                    }
                    else
                    {
                        iftriangle = 0;
                        break;
                    }
                }
                if ( iftriangle == 0 )
                    break;
            }

            if ( iftriangle == 0 )
                continue;
            // the "mean point" of the quad, common to all four triangles
            // -- perhaps not a good thing to do for the light shading

            xm /= 4.; ym /= 4.; zm /= 4.;

            // now draw the quad as four triangles

            for ( i = 1; i < 3; i++ )
            {
                for ( j = 0; j < 4; j += 3 )
                {
                    shade_triangle( px[j], py[j], pz[j], xm, ym, zm, px[i], py[i], pz[i] );
                }
            }

            // After shading completed for a quad, render surface contours.
            if ( clevel != NULL && ( opt & SURF_CONT ) )
            {
                for ( i = 1; i < 3; i++ )
                {
                    for ( j = 0; j < 4; j += 3 )
#define min3( a, b, c )    ( MIN( ( MIN( a, b ) ), c ) )
#define max3( a, b, c )    ( MAX( ( MAX( a, b ) ), c ) )

                    {
                        for ( k = 0; k < nlevel; k++ )
                        {
                            if ( clevel[k] >= min3( pz[i], zm, pz[j] ) && clevel[k] < max3( pz[i], zm, pz[j] ) )
                            {
                                ct = 0;
                                if ( clevel[k] >= MIN( pz[i], zm ) && clevel[k] < MAX( pz[i], zm ) )     // p0-pm
                                {
                                    xx[ct] = ( ( clevel[k] - pz[i] ) * ( xm - px[i] ) ) / ( zm - pz[i] ) + px[i];
                                    yy[ct] = ( ( clevel[k] - pz[i] ) * ( ym - py[i] ) ) / ( zm - pz[i] ) + py[i];
                                    ct++;
                                }

                                if ( clevel[k] >= MIN( pz[i], pz[j] ) && clevel[k] < MAX( pz[i], pz[j] ) )     // p0-p1
                                {
                                    xx[ct] = ( ( clevel[k] - pz[i] ) * ( px[j] - px[i] ) ) / ( pz[j] - pz[i] ) + px[i];
                                    yy[ct] = ( ( clevel[k] - pz[i] ) * ( py[j] - py[i] ) ) / ( pz[j] - pz[i] ) + py[i];
                                    ct++;
                                }

                                if ( clevel[k] >= MIN( pz[j], zm ) && clevel[k] < MAX( pz[j], zm ) )     // p1-pm
                                {
                                    xx[ct] = ( ( clevel[k] - pz[j] ) * ( xm - px[j] ) ) / ( zm - pz[j] ) + px[j];
                                    yy[ct] = ( ( clevel[k] - pz[j] ) * ( ym - py[j] ) ) / ( zm - pz[j] ) + py[j];
                                    ct++;
                                }

                                if ( ct == 2 )
                                {
                                    // yes, xx and yy are the intersection points of the triangle with
                                    // the contour line -- draw a straight line betweeen the points
                                    // -- at the end this will make up the contour line

                                    // surface contour with color set by user
                                    plcol0( color );
                                    zz[0] = zz[1] = clevel[k];
                                    plline3( 2, xx, yy, zz );

                                    // don't break; one triangle can span various contour levels
                                }
                                else
                                    plwarn( "plsurf3dl: ***ERROR***\n" );
                            }
                        }
                    }
                }
            }
        }
    }

    if ( opt & FACETED )
    {
        plcol0( 0 );
        plfplot3dcl( x, y, zops, zp, nx, ny, MESH | DRAW_LINEXY, NULL, 0,
            indexxmin, indexxmax, indexymin, indexymax );
    }

    if ( opt & DRAW_SIDES ) // the sides look ugly !!!
    {                       // draw one more row with all the Z's set to zmin
        plP_grange( &zscale, &zmin, &zmax );

        iSlow      = nSlow - 1;
        iftriangle = 1;
        for ( iFast = 0; iFast < nFast - 1; iFast++ )
        {
            for ( i = 0; i < 2; i++ )
            {
                ix = ixFast * ( iFast + i ) + ixSlow * iSlow + ixOrigin;
                iy = iyFast * ( iFast + i ) + iySlow * iSlow + iyOrigin;
                if ( indexxmin <= ix && ix < indexxmax &&
                     indexymin[ix] <= iy && iy < indexymax[ix] )
                {
                    px[2 * i] = x[ix];
                    py[2 * i] = y[iy];
                    pz[2 * i] = getz( zp, ix, iy );
                }
                else
                {
                    iftriangle = 0;
                    break;
                }
            }
            if ( iftriangle == 0 )
                break;
            // now draw the quad as two triangles (4 might be better)

            shade_triangle( px[0], py[0], pz[0], px[2], py[2], pz[2], px[0], py[0], zmin );
            shade_triangle( px[2], py[2], pz[2], px[2], py[2], zmin, px[0], py[0], zmin );
        }

        iFast      = nFast - 1;
        iftriangle = 1;
        for ( iSlow = 0; iSlow < nSlow - 1; iSlow++ )
        {
            for ( i = 0; i < 2; i++ )
            {
                ix = ixFast * iFast + ixSlow * ( iSlow + i ) + ixOrigin;
                iy = iyFast * iFast + iySlow * ( iSlow + i ) + iyOrigin;
                if ( indexxmin <= ix && ix < indexxmax &&
                     indexymin[ix] <= iy && iy < indexymax[ix] )
                {
                    px[2 * i] = x[ix];
                    py[2 * i] = y[iy];
                    pz[2 * i] = getz( zp, ix, iy );
                }
                else
                {
                    iftriangle = 0;
                    break;
                }
            }
            if ( iftriangle == 0 )
                break;

            // now draw the quad as two triangles (4 might be better)
            shade_triangle( px[0], py[0], pz[0], px[2], py[2], pz[2], px[0], py[0], zmin );
            shade_triangle( px[2], py[2], pz[2], px[2], py[2], zmin, px[0], py[0], zmin );
        }
    }
}

//--------------------------------------------------------------------------
// void plot3d(x, y, z, nx, ny, opt, side)
//
// Plots a 3-d representation of the function z[x][y]. The x values
// are stored as x[0..nx-1], the y values as y[0..ny-1], and the z
// values are in the 2-d array z[][]. The integer "opt" specifies:
// see plot3dcl() below
//--------------------------------------------------------------------------

void
c_plot3d( PLFLT_VECTOR x, PLFLT_VECTOR y, PLFLT_MATRIX z,
          PLINT nx, PLINT ny, PLINT opt, PLBOOL side )
{
    plfplot3dc( x, y, plf2ops_c(), (PLPointer) z, nx, ny, opt | ( side != 0 ? DRAW_SIDES : 0 ), NULL, 0 );
}

void
plfplot3d( PLFLT_VECTOR x, PLFLT_VECTOR y, PLF2OPS zops, PLPointer zp,
           PLINT nx, PLINT ny, PLINT opt, PLBOOL side )
{
    plfplot3dc( x, y, zops, zp, nx, ny, opt | ( side != 0 ? DRAW_SIDES : 0 ), NULL, 0 );
}

//--------------------------------------------------------------------------
// void plot3dc(x, y, z, nx, ny, opt, clevel, nlevel)
//
// Plots a 3-d representation of the function z[x][y]. The x values
// are stored as x[0..nx-1], the y values as y[0..ny-1], and the z
// values are in the 2-d array z[][]. The integer "opt" specifies:
// see plot3dcl() below
//--------------------------------------------------------------------------

void
c_plot3dc( PLFLT_VECTOR x, PLFLT_VECTOR y, PLFLT_MATRIX z,
           PLINT nx, PLINT ny, PLINT opt,
           PLFLT_VECTOR clevel, PLINT nlevel )
{
    plfplot3dcl( x, y, plf2ops_c(), (PLPointer) z, nx, ny, opt, clevel, nlevel, 0, 0, NULL, NULL );
}

void
plfplot3dc( PLFLT_VECTOR x, PLFLT_VECTOR y, PLF2OPS zops, PLPointer zp,
            PLINT nx, PLINT ny, PLINT opt, PLFLT_VECTOR clevel, PLINT nlevel )
{
    plfplot3dcl( x, y, zops, zp, nx, ny, opt, clevel, nlevel, 0, 0, NULL, NULL );
}

//--------------------------------------------------------------------------
// void plot3dcl(x, y, z, nx, ny, opt, clevel, nlevel,
//       indexxmin, indexxmax, indexymin, indexymax)
//
// Plots a 3-d representation of the function z[x][y]. The x values
// are stored as x[0..nx-1], the y values as y[0..ny-1], and the z
// values are in the 2-d array z[][]. The integer "opt" specifies:
//
//  DRAW_LINEX :  Draw lines parallel to x-axis
//  DRAW_LINEY :  Draw lines parallel to y-axis
//  DRAW_LINEXY:  Draw lines parallel to both axes
//  MAG_COLOR:    Magnitude coloring of wire frame
//  BASE_CONT:    Draw contour at bottom xy plane
//  TOP_CONT:     Draw contour at top xy plane (not yet)
//  DRAW_SIDES:   Draw sides around the plot
//  MESH:       Draw the "under" side of the plot
//
// or any bitwise combination, e.g. "MAG_COLOR | DRAW_LINEX"
// indexymin and indexymax are arrays which specify the y index limits
// (following the convention that the upper range limit is one more than
// actual index limit) for an x index range of indexxmin, indexxmax.
//--------------------------------------------------------------------------

void
c_plot3dcl( PLFLT_VECTOR x, PLFLT_VECTOR y, PLFLT_MATRIX z,
            PLINT nx, PLINT ny, PLINT opt,
            PLFLT_VECTOR clevel, PLINT nlevel,
            PLINT indexxmin, PLINT indexxmax, PLINT_VECTOR indexymin, PLINT_VECTOR indexymax )
{
    plfplot3dcl( x, y, plf2ops_c(), (PLPointer) z, nx, ny,
        opt, clevel, nlevel, indexxmin, indexxmax, indexymin, indexymax );
}

//--------------------------------------------------------------------------
//! Plots a 3-d representation of the virtual function z, which is represented
//! by zops and zp.
//!
//! @param x The x values are stored as x[0..nx-1]
//! @param y The y values are stored as y[0..ny-1]
//! @param zops Pointer to a plf2ops_t structure of function pointers (see
//! plf2opts_t in plplot.h) which define how to perform various manipulations
//! (including retrieval) on the elements of the the 2D data field pointed to
//! by zp.  Pointers suitable for passing as zops can be obtained for some
//! predefined types of 2-d data storage by calling one of the plf2ops_*()
//! functions (see plf2ops.c) or you can create your own set for arbitrary 2-d
//! storage formats.
//! @param zp Pointer to 2D data field.  This pointer is passed to the
//! functions of zops whenever the 2D field needs to be manipulated.  The
//! combination of zops and zp provides total flexibility in how the underlying
//! data values are managed.
//! @param nx The number of values in x.
//! @param ny The number of values in y.
//! @param opt Specifies options for the plot.  It can be a bitwise OR-ing of
//! these:
//!   - DRAW_LINEX :  Draw lines parallel to x-axis
//!   - DRAW_LINEY :  Draw lines parallel to y-axis
//!   - DRAW_LINEXY:  Draw lines parallel to both axes
//!   - MAG_COLOR:    Magnitude coloring of wire frame
//!   - BASE_CONT:    Draw contour at bottom xy plane
//!   - TOP_CONT:     Draw contour at top xy plane (not yet)
//!   - DRAW_SIDES:   Draw sides around the plot
//!   - MESH:         Draw the "under" side of the plot
//! or any bitwise OR'd combination, e.g. "MAG_COLOR | DRAW_LINEX"
//! @param clevel z values at which to draw contours
//! @param nlevel Number of values in clevels
//! @param PL_UNUSED( indexxmin ) Not used.
//! @param PL_UNUSED( indexxmax ) Not used.
//! @param PL_UNUSED( indexymin ) Not used.
//! @param PL_UNUSED( indexymax ) Not used.
//!
//--------------------------------------------------------------------------

void
plfplot3dcl( PLFLT_VECTOR x, PLFLT_VECTOR y, PLF2OPS zops, PLPointer zp,
             PLINT nx, PLINT ny, PLINT opt,
             PLFLT_VECTOR clevel, PLINT nlevel,
             PLINT PL_UNUSED( indexxmin ), PLINT PL_UNUSED( indexxmax ), PLINT_VECTOR PL_UNUSED( indexymin ), PLINT_VECTOR PL_UNUSED( indexymax ) )
{
    PLFLT cxx, cxy, cyx, cyy, cyz;
    PLINT init, ix, iy, color;
    PLFLT width;
    PLFLT xmin, xmax, ymin, ymax, zmin, zmax, zscale;
    PLINT ixmin   = 0, ixmax = nx - 1, iymin = 0, iymax = ny - 1;
    PLINT clipped = 0, base_cont = 0, side = 0;
    PLFLT ( *getz )( PLPointer, PLINT, PLINT ) = zops->get;
    PLFLT *_x = NULL, *_y = NULL, **_z = NULL;
    PLFLT_VECTOR x_modified, y_modified;
    int i;

    pl3mode = 0;

    if ( plsc->level < 3 )
    {
        myabort( "plot3dcl: Please set up window first" );
        return;
    }

    if ( ( opt & 3 ) == 0 )
    {
        myabort( "plot3dcl: Bad option" );
        return;
    }

    if ( nx <= 0 || ny <= 0 )
    {
        myabort( "plot3dcl: Bad array dimensions." );
        return;
    }

    plP_gdom( &xmin, &xmax, &ymin, &ymax );
    plP_grange( &zscale, &zmin, &zmax );
    if ( zmin > zmax )
    {
        PLFLT t = zmin;
        zmin = zmax;
        zmax = t;
    }

// Check that points in x and in y are strictly increasing

    for ( i = 0; i < nx - 1; i++ )
    {
        if ( x[i] >= x[i + 1] )
        {
            myabort( "plot3dcl: X array must be strictly increasing" );
            return;
        }
    }
    for ( i = 0; i < ny - 1; i++ )
    {
        if ( y[i] >= y[i + 1] )
        {
            myabort( "plot3dcl: Y array must be strictly increasing" );
            return;
        }
    }

    if ( opt & MESH )
        pl3mode = 1;

    if ( opt & DRAW_SIDES )
        side = 1;

    // figure out the part of the data to use
    if ( xmin < x[0] )
        xmin = x[0];
    if ( xmax > x[nx - 1] )
        xmax = x[nx - 1];
    if ( ymin < y[0] )
        ymin = y[0];
    if ( ymax > y[ny - 1] )
        ymax = y[ny - 1];
    for ( ixmin = 0; ixmin < nx - 1 && x[ixmin + 1] < xmin; ixmin++ )
    {
    }
    for ( ixmax = nx - 1; ixmax > 0 && x[ixmax - 1] > xmax; ixmax-- )
    {
    }
    for ( iymin = 0; iymin < ny - 1 && y[iymin + 1] < ymin; iymin++ )
    {
    }
    for ( iymax = ny - 1; iymax > 0 && y[iymax - 1] > ymax; iymax-- )
    {
    }
    //fprintf(stderr, "(%d,%d) %d %d %d %d\n", nx, ny, ixmin, ixmax, iymin, iymax);
    // do we need to clip?
    if ( ixmin > 0 || ixmax < nx - 1 || iymin > 0 || iymax < ny - 1 )
    {
        // adjust the input so it stays within bounds
        int _nx = ixmax - ixmin + 1;
        int _ny = iymax - iymin + 1;
        PLFLT ty0, ty1, tx0, tx1;
        int j;

        if ( _nx <= 1 || _ny <= 1 )
        {
            myabort( "plot3dcl: selected x or y range has no data" );
            return;
        }

        // allocate storage for new versions of the input vectors
        if ( ( ( _x = (PLFLT *) malloc( (size_t) _nx * sizeof ( PLFLT ) ) ) == NULL ) ||
             ( ( _y = (PLFLT *) malloc( (size_t) _ny * sizeof ( PLFLT ) ) ) == NULL ) ||
             ( ( _z = (PLFLT **) malloc( (size_t) _nx * sizeof ( PLFLT* ) ) ) == NULL ) )
        {
            plexit( "c_plot3dcl: Insufficient memory" );
        }

        clipped = 1;

        // copy over the independent variables
        _x[0]       = xmin;
        _x[_nx - 1] = xmax;
        for ( i = 1; i < _nx - 1; i++ )
            _x[i] = x[ixmin + i];
        _y[0]       = ymin;
        _y[_ny - 1] = ymax;
        for ( i = 1; i < _ny - 1; i++ )
            _y[i] = y[iymin + i];

        // copy the data array so we can interpolate around the edges
        for ( i = 0; i < _nx; i++ )
        {
            if ( ( _z[i] = (PLFLT *) malloc( (size_t) _ny * sizeof ( PLFLT ) ) ) == NULL )
            {
                plexit( "c_plot3dcl: Insufficient memory" );
            }
        }

        // interpolation factors for the 4 edges
        ty0 = ( _y[0] - y[iymin] ) / ( y[iymin + 1] - y[iymin] );
        ty1 = ( _y[_ny - 1] - y[iymax - 1] ) / ( y[iymax] - y[iymax - 1] );
        tx0 = ( _x[0] - x[ixmin] ) / ( x[ixmin + 1] - x[ixmin] );
        tx1 = ( _x[_nx - 1] - x[ixmax - 1] ) / ( x[ixmax] - x[ixmax - 1] );
        for ( i = 0; i < _nx; i++ )
        {
            if ( i == 0 )
            {
                _z[i][0] = getz( zp, ixmin, iymin ) * ( 1 - ty0 ) * ( 1 - tx0 ) + getz( zp, ixmin, iymin + 1 ) * ( 1 - tx0 ) * ty0
                           + getz( zp, ixmin + 1, iymin ) * tx0 * ( 1 - ty0 ) + getz( zp, ixmin + 1, iymin + 1 ) * tx0 * ty0;
                for ( j = 1; j < _ny - 1; j++ )
                    _z[i][j] = getz( zp, ixmin, iymin + j ) * ( 1 - tx0 ) + getz( zp, ixmin + 1, iymin + j ) * tx0;
                _z[i][_ny - 1] = getz( zp, ixmin, iymax - 1 ) * ( 1 - tx0 ) * ( 1 - ty1 ) + getz( zp, ixmin, iymax ) * ( 1 - tx0 ) * ty1
                                 + getz( zp, ixmin + 1, iymax - 1 ) * tx0 * ( 1 - ty1 ) + getz( zp, ixmin + 1, iymax ) * tx0 * ty1;
            }
            else if ( i == _nx - 1 )
            {
                _z[i][0] = getz( zp, ixmax - 1, iymin ) * ( 1 - tx1 ) * ( 1 - ty0 ) + getz( zp, ixmax - 1, iymin + 1 ) * ( 1 - tx1 ) * ty0
                           + getz( zp, ixmax, iymin ) * tx1 * ( 1 - ty0 ) + getz( zp, ixmax, iymin + 1 ) * tx1 * ty0;
                for ( j = 1; j < _ny - 1; j++ )
                    _z[i][j] = getz( zp, ixmax - 1, iymin + j ) * ( 1 - tx1 ) + getz( zp, ixmax, iymin + j ) * tx1;
                _z[i][_ny - 1] = getz( zp, ixmax - 1, iymax - 1 ) * ( 1 - tx1 ) * ( 1 - ty1 ) + getz( zp, ixmax, iymax ) * ( 1 - tx1 ) * ty1
                                 + getz( zp, ixmax, iymax - 1 ) * tx1 * ( 1 - ty1 ) + getz( zp, ixmax, iymax ) * tx1 * ty1;
            }
            else
            {
                _z[i][0] = getz( zp, ixmin + i, iymin ) * ( 1 - ty0 ) + getz( zp, ixmin + i, iymin + 1 ) * ty0;
                for ( j = 1; j < _ny - 1; j++ )
                    _z[i][j] = getz( zp, ixmin + i, iymin + j );
                _z[i][_ny - 1] = getz( zp, ixmin + i, iymax - 1 ) * ( 1 - ty1 ) + getz( zp, ixmin + i, iymax ) * ty1;
            }
            for ( j = 0; j < _ny; j++ )
            {
                if ( _z[i][j] < zmin )
                    _z[i][j] = zmin;
                else if ( _z[i][j] > zmax )
                    _z[i][j] = zmax;
            }
        }
        // replace the input with our clipped versions
        zp   = (PLPointer) _z;
        getz = plf2ops_c()->get;
        nx   = _nx;
        ny   = _ny;
        // Do not want to modify input x and y (const modifier)
        x_modified = (PLFLT_VECTOR) _x;
        y_modified = (PLFLT_VECTOR) _y;
    }
    else
    {
        x_modified = x;
        y_modified = y;
    }

    // From here on must use x_modified and y_modified rather than
    // x and y.
    if ( ( opt & BASE_CONT ) || ( opt & TOP_CONT ) || ( opt & MAG_COLOR ) )
    {
        //
        // Don't use the data z value to scale the color, use the z axis
        // values set by plw3d()
        //
        // plMinMax2dGrid(z, nx, ny, &fc_maxz, &fc_minz);
        //

        fc_minz = plsc->ranmi;
        fc_maxz = plsc->ranma;

        if ( fc_maxz == fc_minz )
        {
            plwarn( "plot3dcl: Maximum and minimum Z values are equal! \"fixing\"..." );
            fc_maxz = fc_minz + 1e-6;
        }
    }

    if ( opt & BASE_CONT )    // If enabled, draw the contour at the base.
    {
        if ( clevel != NULL && nlevel != 0 )
        {
            base_cont = 1;
            // even if MESH is not set, "set it",
            // as the base contour can only be done in this case
            pl3mode = 1;
        }
    }

    if ( opt & MAG_COLOR )    // If enabled, use magnitude colored wireframe
    {
        if ( ( ctmp = (PLFLT *) malloc( (size_t) ( 2 * MAX( nx, ny ) ) * sizeof ( PLFLT ) ) ) == NULL )
        {
            plexit( "c_plot3dcl: Insufficient memory" );
        }
    }
    else
        ctmp = NULL;

    // next logic only knows opt = 1 | 2 | 3, make sure that it only gets that
    opt &= DRAW_LINEXY;

    // Allocate work arrays

    utmp = (PLINT *) malloc( (size_t) ( 2 * MAX( nx, ny ) ) * sizeof ( PLINT ) );
    vtmp = (PLINT *) malloc( (size_t) ( 2 * MAX( nx, ny ) ) * sizeof ( PLINT ) );

    if ( !utmp || !vtmp )
        myexit( "plot3dcl: Out of memory." );

    plP_gw3wc( &cxx, &cxy, &cyx, &cyy, &cyz );
    init = 1;
// Call 3d line plotter.  Each viewing quadrant
// (perpendicular to x-y plane) must be handled separately.
    if ( cxx >= 0.0 && cxy <= 0.0 )
    {
        if ( opt == DRAW_LINEY )
            plt3zz( 1, ny, 1, -1, -opt, &init, x_modified, y_modified, zops, zp, nx, ny, utmp, vtmp, ctmp );
        else
        {
            for ( iy = 2; iy <= ny; iy++ )
                plt3zz( 1, iy, 1, -1, -opt, &init, x_modified, y_modified, zops, zp, nx, ny, utmp, vtmp, ctmp );
        }
        if ( opt == DRAW_LINEX )
            plt3zz( 1, ny, 1, -1, opt, &init, x_modified, y_modified, zops, zp, nx, ny, utmp, vtmp, ctmp );
        else
        {
            for ( ix = 1; ix <= nx - 1; ix++ )
                plt3zz( ix, ny, 1, -1, opt, &init, x_modified, y_modified, zops, zp, nx, ny, utmp, vtmp, ctmp );
        }
    }

    else if ( cxx <= 0.0 && cxy <= 0.0 )
    {
        if ( opt == DRAW_LINEX )
            plt3zz( nx, ny, -1, -1, opt, &init, x_modified, y_modified, zops, zp, nx, ny, utmp, vtmp, ctmp );
        else
        {
            for ( ix = 2; ix <= nx; ix++ )
                plt3zz( ix, ny, -1, -1, opt, &init, x_modified, y_modified, zops, zp, nx, ny, utmp, vtmp, ctmp );
        }
        if ( opt == DRAW_LINEY )
            plt3zz( nx, ny, -1, -1, -opt, &init, x_modified, y_modified, zops, zp, nx, ny, utmp, vtmp, ctmp );
        else
        {
            for ( iy = ny; iy >= 2; iy-- )
                plt3zz( nx, iy, -1, -1, -opt, &init, x_modified, y_modified, zops, zp, nx, ny, utmp, vtmp, ctmp );
        }
    }

    else if ( cxx <= 0.0 && cxy >= 0.0 )
    {
        if ( opt == DRAW_LINEY )
            plt3zz( nx, 1, -1, 1, -opt, &init, x_modified, y_modified, zops, zp, nx, ny, utmp, vtmp, ctmp );
        else
        {
            for ( iy = ny - 1; iy >= 1; iy-- )
                plt3zz( nx, iy, -1, 1, -opt, &init, x_modified, y_modified, zops, zp, nx, ny, utmp, vtmp, ctmp );
        }
        if ( opt == DRAW_LINEX )
            plt3zz( nx, 1, -1, 1, opt, &init, x_modified, y_modified, zops, zp, nx, ny, utmp, vtmp, ctmp );
        else
        {
            for ( ix = nx; ix >= 2; ix-- )
                plt3zz( ix, 1, -1, 1, opt, &init, x_modified, y_modified, zops, zp, nx, ny, utmp, vtmp, ctmp );
        }
    }

    else if ( cxx >= 0.0 && cxy >= 0.0 )
    {
        if ( opt == DRAW_LINEX )
            plt3zz( 1, 1, 1, 1, opt, &init, x_modified, y_modified, zops, zp, nx, ny, utmp, vtmp, ctmp );
        else
        {
            for ( ix = nx - 1; ix >= 1; ix-- )
                plt3zz( ix, 1, 1, 1, opt, &init, x_modified, y_modified, zops, zp, nx, ny, utmp, vtmp, ctmp );
        }
        if ( opt == DRAW_LINEY )
            plt3zz( 1, 1, 1, 1, -opt, &init, x_modified, y_modified, zops, zp, nx, ny, utmp, vtmp, ctmp );
        else
        {
            for ( iy = 1; iy <= ny - 1; iy++ )
                plt3zz( 1, iy, 1, 1, -opt, &init, x_modified, y_modified, zops, zp, nx, ny, utmp, vtmp, ctmp );
        }
    }

    // draw contour at the base. Not 100%! Why?

    if ( base_cont )
    {
        int np = NPTS, j;
        CONT_LEVEL *cont, *clev;
        CONT_LINE *cline;

        PLINT *uu = (PLINT *) malloc( (size_t) NPTS * sizeof ( PLINT ) );
        PLINT *vv = (PLINT *) malloc( (size_t) NPTS * sizeof ( PLINT ) );
        // prepare cont_store input
        PLFLT **zstore;
        PLcGrid2 cgrid2;

        if ( ( uu == NULL ) || ( vv == NULL ) )
        {
            plexit( "c_plot3dcl: Insufficient memory" );
        }

        cgrid2.nx = nx;
        cgrid2.ny = ny;
        plAlloc2dGrid( &cgrid2.xg, nx, ny );
        plAlloc2dGrid( &cgrid2.yg, nx, ny );
        plAlloc2dGrid( &zstore, nx, ny );

        for ( i = 0; i < nx; i++ )
        {
            for ( j = 0; j < ny; j++ )
            {
                cgrid2.xg[i][j] = x_modified[i];
                cgrid2.yg[i][j] = y_modified[j];
                zstore[i][j]    = getz( zp, i, j );
            }
        }

        pl3upv = 0;

        // Fill cont structure with contours.
        cont_store( (PLFLT_MATRIX) zstore, nx, ny, 1, nx, 1, ny,
            clevel, nlevel, pltr2, (void *) &cgrid2, &cont );

        // Free the 2D input arrays to cont_store since not needed any more.
        plFree2dGrid( zstore, nx, ny );
        plFree2dGrid( cgrid2.xg, nx, ny );
        plFree2dGrid( cgrid2.yg, nx, ny );

        // follow the contour levels and lines
        clev = cont;
        do  // for each contour level
        {
            cline = clev->line;
            do  // there are several lines that make up each contour
            {
                int cx, k, l, m, start, end;
                PLFLT tx, ty;
                if ( cline->npts > np )
                {
                    np = cline->npts;
                    if ( ( ( uu = (PLINT *) realloc( uu, (size_t) np * sizeof ( PLINT ) ) ) == NULL ) ||
                         ( ( vv = (PLINT *) realloc( vv, (size_t) np * sizeof ( PLINT ) ) ) == NULL ) )
                    {
                        plexit( "c_plot3dcl: Insufficient memory" );
                    }
                }

                // the hidden line plotter plnxtv() only works OK if the x points are in increasing order.
                // As this does not always happens, the situation must be detected and the line segment
                // must be reversed before being plotted
                i = 0;
                if ( cline->npts > 0 )
                {
                    do
                    {
                        plcol1( ( clev->level - fc_minz ) / ( fc_maxz - fc_minz ) );
                        cx = plP_wcpcx( plP_w3wcx( cline->x[i], cline->y[i], plsc->ranmi ) );
                        for ( j = i; j < cline->npts; j++ ) // convert to 2D coordinates
                        {
                            uu[j] = plP_wcpcx( plP_w3wcx( cline->x[j], cline->y[j], plsc->ranmi ) );
                            vv[j] = plP_wcpcy( plP_w3wcy( cline->x[j], cline->y[j], plsc->ranmi ) );
                            if ( uu[j] < cx ) // find turn back point
                                break;
                            else
                                cx = uu[j];
                        }
                        plnxtv( &uu[i], &vv[i], NULL, j - i, 0 ); // plot line with increasing x

                        if ( j < cline->npts )                    // line not yet finished,
                        {
                            start = j - 1;
                            for ( i = start; i < cline->npts; i++ ) // search turn forward point
                            {
                                uu[i] = plP_wcpcx( plP_w3wcx( cline->x[i], cline->y[i], plsc->ranmi ) );
                                if ( uu[i] > cx )
                                    break;
                                else
                                    cx = uu[i];
                            }
                            end = i - 1;

                            for ( k = 0; k < ( end - start + 1 ) / 2; k++ ) // reverse line segment
                            {
                                l           = start + k;
                                m           = end - k;
                                tx          = cline->x[l];
                                ty          = cline->y[l];
                                cline->x[l] = cline->x[m];
                                cline->y[l] = cline->y[m];
                                cline->x[m] = tx;
                                cline->y[m] = ty;
                            }

                            // convert to 2D coordinates
                            for ( j = start; j <= end; j++ )
                            {
                                uu[j] = plP_wcpcx( plP_w3wcx( cline->x[j], cline->y[j], plsc->ranmi ) );
                                vv[j] = plP_wcpcy( plP_w3wcy( cline->x[j], cline->y[j], plsc->ranmi ) );
                            }
                            plnxtv( &uu[start], &vv[start], NULL, end - start + 1, 0 ); // and plot it

                            cline->x[end] = cline->x[start];
                            cline->y[end] = cline->y[start];
                            i             = end; // restart where it was left
                        }
                    } while ( j < cline->npts && i < cline->npts );
                }
                cline = cline->next;
            }
            while ( cline != NULL );
            clev = clev->next;
        }
        while ( clev != NULL );

        cont_clean_store( cont ); // now release the contour memory
        pl3upv = 1;
        free( uu );
        free( vv );
    }

// Finish up by drawing sides, background grid (both are optional)

    if ( side )
        plside3( x_modified, y_modified, zops, zp, nx, ny, opt );

    if ( zbflg )
    {
        color = plsc->icol0;
        width = plsc->width;
        plwidth( zbwidth );
        plcol0( zbcol );
        plgrid3( zbtck );
        plwidth( width );
        plcol0( color );
    }

    freework();

    if ( clipped )
    {
        free( _x );
        free( _y );
        for ( i = 0; i < nx; i++ )
            free( _z[i] );
        free( _z );
    }
}

//--------------------------------------------------------------------------
// void plxyindexlimits()
//
// Transform from y array limits to corresponding x array limits (or vice
// versa).
//
// N.B. we follow the convention here that all upper range limits are one
// more than the actual last index.
// instart (>= 0) through inn is the index range where
// the input inarray_min and inarray_max arrays are defined.
//
// outstart (>= 0), through outn (with outn <= outnmax) is the index range
// where the output outarray_min and outarray_max arrays are defined.
//
// In order to assure the transformation from y array limits to x array limits
// (or vice versa) is single-valued, this programme plaborts if the
// inarray_min array has a maximum or inarray_max array has a minimum.
//--------------------------------------------------------------------------

//static void
//plxyindexlimits( PLINT instart, PLINT inn,
//                 PLINT *inarray_min, PLINT *inarray_max,
//                 PLINT *outstart, PLINT *outn, PLINT outnmax,
//                 PLINT *outarray_min, PLINT *outarray_max )
//{
//    PLINT i, j;
//    if ( inn < 0 )
//    {
//        myabort( "plxyindexlimits: Must have instart >= 0" );
//        return;
//    }
//    if ( inn - instart <= 0 )
//    {
//        myabort( "plxyindexlimits: Must have at least 1 defined point" );
//        return;
//    }
//    *outstart = inarray_min[instart];
//    *outn     = inarray_max[instart];
//    for ( i = instart; i < inn; i++ )
//    {
//        *outstart = MIN( *outstart, inarray_min[i] );
//        *outn     = MAX( *outn, inarray_max[i] );
//        if ( i + 2 < inn )
//        {
//            if ( inarray_min[i] < inarray_min[i + 1] &&
//                 inarray_min[i + 1] > inarray_min[i + 2] )
//            {
//                myabort( "plxyindexlimits: inarray_min must not have a maximum" );
//                return;
//            }
//            if ( inarray_max[i] > inarray_max[i + 1] &&
//                 inarray_max[i + 1] < inarray_max[i + 2] )
//            {
//                myabort( "plxyindexlimits: inarray_max must not have a minimum" );
//                return;
//            }
//        }
//    }
//    if ( *outstart < 0 )
//    {
//        myabort( "plxyindexlimits: Must have all elements of inarray_min >= 0" );
//        return;
//    }
//    if ( *outn > outnmax )
//    {
//        myabort( "plxyindexlimits: Must have all elements of inarray_max <= outnmax" );
//        return;
//    }
//    for ( j = *outstart; j < *outn; j++ )
//    {
//        i = instart;
//        // Find first valid i for this j.
//        while ( i < inn && !( inarray_min[i] <= j && j < inarray_max[i] ) )
//            i++;
//        if ( i < inn )
//            outarray_min[j] = i;
//        else
//        {
//            myabort( "plxyindexlimits: bad logic; invalid i should never happen" );
//            return;
//        }
//        // Find next invalid i for this j.
//        while ( i < inn && ( inarray_min[i] <= j && j < inarray_max[i] ) )
//            i++;
//        outarray_max[j] = i;
//    }
//}

//--------------------------------------------------------------------------
// void plP_gzback()
//
// Get background parameters for 3d plot.
//--------------------------------------------------------------------------

void
plP_gzback( PLINT **zbf, PLINT **zbc, PLFLT **zbt, PLFLT **zbw )
{
    *zbf = &zbflg;
    *zbc = &zbcol;
    *zbt = &zbtck;
    *zbw = &zbwidth;
}

//--------------------------------------------------------------------------
// PLFLT plGetAngleToLight()
//
// Gets cos of angle between normal to a polygon and a light source.
// Requires at least 3 elements, forming non-parallel lines
// in the arrays.
//--------------------------------------------------------------------------

static PLFLT
plGetAngleToLight( PLFLT* x, PLFLT* y, PLFLT* z )
{
    PLFLT vx1, vx2, vy1, vy2, vz1, vz2;
    PLFLT px, py, pz;
    PLFLT vlx, vly, vlz;
    PLFLT mag1, mag2;
    PLFLT cosangle;

    vx1 = x[1] - x[0];
    vx2 = x[2] - x[1];
    vy1 = y[1] - y[0];
    vy2 = y[2] - y[1];
    vz1 = z[1] - z[0];
    vz2 = z[2] - z[1];

// Find vector perpendicular to the face
    px   = vy1 * vz2 - vz1 * vy2;
    py   = vz1 * vx2 - vx1 * vz2;
    pz   = vx1 * vy2 - vy1 * vx2;
    mag1 = px * px + py * py + pz * pz;

// Vectors were parallel!
    if ( mag1 == 0 )
        return 1;

    vlx  = xlight - x[0];
    vly  = ylight - y[0];
    vlz  = zlight - z[0];
    mag2 = vlx * vlx + vly * vly + vlz * vlz;
    if ( mag2 == 0 )
        return 1;

// Now have 3 vectors going through the first point on the given surface
    cosangle = fabs( ( vlx * px + vly * py + vlz * pz ) / sqrt( mag1 * mag2 ) );

// In case of numerical rounding
    if ( cosangle > 1 )
        cosangle = 1;
    return cosangle;
}

//--------------------------------------------------------------------------
// void plt3zz()
//
// Draws the next zig-zag line for a 3-d plot.  The data is stored in array
// z[][] as a function of x[] and y[], and is plotted out starting at index
// (x0,y0).
//
// Depending on the state of "flag", the sequence of data points sent to
// plnxtv is altered so as to allow cross-hatch plotting, or plotting
// parallel to either the x-axis or the y-axis.
//--------------------------------------------------------------------------

static void
plt3zz( PLINT x0, PLINT y0, PLINT dx, PLINT dy, PLINT flag, PLINT *init,
        PLFLT_VECTOR x, PLFLT_VECTOR y, PLF2OPS zops, PLPointer zp, PLINT nx, PLINT ny,
        PLINT *u, PLINT *v, PLFLT* c )
{
    PLINT n = 0;
    PLFLT x2d, y2d;
    PLFLT ( *getz )( PLPointer, PLINT, PLINT ) = zops->get;

    while ( 1 <= x0 && x0 <= nx && 1 <= y0 && y0 <= ny )
    {
        x2d  = plP_w3wcx( x[x0 - 1], y[y0 - 1], getz( zp, x0 - 1, y0 - 1 ) );
        y2d  = plP_w3wcy( x[x0 - 1], y[y0 - 1], getz( zp, x0 - 1, y0 - 1 ) );
        u[n] = plP_wcpcx( x2d );
        v[n] = plP_wcpcy( y2d );
        if ( c != NULL )
            c[n] = ( getz( zp, x0 - 1, y0 - 1 ) - fc_minz ) / ( fc_maxz - fc_minz );

        switch ( flag )
        {
        case -3:
            y0  += dy;
            flag = -flag;
            break;
        case -2:
            y0 += dy;
            break;
        case -1:
            y0  += dy;
            flag = -flag;
            break;
        case 1:
            x0 += dx;
            break;
        case 2:
            x0  += dx;
            flag = -flag;
            break;
        case 3:
            x0  += dx;
            flag = -flag;
            break;
        }
        n++;
    }

    if ( flag == 1 || flag == -2 )
    {
        if ( flag == 1 )
        {
            x0 -= dx;
            y0 += dy;
        }
        else if ( flag == -2 )
        {
            y0 -= dy;
            x0 += dx;
        }
        if ( 1 <= x0 && x0 <= nx && 1 <= y0 && y0 <= ny )
        {
            x2d  = plP_w3wcx( x[x0 - 1], y[y0 - 1], getz( zp, x0 - 1, y0 - 1 ) );
            y2d  = plP_w3wcy( x[x0 - 1], y[y0 - 1], getz( zp, x0 - 1, y0 - 1 ) );
            u[n] = plP_wcpcx( x2d );
            v[n] = plP_wcpcy( y2d );
            if ( c != NULL )
                c[n] = ( getz( zp, x0 - 1, y0 - 1 ) - fc_minz ) / ( fc_maxz - fc_minz );
            n++;
        }
    }

// All the setup is done.  Time to do the work.

    plnxtv( u, v, c, n, *init );
    *init = 0;
}

//--------------------------------------------------------------------------
// void plside3()
//
// This routine draws sides around the front of the 3d plot so that
// it does not appear to float.
//--------------------------------------------------------------------------

static void
plside3( PLFLT_VECTOR x, PLFLT_VECTOR y, PLF2OPS zops, PLPointer zp, PLINT nx, PLINT ny, PLINT opt )
{
    PLINT i;
    PLFLT cxx, cxy, cyx, cyy, cyz;
    PLFLT xmin, ymin, zmin, xmax, ymax, zmax, zscale;
    PLFLT tx, ty, ux, uy;
    PLFLT ( *getz )( PLPointer, PLINT, PLINT ) = zops->get;

    plP_gw3wc( &cxx, &cxy, &cyx, &cyy, &cyz );
    plP_gdom( &xmin, &xmax, &ymin, &ymax );
    plP_grange( &zscale, &zmin, &zmax );

// Get x, y coordinates of legs and plot

    if ( cxx >= 0.0 && cxy <= 0.0 )
    {
        if ( opt != 1 )
        {
            for ( i = 0; i < nx; i++ )
            {
                tx = plP_w3wcx( x[i], y[0], zmin );
                ty = plP_w3wcy( x[i], y[0], zmin );
                ux = plP_w3wcx( x[i], y[0], getz( zp, i, 0 ) );
                uy = plP_w3wcy( x[i], y[0], getz( zp, i, 0 ) );
                pljoin( tx, ty, ux, uy );
            }
        }
        if ( opt != 2 )
        {
            for ( i = 0; i < ny; i++ )
            {
                tx = plP_w3wcx( x[0], y[i], zmin );
                ty = plP_w3wcy( x[0], y[i], zmin );
                ux = plP_w3wcx( x[0], y[i], getz( zp, 0, i ) );
                uy = plP_w3wcy( x[0], y[i], getz( zp, 0, i ) );
                pljoin( tx, ty, ux, uy );
            }
        }
    }
    else if ( cxx <= 0.0 && cxy <= 0.0 )
    {
        if ( opt != 1 )
        {
            for ( i = 0; i < nx; i++ )
            {
                tx = plP_w3wcx( x[i], y[ny - 1], zmin );
                ty = plP_w3wcy( x[i], y[ny - 1], zmin );
                ux = plP_w3wcx( x[i], y[ny - 1], getz( zp, i, ny - 1 ) );
                uy = plP_w3wcy( x[i], y[ny - 1], getz( zp, i, ny - 1 ) );
                pljoin( tx, ty, ux, uy );
            }
        }
        if ( opt != 2 )
        {
            for ( i = 0; i < ny; i++ )
            {
                tx = plP_w3wcx( x[0], y[i], zmin );
                ty = plP_w3wcy( x[0], y[i], zmin );
                ux = plP_w3wcx( x[0], y[i], getz( zp, 0, i ) );
                uy = plP_w3wcy( x[0], y[i], getz( zp, 0, i ) );
                pljoin( tx, ty, ux, uy );
            }
        }
    }
    else if ( cxx <= 0.0 && cxy >= 0.0 )
    {
        if ( opt != 1 )
        {
            for ( i = 0; i < nx; i++ )
            {
                tx = plP_w3wcx( x[i], y[ny - 1], zmin );
                ty = plP_w3wcy( x[i], y[ny - 1], zmin );
                ux = plP_w3wcx( x[i], y[ny - 1], getz( zp, i, ny - 1 ) );
                uy = plP_w3wcy( x[i], y[ny - 1], getz( zp, i, ny - 1 ) );
                pljoin( tx, ty, ux, uy );
            }
        }
        if ( opt != 2 )
        {
            for ( i = 0; i < ny; i++ )
            {
                tx = plP_w3wcx( x[nx - 1], y[i], zmin );
                ty = plP_w3wcy( x[nx - 1], y[i], zmin );
                ux = plP_w3wcx( x[nx - 1], y[i], getz( zp, nx - 1, i ) );
                uy = plP_w3wcy( x[nx - 1], y[i], getz( zp, nx - 1, i ) );
                pljoin( tx, ty, ux, uy );
            }
        }
    }
    else if ( cxx >= 0.0 && cxy >= 0.0 )
    {
        if ( opt != 1 )
        {
            for ( i = 0; i < nx; i++ )
            {
                tx = plP_w3wcx( x[i], y[0], zmin );
                ty = plP_w3wcy( x[i], y[0], zmin );
                ux = plP_w3wcx( x[i], y[0], getz( zp, i, 0 ) );
                uy = plP_w3wcy( x[i], y[0], getz( zp, i, 0 ) );
                pljoin( tx, ty, ux, uy );
            }
        }
        if ( opt != 2 )
        {
            for ( i = 0; i < ny; i++ )
            {
                tx = plP_w3wcx( x[nx - 1], y[i], zmin );
                ty = plP_w3wcy( x[nx - 1], y[i], zmin );
                ux = plP_w3wcx( x[nx - 1], y[i], getz( zp, nx - 1, i ) );
                uy = plP_w3wcy( x[nx - 1], y[i], getz( zp, nx - 1, i ) );
                pljoin( tx, ty, ux, uy );
            }
        }
    }
}

//--------------------------------------------------------------------------
// void plgrid3()
//
// This routine draws a grid around the back side of the 3d plot with
// hidden line removal.
//--------------------------------------------------------------------------

static void
plgrid3( PLFLT tick )
{
    PLFLT xmin, ymin, zmin, xmax, ymax, zmax, zscale;
    PLFLT cxx, cxy, cyx, cyy, cyz, zmin_in, zmax_in;
    PLINT u[3], v[3];
    PLINT nsub = 0;
    PLFLT tp;

    plP_gw3wc( &cxx, &cxy, &cyx, &cyy, &cyz );
    plP_gdom( &xmin, &xmax, &ymin, &ymax );
    plP_grange( &zscale, &zmin_in, &zmax_in );
    zmin = ( zmax_in > zmin_in ) ? zmin_in : zmax_in;
    zmax = ( zmax_in > zmin_in ) ? zmax_in : zmin_in;

    pldtik( zmin, zmax, &tick, &nsub, FALSE );
    tp     = tick * floor( zmin / tick ) + tick;
    pl3upv = 0;

    if ( cxx >= 0.0 && cxy <= 0.0 )
    {
        while ( tp <= zmax )
        {
            u[0] = plP_wcpcx( plP_w3wcx( xmin, ymax, tp ) );
            v[0] = plP_wcpcy( plP_w3wcy( xmin, ymax, tp ) );
            u[1] = plP_wcpcx( plP_w3wcx( xmax, ymax, tp ) );
            v[1] = plP_wcpcy( plP_w3wcy( xmax, ymax, tp ) );
            u[2] = plP_wcpcx( plP_w3wcx( xmax, ymin, tp ) );
            v[2] = plP_wcpcy( plP_w3wcy( xmax, ymin, tp ) );
            plnxtv( u, v, 0, 3, 0 );

            tp += tick;
        }
        u[0] = plP_wcpcx( plP_w3wcx( xmax, ymax, zmin ) );
        v[0] = plP_wcpcy( plP_w3wcy( xmax, ymax, zmin ) );
        u[1] = plP_wcpcx( plP_w3wcx( xmax, ymax, zmax ) );
        v[1] = plP_wcpcy( plP_w3wcy( xmax, ymax, zmax ) );
        plnxtv( u, v, 0, 2, 0 );
    }
    else if ( cxx <= 0.0 && cxy <= 0.0 )
    {
        while ( tp <= zmax )
        {
            u[0] = plP_wcpcx( plP_w3wcx( xmax, ymax, tp ) );
            v[0] = plP_wcpcy( plP_w3wcy( xmax, ymax, tp ) );
            u[1] = plP_wcpcx( plP_w3wcx( xmax, ymin, tp ) );
            v[1] = plP_wcpcy( plP_w3wcy( xmax, ymin, tp ) );
            u[2] = plP_wcpcx( plP_w3wcx( xmin, ymin, tp ) );
            v[2] = plP_wcpcy( plP_w3wcy( xmin, ymin, tp ) );
            plnxtv( u, v, 0, 3, 0 );

            tp += tick;
        }
        u[0] = plP_wcpcx( plP_w3wcx( xmax, ymin, zmin ) );
        v[0] = plP_wcpcy( plP_w3wcy( xmax, ymin, zmin ) );
        u[1] = plP_wcpcx( plP_w3wcx( xmax, ymin, zmax ) );
        v[1] = plP_wcpcy( plP_w3wcy( xmax, ymin, zmax ) );
        plnxtv( u, v, 0, 2, 0 );
    }
    else if ( cxx <= 0.0 && cxy >= 0.0 )
    {
        while ( tp <= zmax )
        {
            u[0] = plP_wcpcx( plP_w3wcx( xmax, ymin, tp ) );
            v[0] = plP_wcpcy( plP_w3wcy( xmax, ymin, tp ) );
            u[1] = plP_wcpcx( plP_w3wcx( xmin, ymin, tp ) );
            v[1] = plP_wcpcy( plP_w3wcy( xmin, ymin, tp ) );
            u[2] = plP_wcpcx( plP_w3wcx( xmin, ymax, tp ) );
            v[2] = plP_wcpcy( plP_w3wcy( xmin, ymax, tp ) );
            plnxtv( u, v, 0, 3, 0 );

            tp += tick;
        }
        u[0] = plP_wcpcx( plP_w3wcx( xmin, ymin, zmin ) );
        v[0] = plP_wcpcy( plP_w3wcy( xmin, ymin, zmin ) );
        u[1] = plP_wcpcx( plP_w3wcx( xmin, ymin, zmax ) );
        v[1] = plP_wcpcy( plP_w3wcy( xmin, ymin, zmax ) );
        plnxtv( u, v, 0, 2, 0 );
    }
    else if ( cxx >= 0.0 && cxy >= 0.0 )
    {
        while ( tp <= zmax )
        {
            u[0] = plP_wcpcx( plP_w3wcx( xmin, ymin, tp ) );
            v[0] = plP_wcpcy( plP_w3wcy( xmin, ymin, tp ) );
            u[1] = plP_wcpcx( plP_w3wcx( xmin, ymax, tp ) );
            v[1] = plP_wcpcy( plP_w3wcy( xmin, ymax, tp ) );
            u[2] = plP_wcpcx( plP_w3wcx( xmax, ymax, tp ) );
            v[2] = plP_wcpcy( plP_w3wcy( xmax, ymax, tp ) );
            plnxtv( u, v, 0, 3, 0 );

            tp += tick;
        }
        u[0] = plP_wcpcx( plP_w3wcx( xmin, ymax, zmin ) );
        v[0] = plP_wcpcy( plP_w3wcy( xmin, ymax, zmin ) );
        u[1] = plP_wcpcx( plP_w3wcx( xmin, ymax, zmax ) );
        v[1] = plP_wcpcy( plP_w3wcy( xmin, ymax, zmax ) );
        plnxtv( u, v, 0, 2, 0 );
    }
    pl3upv = 1;
}

//--------------------------------------------------------------------------
// void plnxtv()
//
// Draw the next view of a 3-d plot. The physical coordinates of the
// points for the next view are placed in the n points of arrays u and
// v. The silhouette found so far is stored in the heap as a set of m peak
// points.
//
// These routines dynamically allocate memory for hidden line removal.
// Memory is allocated in blocks of 2*BINC*sizeof(PLINT) bytes.  Large
// values of BINC give better performance but also allocate more memory
// than is needed. If your 3D plots are very "spiky" or you are working
// with very large matrices then you will probably want to increase BINC.
//--------------------------------------------------------------------------

static void
plnxtv( PLINT *u, PLINT *v, PLFLT* c, PLINT n, PLINT init )
{
    plnxtvhi( u, v, c, n, init );

    if ( pl3mode )
        plnxtvlo( u, v, c, n, init );
}

//--------------------------------------------------------------------------
// void plnxtvhi()
//
// Draw the top side of the 3-d plot.
//--------------------------------------------------------------------------

static void
plnxtvhi( PLINT *u, PLINT *v, PLFLT* c, PLINT n, PLINT init )
{
    //
    // For the initial set of points, just display them and store them as the
    // peak points.
    //
    if ( init == 1 )
    {
        int i;
        oldhiview = (PLINT *) malloc( (size_t) ( 2 * n ) * sizeof ( PLINT ) );
        if ( !oldhiview )
            myexit( "plnxtvhi: Out of memory." );

        oldhiview[0] = u[0];
        oldhiview[1] = v[0];
        plP_draw3d( u[0], v[0], c, 0, 1 );
        for ( i = 1; i < n; i++ )
        {
            oldhiview[2 * i]     = u[i];
            oldhiview[2 * i + 1] = v[i];
            plP_draw3d( u[i], v[i], c, i, 0 );
        }
        mhi = n;
        return;
    }

    //
    // Otherwise, we need to consider hidden-line removal problem. We scan
    // over the points in both the old (i.e. oldhiview[]) and new (i.e. u[]
    // and v[]) arrays in order of increasing x coordinate.  At each stage, we
    // find the line segment in the other array (if one exists) that straddles
    // the x coordinate of the point. We have to determine if the point lies
    // above or below the line segment, and to check if the below/above status
    // has changed since the last point.
    //
    // If pl3upv = 0 we do not update the view, this is useful for drawing
    // lines on the graph after we are done plotting points.  Hidden line
    // removal is still done, but the view is not updated.
    //
    xxhi = 0;
    if ( pl3upv != 0 )
    {
        newhisize = 2 * ( mhi + BINC );
        if ( newhiview != NULL )
        {
            newhiview =
                (PLINT *) realloc( (void *) newhiview,
                    (size_t) newhisize * sizeof ( PLINT ) );
        }
        else
        {
            newhiview =
                (PLINT *) malloc( (size_t) newhisize * sizeof ( PLINT ) );
        }
        if ( !newhiview )
            myexit( "plnxtvhi: Out of memory." );
    }

    // Do the draw or shading with hidden line removal

    plnxtvhi_draw( u, v, c, n );

    // Set oldhiview

    swaphiview();
}

//--------------------------------------------------------------------------
// void plnxtvhi_draw()
//
// Draw the top side of the 3-d plot.
//--------------------------------------------------------------------------

static void
plnxtvhi_draw( PLINT *u, PLINT *v, PLFLT* c, PLINT n )
{
    PLINT i   = 0, j = 0, first = 1;
    PLINT sx1 = 0, sx2 = 0, sy1 = 0, sy2 = 0;
    PLINT su1, su2, sv1, sv2;
    PLINT cx, cy, px, py;
    PLINT seg, ptold, lstold = 0, pthi, pnewhi = 0, newhi, change, ochange = 0;

//
// (oldhiview[2*i], oldhiview[2*i]) is the i'th point in the old array
// (u[j], v[j]) is the j'th point in the new array
//

//
// First attempt at 3d shading.  It works ok for simple plots, but
// will just not draw faces, or draw them overlapping for very
// jagged plots
//

    while ( i < mhi || j < n )
    {
        //
        // The coordinates of the point under consideration are (px,py).  The
        // line segment joins (sx1,sy1) to (sx2,sy2).  "ptold" is true if the
        // point lies in the old array. We set it by comparing the x coordinates
        // of the i'th old point and the j'th new point, being careful if we
        // have fallen past the edges. Having found the point, load up the point
        // and segment coordinates appropriately.
        //

        ptold = ( j >= n || ( i < mhi && oldhiview[2 * i] < u[j] ) );
        if ( ptold )
        {
            px  = oldhiview[2 * i];
            py  = oldhiview[2 * i + 1];
            seg = j > 0 && j < n;
            if ( seg )
            {
                sx1 = u[j - 1];
                sy1 = v[j - 1];
                sx2 = u[j];
                sy2 = v[j];
            }
        }
        else
        {
            px  = u[j];
            py  = v[j];
            seg = i > 0 && i < mhi;
            if ( seg )
            {
                sx1 = oldhiview[2 * ( i - 1 )];
                sy1 = oldhiview[2 * ( i - 1 ) + 1];
                sx2 = oldhiview[2 * i];
                sy2 = oldhiview[2 * i + 1];
            }
        }

        //
        // Now determine if the point is higher than the segment, using the
        // logical function "above". We also need to know if it is the old view
        // or the new view that is higher. "newhi" is set true if the new view
        // is higher than the old.
        //
        if ( seg )
            pthi = plabv( px, py, sx1, sy1, sx2, sy2 );
        else
            pthi = 1;

        newhi = ( ptold && !pthi ) || ( !ptold && pthi );
        //
        // The last point and this point lie on different sides of
        // the current silouette
        //
        change = ( newhi && !pnewhi ) || ( !newhi && pnewhi );

        //
        // There is a new intersection point to put in the peak array if the
        // state of "newhi" changes.
        //
        if ( first )
        {
            plP_draw3d( px, py, c, j, 1 );
            first  = 0;
            lstold = ptold;
            savehipoint( px, py );
            pthi    = 0;
            ochange = 0;
        }
        else if ( change )
        {
            //
            // Take care of special cases at end of arrays.  If pl3upv is 0 the
            // endpoints are not connected to the old view.
            //
            if ( pl3upv == 0 && ( ( !ptold && j == 0 ) || ( ptold && i == 0 ) ) )
            {
                plP_draw3d( px, py, c, j, 1 );
                lstold  = ptold;
                pthi    = 0;
                ochange = 0;
            }
            else if ( pl3upv == 0 &&
                      ( ( !ptold && i >= mhi ) || ( ptold && j >= n ) ) )
            {
                plP_draw3d( px, py, c, j, 1 );
                lstold  = ptold;
                pthi    = 0;
                ochange = 0;
            }
            else
            {
                //
                // If pl3upv is not 0 then we do want to connect the current line
                // with the previous view at the endpoints.  Also find intersection
                // point with old view.
                //
                if ( i == 0 )
                {
                    sx1 = oldhiview[0];
                    sy1 = -1;
                    sx2 = oldhiview[0];
                    sy2 = oldhiview[1];
                }
                else if ( i >= mhi )
                {
                    sx1 = oldhiview[2 * ( mhi - 1 )];
                    sy1 = oldhiview[2 * ( mhi - 1 ) + 1];
                    sx2 = oldhiview[2 * ( mhi - 1 )];
                    sy2 = -1;
                }
                else
                {
                    sx1 = oldhiview[2 * ( i - 1 )];
                    sy1 = oldhiview[2 * ( i - 1 ) + 1];
                    sx2 = oldhiview[2 * i];
                    sy2 = oldhiview[2 * i + 1];
                }

                if ( j == 0 )
                {
                    su1 = u[0];
                    sv1 = -1;
                    su2 = u[0];
                    sv2 = v[0];
                }
                else if ( j >= n )
                {
                    su1 = u[n - 1];
                    sv1 = v[n - 1];
                    su2 = u[n - 1];
                    sv2 = -1;
                }
                else
                {
                    su1 = u[j - 1];
                    sv1 = v[j - 1];
                    su2 = u[j];
                    sv2 = v[j];
                }

                // Determine the intersection

                pl3cut( sx1, sy1, sx2, sy2, su1, sv1, su2, sv2, &cx, &cy );
                if ( cx == px && cy == py )
                {
                    if ( lstold && !ochange )
                        plP_draw3d( px, py, c, j, 1 );
                    else
                        plP_draw3d( px, py, c, j, 0 );

                    savehipoint( px, py );
                    lstold = 1;
                    pthi   = 0;
                }
                else
                {
                    if ( lstold && !ochange )
                        plP_draw3d( cx, cy, c, j, 1 );
                    else
                        plP_draw3d( cx, cy, c, j, 0 );

                    lstold = 1;
                    savehipoint( cx, cy );
                }
                ochange = 1;
            }
        }

        // If point is high then draw plot to point and update view.

        if ( pthi )
        {
            if ( lstold && ptold )
                plP_draw3d( px, py, c, j, 1 );
            else
                plP_draw3d( px, py, c, j, 0 );

            savehipoint( px, py );
            lstold  = ptold;
            ochange = 0;
        }
        pnewhi = newhi;

        if ( ptold )
            i++;
        else
            j++;
    }
}

//--------------------------------------------------------------------------
// void  plP_draw3d()
//
// Does a simple move or line draw.
//--------------------------------------------------------------------------

static void
plP_draw3d( PLINT x, PLINT y, PLFLT *c, PLINT j, PLINT move )
{
    if ( move )
        plP_movphy( x, y );
    else
    {
        if ( c != NULL )
            plcol1( c[j - 1] );
        plP_draphy( x, y );
    }
}

//--------------------------------------------------------------------------
// void plnxtvlo()
//
// Draw the bottom side of the 3-d plot.
//--------------------------------------------------------------------------

static void
plnxtvlo( PLINT *u, PLINT *v, PLFLT*c, PLINT n, PLINT init )
{
    PLINT i, j, first;
    PLINT sx1 = 0, sx2 = 0, sy1 = 0, sy2 = 0;
    PLINT su1, su2, sv1, sv2;
    PLINT cx, cy, px, py;
    PLINT seg, ptold, lstold = 0, ptlo, pnewlo, newlo, change, ochange = 0;

    first  = 1;
    pnewlo = 0;

    //
    // For the initial set of points, just display them and store them as the
    // peak points.
    //
    if ( init == 1 )
    {
        oldloview = (PLINT *) malloc( (size_t) ( 2 * n ) * sizeof ( PLINT ) );
        if ( !oldloview )
            myexit( "\nplnxtvlo: Out of memory." );

        plP_draw3d( u[0], v[0], c, 0, 1 );
        oldloview[0] = u[0];
        oldloview[1] = v[0];
        for ( i = 1; i < n; i++ )
        {
            plP_draw3d( u[i], v[i], c, i, 0 );
            oldloview[2 * i]     = u[i];
            oldloview[2 * i + 1] = v[i];
        }
        mlo = n;
        return;
    }

    //
    // Otherwise, we need to consider hidden-line removal problem. We scan
    // over the points in both the old (i.e. oldloview[]) and new (i.e. u[]
    // and v[]) arrays in order of increasing x coordinate.  At each stage, we
    // find the line segment in the other array (if one exists) that straddles
    // the x coordinate of the point. We have to determine if the point lies
    // above or below the line segment, and to check if the below/above status
    // has changed since the last point.
    //
    // If pl3upv = 0 we do not update the view, this is useful for drawing
    // lines on the graph after we are done plotting points.  Hidden line
    // removal is still done, but the view is not updated.
    //
    xxlo = 0;
    i    = 0;
    j    = 0;
    if ( pl3upv != 0 )
    {
        newlosize = 2 * ( mlo + BINC );
        if ( newloview != NULL )
        {
            newloview =
                (PLINT *) realloc( (void *) newloview,
                    (size_t) newlosize * sizeof ( PLINT ) );
        }
        else
        {
            newloview =
                (PLINT *) malloc( (size_t) newlosize * sizeof ( PLINT ) );
        }
        if ( !newloview )
            myexit( "plnxtvlo: Out of memory." );
    }

    //
    // (oldloview[2*i], oldloview[2*i]) is the i'th point in the old array
    // (u[j], v[j]) is the j'th point in the new array.
    //
    while ( i < mlo || j < n )
    {
        //
        // The coordinates of the point under consideration are (px,py).  The
        // line segment joins (sx1,sy1) to (sx2,sy2).  "ptold" is true if the
        // point lies in the old array. We set it by comparing the x coordinates
        // of the i'th old point and the j'th new point, being careful if we
        // have fallen past the edges. Having found the point, load up the point
        // and segment coordinates appropriately.
        //
        ptold = ( j >= n || ( i < mlo && oldloview[2 * i] < u[j] ) );
        if ( ptold )
        {
            px  = oldloview[2 * i];
            py  = oldloview[2 * i + 1];
            seg = j > 0 && j < n;
            if ( seg )
            {
                sx1 = u[j - 1];
                sy1 = v[j - 1];
                sx2 = u[j];
                sy2 = v[j];
            }
        }
        else
        {
            px  = u[j];
            py  = v[j];
            seg = i > 0 && i < mlo;
            if ( seg )
            {
                sx1 = oldloview[2 * ( i - 1 )];
                sy1 = oldloview[2 * ( i - 1 ) + 1];
                sx2 = oldloview[2 * i];
                sy2 = oldloview[2 * i + 1];
            }
        }

        //
        // Now determine if the point is lower than the segment, using the
        // logical function "above". We also need to know if it is the old view
        // or the new view that is lower. "newlo" is set true if the new view is
        // lower than the old.
        //
        if ( seg )
            ptlo = !plabv( px, py, sx1, sy1, sx2, sy2 );
        else
            ptlo = 1;

        newlo  = ( ptold && !ptlo ) || ( !ptold && ptlo );
        change = ( newlo && !pnewlo ) || ( !newlo && pnewlo );

        //
        // There is a new intersection point to put in the peak array if the
        // state of "newlo" changes.
        //
        if ( first )
        {
            plP_draw3d( px, py, c, j, 1 );
            first  = 0;
            lstold = ptold;
            savelopoint( px, py );
            ptlo    = 0;
            ochange = 0;
        }
        else if ( change )
        {
            //
            // Take care of special cases at end of arrays.  If pl3upv is 0 the
            // endpoints are not connected to the old view.
            //
            if ( pl3upv == 0 && ( ( !ptold && j == 0 ) || ( ptold && i == 0 ) ) )
            {
                plP_draw3d( px, py, c, j, 1 );
                lstold  = ptold;
                ptlo    = 0;
                ochange = 0;
            }
            else if ( pl3upv == 0 &&
                      ( ( !ptold && i >= mlo ) || ( ptold && j >= n ) ) )
            {
                plP_draw3d( px, py, c, j, 1 );
                lstold  = ptold;
                ptlo    = 0;
                ochange = 0;
            }

            //
            // If pl3upv is not 0 then we do want to connect the current line
            // with the previous view at the endpoints.  Also find intersection
            // point with old view.
            //
            else
            {
                if ( i == 0 )
                {
                    sx1 = oldloview[0];
                    sy1 = 100000;
                    sx2 = oldloview[0];
                    sy2 = oldloview[1];
                }
                else if ( i >= mlo )
                {
                    sx1 = oldloview[2 * ( mlo - 1 )];
                    sy1 = oldloview[2 * ( mlo - 1 ) + 1];
                    sx2 = oldloview[2 * ( mlo - 1 )];
                    sy2 = 100000;
                }
                else
                {
                    sx1 = oldloview[2 * ( i - 1 )];
                    sy1 = oldloview[2 * ( i - 1 ) + 1];
                    sx2 = oldloview[2 * i];
                    sy2 = oldloview[2 * i + 1];
                }

                if ( j == 0 )
                {
                    su1 = u[0];
                    sv1 = 100000;
                    su2 = u[0];
                    sv2 = v[0];
                }
                else if ( j >= n )
                {
                    su1 = u[n - 1];
                    sv1 = v[n - 1];
                    su2 = u[n];
                    sv2 = 100000;
                }
                else
                {
                    su1 = u[j - 1];
                    sv1 = v[j - 1];
                    su2 = u[j];
                    sv2 = v[j];
                }

                // Determine the intersection

                pl3cut( sx1, sy1, sx2, sy2, su1, sv1, su2, sv2, &cx, &cy );
                if ( cx == px && cy == py )
                {
                    if ( lstold && !ochange )
                        plP_draw3d( px, py, c, j, 1 );
                    else
                        plP_draw3d( px, py, c, j, 0 );

                    savelopoint( px, py );
                    lstold = 1;
                    ptlo   = 0;
                }
                else
                {
                    if ( lstold && !ochange )
                        plP_draw3d( cx, cy, c, j, 1 );
                    else
                        plP_draw3d( cx, cy, c, j, 0 );

                    lstold = 1;
                    savelopoint( cx, cy );
                }
                ochange = 1;
            }
        }

        // If point is low then draw plot to point and update view.

        if ( ptlo )
        {
            if ( lstold && ptold )
                plP_draw3d( px, py, c, j, 1 );
            else
                plP_draw3d( px, py, c, j, 0 );

            savelopoint( px, py );
            lstold  = ptold;
            ochange = 0;
        }

        pnewlo = newlo;

        if ( ptold )
            i = i + 1;
        else
            j = j + 1;
    }

    // Set oldloview

    swaploview();
}

//--------------------------------------------------------------------------
// savehipoint
// savelopoint
//
// Add a point to the list of currently visible peaks/valleys, when
// drawing the top/bottom surface, respectively.
//--------------------------------------------------------------------------

static void
savehipoint( PLINT px, PLINT py )
{
    if ( pl3upv == 0 )
        return;

    if ( xxhi >= newhisize )      // allocate additional space
    {
        newhisize += 2 * BINC;
        newhiview  = (PLINT *) realloc( (void *) newhiview,
            (size_t) newhisize * sizeof ( PLINT ) );
        if ( !newhiview )
            myexit( "savehipoint: Out of memory." );
    }

    newhiview[xxhi] = px;
    xxhi++;
    newhiview[xxhi] = py;
    xxhi++;
}

static void
savelopoint( PLINT px, PLINT py )
{
    if ( pl3upv == 0 )
        return;

    if ( xxlo >= newlosize )      // allocate additional space
    {
        newlosize += 2 * BINC;
        newloview  = (PLINT *) realloc( (void *) newloview,
            (size_t) newlosize * sizeof ( PLINT ) );
        if ( !newloview )
            myexit( "savelopoint: Out of memory." );
    }

    newloview[xxlo] = px;
    xxlo++;
    newloview[xxlo] = py;
    xxlo++;
}

//--------------------------------------------------------------------------
// swaphiview
// swaploview
//
// Swaps the top/bottom views.  Need to do a real swap so that the
// memory cleanup routine really frees everything (and only once).
//--------------------------------------------------------------------------

static void
swaphiview( void )
{
    PLINT *tmp;

    if ( pl3upv != 0 )
    {
        mhi       = xxhi / 2;
        tmp       = oldhiview;
        oldhiview = newhiview;
        newhiview = tmp;
    }
}

static void
swaploview( void )
{
    PLINT *tmp;

    if ( pl3upv != 0 )
    {
        mlo       = xxlo / 2;
        tmp       = oldloview;
        oldloview = newloview;
        newloview = tmp;
    }
}

//--------------------------------------------------------------------------
// freework
//
// Frees memory associated with work arrays
//--------------------------------------------------------------------------

static void
freework( void )
{
    free_mem( oldhiview );
    free_mem( oldloview );
    free_mem( newhiview );
    free_mem( newloview );
    free_mem( vtmp );
    free_mem( utmp );
    free_mem( ctmp );
}

//--------------------------------------------------------------------------
// myexit
//
// Calls plexit, cleaning up first.
//--------------------------------------------------------------------------

static void
myexit( PLCHAR_VECTOR msg )
{
    freework();
    plexit( msg );
}

//--------------------------------------------------------------------------
// myabort
//
// Calls plabort, cleaning up first.
// Caller should return to the user program.
//--------------------------------------------------------------------------

static void
myabort( PLCHAR_VECTOR msg )
{
    freework();
    plabort( msg );
}

//--------------------------------------------------------------------------
// int plabv()
//
// Determines if point (px,py) lies above the line joining (sx1,sy1) to
// (sx2,sy2). It only works correctly if sx1 <= px <= sx2.
//--------------------------------------------------------------------------

static int
plabv( PLINT px, PLINT py, PLINT sx1, PLINT sy1, PLINT sx2, PLINT sy2 )
{
    int above;

    if ( py >= sy1 && py >= sy2 )
        above = 1;
    else if ( py < sy1 && py < sy2 )
        above = 0;
    else if ( (double) ( sx2 - sx1 ) * ( py - sy1 ) >=
              (double) ( px - sx1 ) * ( sy2 - sy1 ) )
        above = 1;
    else
        above = 0;

    return above;
}

//--------------------------------------------------------------------------
// void pl3cut()
//
// Determines the point of intersection (cx,cy) between the line joining
// (sx1,sy1) to (sx2,sy2) and the line joining (su1,sv1) to (su2,sv2).
//--------------------------------------------------------------------------

static void
pl3cut( PLINT sx1, PLINT sy1, PLINT sx2, PLINT sy2,
        PLINT su1, PLINT sv1, PLINT su2, PLINT sv2, PLINT *cx, PLINT *cy )
{
    PLINT x21, y21, u21, v21, yv1, xu1, a, b;
    double fa, fb;

    x21 = sx2 - sx1;
    y21 = sy2 - sy1;
    u21 = su2 - su1;
    v21 = sv2 - sv1;
    yv1 = sy1 - sv1;
    xu1 = sx1 - su1;

    a  = x21 * v21 - y21 * u21;
    fa = (double) a;
    if ( a == 0 )
    {
        if ( sx2 < su2 )
        {
            *cx = sx2;
            *cy = sy2;
        }
        else
        {
            *cx = su2;
            *cy = sv2;
        }
        return;
    }
    else
    {
        b   = yv1 * u21 - xu1 * v21;
        fb  = (double) b;
        *cx = (PLINT) ( sx1 + ( fb * x21 ) / fa + .5 );
        *cy = (PLINT) ( sy1 + ( fb * y21 ) / fa + .5 );
    }
}

//--------------------------------------------------------------------------
// void plRotationShear
//
// Calculates the rotation and shear angles from a plplot transformation matrix
//
// N.B. the plot transformation matrix
//
// [xFormMatrix[0] xFormMatrix[2]]
// [xFormMatrix[1] xFormMatrix[3]]
//
// is calculated as
//
// [stride cos(t)    stride sin(t)]
// [sin(p-t)              cos(p-t)]
//
// where t is the rotation angle and p is the shear angle.
// The purpose of this routine is to determine stride, rotation angle,
// and shear angle from xFormMatrix.
//
// For information only, xFormMatrix is the product of the following
// rotation, skew(shear), and scale matrices:
//
//  [stride    0] [1      0] [cos(t)  sin(t)]
//  [0    cos(p)] [tan(p) 1] [-sin(t) cos(t)]
//
//--------------------------------------------------------------------------

void
plRotationShear( PLFLT *xFormMatrix, PLFLT *rotation, PLFLT *shear, PLFLT *stride )
{
    PLFLT smr;
    *stride = sqrt( xFormMatrix[0] * xFormMatrix[0] + xFormMatrix[2] * xFormMatrix[2] );

    // Calculate rotation in range from -pi to pi.
    *rotation = atan2( xFormMatrix[2], xFormMatrix[0] );

    // Calculate shear - rotation in range from -pi to pi.
    smr = atan2( xFormMatrix[1], xFormMatrix[3] );

    // Calculate shear in range from -2 pi to 2 pi.
    *shear = smr + *rotation;

    // Calculate shear in range from -pi to pi.
    if ( *shear < -PI )
        *shear += 2. * PI;
    else if ( *shear > PI )
        *shear -= 2. * PI;

    // Actually must honour some convention to calculate the negative
    // of the shear angle instead of the shear angle. Why??
    *shear = -*shear;
    // Comment out the modified old logic which determines the negative
    // of the shear angle in a more complicated way.  Note, the modification
    // to divide the asin argument by *stride which solved a long-standing
    // bug (as does the above logic in a simpler way).
    //
    //shear = -asin( (xFormMatrix[0] * xFormMatrix[1] +
    //               xFormMatrix[2] * xFormMatrix[3] )/ *stride);
    //

    // Compute the cross product of the vectors [1,0] and [0,1] to
    // determine if we need to make a "quadrant 3,4" adjustment
    // to the shear angle.

    //
    // if ( xFormMatrix[0] * xFormMatrix[3] - xFormMatrix[1] * xFormMatrix[2] < 0.0 )
    // {
    //shear = -( M_PI + *shear );
    // }
    //
}