/* quadric.f -- translated by f2c (version 20061008).
   You must link the resulting object file with libf2c:
	on Microsoft Windows system, link with libf2c.lib;
	on Linux or Unix systems, link with .../path/to/libf2c.a -lm
	or, if you install libf2c.a in a standard place, with -lf2c -lm
	-- in that order, at the end of the command line, as in
		cc *.o -lf2c -lm
	Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,

		http://www.netlib.org/f2c/libf2c.zip
*/

#include "f2c.h"

/* Common Block Declarations */

struct {
    real xcent, ycent, scale, eyepos, sxcent, sycent, tmat[16]	/* was [4][4] 
	    */, tinv[16]	/* was [4][4] */, tinvt[16]	/* was [4][4] 
	    */, srot[16]	/* was [4][4] */, srtinv[16]	/* was [4][4] 
	    */, srtinvt[16]	/* was [4][4] */, rafter[16]	/* was [4][4] 
	    */, tafter[3];
} matrices_;

#define matrices_1 matrices_

struct {
    integer noise;
    logical verbose;
} asscom_;

#define asscom_1 asscom_

struct {
    integer ntx, nty, npx, npy;
} raster_;

#define raster_1 raster_

struct {
    integer kount[65536]	/* was [256][256] */, mount[129600]	/* 
	    was [360][360] */, ttrans[65536]	/* was [256][256] */, istrans;
} lists_;

#define lists_1 lists_

struct {
    real trulim[6]	/* was [3][2] */, zlim[2], frontclip, backclip;
    logical isolation;
} niceties_;

#define niceties_1 niceties_

/* Table of constant values */

static integer c__1 = 1;
static integer c__3 = 3;
static integer c__4 = 4;
static logical c_false = FALSE_;

/* *********************************************************************** */
/*              Support routines for quadric surfaces                   * */
/* *********************************************************************** */
/* EAM Jun 1997	- initial version, supports version 2.4(alpha) of render */
/* EAM May 1998	- additional error checking to go with Parvati/rastep */
/* EAM Jan 1999	- version 2.4i */
/* *********************************************************************** */

/* Quadric surfaces include spheres, cones, ellipsoids, paraboloids, and */
/* hyperboloids.  The motivation for this code was to allow rendering */
/* thermal ellipsoids for atoms, so the other shapes have not been */
/* extensively tested. */
/* A quadric surface is described by 10 parameters (A ... J). */
/* For efficiency during rendering it is also useful to know the center and */
/* a bounding sphere. So a QUADRIC descriptor to render has 17 parameters: */
/* 14 (object type QUADRIC) */
/*    X Y Z RADLIM RED GRN BLU */
/*    A B C D E F G H I J */

/* The surface itself is the set of points for which Q(x,y,z) = 0 */
/* where */
/*      Q(x,y,z)	=  A*x^2 +  B*y^2 +  C*z^2 */
/*         	+ 2D*x*y + 2E*y*z + 2F*z*x */
/*         	+ 2G*x   + 2H*y   + 2I*z */
/*         	+  J */

/* It is convenient to store this information in a matrix QQ */
/* 	| QA QD QF QG |		QA = A    QB = B    QC = C */
/* QQ =	| QD QB QE QH |		QD = D    QE = E    QF = F */
/* 	| QF QE QC QI |		QG = G    QH = H    QI = I */
/* 	| QG QH QI QJ |		QJ = J */

/* Then Q(x,y,x) = XT*QQ*X	where X = (x,y,z,1) */
/* The point of this is that a 4x4 homogeneous transformation T can be */
/* applied to QQ by matrix multiplication:     QQ' = TinvT * QQ * Tinv */

/* The surface normal is easily found by taking the partial derivatives */
/* of Q(x,y,z) at the point of interest. */
/* *********************************************************************** */
/* TO DO: */
/* 	- can we distinguish an ellipsoid from other quadrics? Do we care? */
/* *********************************************************************** */
/* CC	Calculate the matrices for quadric transformation */
/* C	QQ' = TINVT * QQ * TINV      (TINV is inverse of transposed TMAT) */

/* Subroutine */ int qsetup_(void)
{
    /* Format strings */
    static char fmt_901[] = "(\002>>> Warning: Determinant of rotation matri"
	    "x =\002,f7.3,\002 <<<\002)";
    static char fmt_902[] = "(\002>>> Warning: Non-zero cross terms in 4th c"
	    "olumn of\002,\002 TMAT <<<\002)";
    static char fmt_903[] = "(\002>>> Warning: Determinant of shadow matri"
	    "x =\002,f7.3,\002 <<<\002)";

    /* System generated locals */
    real r__1;

    /* Builtin functions */
    integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void);

    /* Local variables */
    static real det, tran[16]	/* was [4][4] */, post[16]	/* was [4][4] 
	    */;
    extern doublereal tinv4_(real *, real *);
    extern /* Subroutine */ int tmul4_(real *, real *, real *), trnsp4_(real *
	    , real *);

    /* Fortran I/O blocks */
    static cilist io___4 = { 0, 0, 0, fmt_901, 0 };
    static cilist io___5 = { 0, 0, 0, fmt_902, 0 };
    static cilist io___6 = { 0, 0, 0, fmt_903, 0 };






/* 	TMAT is a post-multiplier matrix, but unfortunately the */
/* 	quadric surface math was worked out for pre-multipliers */
/* 	So the quadric math uses the transpose of TMAT */
    trnsp4_(tran, matrices_1.tmat);

/* 	Remove translation component from TMAT. */
/* 	I had to make a choice whether the X,Y,Z "center" of the quadric */
/* 	is implicitly coded into the coefficients, or whether it is to be */
/* 	maintained explicitly.  Since all other object types do the latter, */
/* 	I have chosen to keep all translation components out of the quadric */
/* 	coefficients */
    tran[12] = 0.f;
    tran[13] = 0.f;
    tran[14] = 0.f;

/* 	July 1999 - Allow post-hoc rotation matrix also */
    tmul4_(post, matrices_1.rafter, tran);

    det = tinv4_(matrices_1.tinv, post) / matrices_1.tmat[15];
    trnsp4_(matrices_1.tinvt, matrices_1.tinv);

/* 	While we're at it, check for legality of rotation matrix */
    if ((r__1 = 1.f - dabs(det), dabs(r__1)) > .02f) {
	io___4.ciunit = asscom_1.noise;
	s_wsfe(&io___4);
	r__1 = dabs(det);
	do_fio(&c__1, (char *)&r__1, (ftnlen)sizeof(real));
	e_wsfe();
    }
    if (matrices_1.tmat[12] != 0.f || matrices_1.tmat[13] != 0.f || 
	    matrices_1.tmat[14] != 0.f) {
	io___5.ciunit = asscom_1.noise;
	s_wsfe(&io___5);
	e_wsfe();
    }

/* 	Same thing again for shadow rotation matrix */
    trnsp4_(tran, matrices_1.srot);
    det = tinv4_(matrices_1.srtinv, tran);
    trnsp4_(matrices_1.srtinvt, matrices_1.srtinv);
    if ((r__1 = 1.f - dabs(det), dabs(r__1)) > .02f) {
	io___6.ciunit = asscom_1.noise;
	s_wsfe(&io___6);
	r__1 = dabs(det);
	do_fio(&c__1, (char *)&r__1, (ftnlen)sizeof(real));
	e_wsfe();
    }

    return 0;
} /* qsetup_ */

/* CC     Process single object descriptor during input phase */
/* C */

logical qinp_(real *buf, real *detail, logical *shadow, real *sdtail)
{
    /* System generated locals */
    integer i__1, i__2;
    logical ret_val, L__1;

    /* Local variables */
    static real rc, xc, yc, zc, qp[16]	/* was [4][4] */, qq[16]	/* 
	    was [4][4] */;
    static integer ix, iy;
    static real qt[16]	/* was [4][4] */, rs, xq, yq, zq, xr, yr, zr, xs, ys, 
	    zs, red, blu, grn, pfac;
    static integer ixhi, iyhi, ixlo, iylo;
    extern /* Subroutine */ int tmul4_(real *, real *, real *);
    extern doublereal persp_(real *);
    static real radlim;
    extern /* Subroutine */ int transf_(real *, real *, real *), assert_(
	    logical *, char *, ftnlen);




/*     Array sizes */

/*     $$$$$$$$$$$$$   ARRAY SIZE LIMITS START HERE   $$$$$$$$$$$$$$ */

/*     Maximum number of tiles */

/*     Number of shadow tiles */
/* **   (One of these can fail to be enough when the aspect ratio is */
/* **   extreme or when the model is far from being "centred" near z=0. */
/* **   Keep them well ahead of MAXNTX, MAXNTY to be on the safe side) */
/* **   EAM - Allow soft failure and monitor required values in NSXMAX,NSYMAX */

/*     Maximum number of pixels per tile */

/*     Maximum number of objects */
/* **   PARAMETER (MAXOBJ =   7500) */

/*     Array elements available for object details */
/*     Should be roughly 10*MAXOBJ */
/* **   PARAMETER (MAXDET =   150 000, MAXSDT =   150 000) */
/* **   PARAMETER (MAXDET = 2 000 000, MAXSDT = 2 000 000) */

/*     Array elements available for sorted lists ("short" lists) */
/*     Increased requirements as more objects are stacked behind each other */
/* **   PARAMETER (MAXSHR =   150 000, MAXSSL =   150 000) */

/*     Maximum number of MATERIAL definitions (object type 8) */

/*     Maximum number of stacked transparent objects at any single pixel */
/*     (any further further stacking is ignored) */

/*     Maximum number of non-shadowing lights (object type 13) */

/*     Maximum levels of file indirection in input stream */

/*     $$$$$$$$$$$$$$$$$  END OF LIMITS  $$$$$$$$$$$$$$$$$$$$$$$ */

/*     Other possibly platform-dependent stuff */
/*     Slop is related to the accuracy (in pixels) to which we must predict */
/*     shadow edges. Too low a value causes whole triangles to be spuriously */
/*     in shadow; too high a value may cause shadows to be missed altogether. */
/*     Perfect accuracy in floating point calculations would allow SLOP << 1. */
/*     Edgeslop is similarly a kludge for dealing with triangles whose explicit */
/*     normals describe wrapping around from front-facing to back-facing. */
/*     Ribbonslop is a kludge so that distortion due to perspective doesn't */
/*     prevent us from identifying ribbon triangles */





/* Assume this is legitimate */
    /* Parameter adjustments */
    --sdtail;
    --detail;
    --buf;

    /* Function Body */
    ret_val = TRUE_;

/* Update limits (have to trust the center coords) */

    xq = buf[1];
    yq = buf[2];
    zq = buf[3];
    radlim = buf[4];
    L__1 = radlim >= 0.f;
    assert_(&L__1, "limiting radius < 0 in quadric", (ftnlen)30);
    if (radlim > 0.f) {
	niceties_1.trulim[0] = dmin(niceties_1.trulim[0],xq);
	niceties_1.trulim[3] = dmax(niceties_1.trulim[3],xq);
	niceties_1.trulim[1] = dmin(niceties_1.trulim[1],yq);
	niceties_1.trulim[4] = dmax(niceties_1.trulim[4],yq);
	niceties_1.trulim[2] = dmin(niceties_1.trulim[2],zq);
	niceties_1.trulim[5] = dmax(niceties_1.trulim[5],zq);
    }

/* Standard color checks */

    red = buf[5];
    grn = buf[6];
    blu = buf[7];
    L__1 = red >= 0.f;
    assert_(&L__1, "red < 0 in quadric", (ftnlen)18);
    L__1 = red <= 1.f;
    assert_(&L__1, "red > 1 in quadric", (ftnlen)18);
    L__1 = grn >= 0.f;
    assert_(&L__1, "grn < 0 in quadric", (ftnlen)18);
    L__1 = grn <= 1.f;
    assert_(&L__1, "grn > 1 in quadric", (ftnlen)18);
    L__1 = blu >= 0.f;
    assert_(&L__1, "blu < 0 in quadric", (ftnlen)18);
    L__1 = blu <= 1.f;
    assert_(&L__1, "blu > 1 in quadric", (ftnlen)18);

/* Transform center before saving */
/* (But can we deal with perspective????) */

    transf_(&xq, &yq, &zq);
    radlim /= matrices_1.tmat[15];
    if (matrices_1.eyepos > 0.f) {
/* 	    pfac = 1./(1.-zq/eyepos) */
	pfac = persp_(&zq);
	xq *= pfac;
	yq *= pfac;
	zq *= pfac;
	radlim *= pfac;
    }
    xc = xq * matrices_1.scale + matrices_1.xcent;
    yc = yq * matrices_1.scale + matrices_1.ycent;
    zc = zq * matrices_1.scale;
    rc = radlim * matrices_1.scale;
/* 	save transformed Z limits */
    niceties_1.zlim[0] = dmin(niceties_1.zlim[0],zc);
    niceties_1.zlim[1] = dmax(niceties_1.zlim[1],zc);

/* 	check for Z-clipping */
    if (zc > niceties_1.frontclip || zc < niceties_1.backclip) {
	ret_val = FALSE_;
	return ret_val;
    }

    detail[1] = xc;
    detail[2] = yc;
    detail[3] = zc;
    detail[4] = rc;
    detail[5] = red;
    detail[6] = grn;
    detail[7] = blu;

/* This is a terrible kludge, but necessary if called from normal3d -size BIGxBIG */
/* Thes test should really be if we are called from normal3d but no flag for that */
    if (raster_1.ntx > 256 || raster_1.nty > 256) {
	goto L101;
    }

/* Tally for tiles the object might impinge on */
/* Again we are relying on the correctness of the center coordinates */

    ixlo = (xc - rc) / raster_1.npx + 1;
    ixhi = (xc + rc) / raster_1.npx + 1;
    iylo = (yc - rc) / raster_1.npy + 1;
    iyhi = (yc + rc) / raster_1.npy + 1;
    if (ixlo < 1) {
	ixlo = 1;
    }
    if (ixlo > raster_1.ntx) {
	goto L101;
    }
    if (ixhi < 1) {
	goto L101;
    }
    if (ixhi > raster_1.ntx) {
	ixhi = raster_1.ntx;
    }
    if (iylo < 1) {
	iylo = 1;
    }
    if (iylo > raster_1.nty) {
	goto L101;
    }
    if (iyhi < 1) {
	goto L101;
    }
    if (iyhi > raster_1.nty) {
	iyhi = raster_1.nty;
    }
    i__1 = iyhi;
    for (iy = iylo; iy <= i__1; ++iy) {
	i__2 = ixhi;
	for (ix = ixlo; ix <= i__2; ++ix) {
	    ++lists_1.kount[ix + (iy << 8) - 257];
	    lists_1.ttrans[ix + (iy << 8) - 257] += lists_1.istrans;
	}
    }
L101:

/* build matrix from coeffients describing quadric surface in standard form */

    qq[0] = buf[8];
    qq[5] = buf[9];
    qq[10] = buf[10];
    qq[4] = buf[11];
    qq[1] = buf[11];
    qq[9] = buf[12];
    qq[6] = buf[12];
    qq[8] = buf[13];
    qq[2] = buf[13];
    qq[12] = buf[14];
    qq[3] = buf[14];
    qq[13] = buf[15];
    qq[7] = buf[15];
    qq[14] = buf[16];
    qq[11] = buf[16];
    qq[15] = buf[17];

/* Transformed matrix QP = TINV(Transpose) * QQ * TINV */
/* where TINV is the inverse of TMAT */

    tmul4_(qt, qq, matrices_1.tinv);
    tmul4_(qp, matrices_1.tinvt, qt);
/* D	noise = 0 */
/* D	write (noise,191) 'QT  ',((QT(i,j),j=1,4),i=1,4) */
/* D 191	format(a,4(/,4f8.4)) */
/* D	write (noise,191) 'QP  ',((QP(i,j),j=1,4),i=1,4) */

/* Save components of quadric surface built from transformed matrix */
/* for use during rendering */

    detail[8] = qp[0];
    detail[9] = qp[5];
    detail[10] = qp[10];
    detail[11] = qp[4];
    detail[12] = qp[9];
    detail[13] = qp[8];
    detail[14] = qp[12];
    detail[15] = qp[13];
    detail[16] = qp[14];
    detail[17] = qp[15];

/* Do it all over again for the shadow buffers NOT TESTED YET! */
/* (since I'm a little confused about what transformation */
/* I need to apply to the QQ matrix in shadow space) */

    if (! (*shadow)) {
	return ret_val;
    }
/* 	first transform center and limiting sphere */
    xr = matrices_1.srot[0] * xq + matrices_1.srot[4] * yq + matrices_1.srot[
	    8] * zq;
    yr = matrices_1.srot[1] * xq + matrices_1.srot[5] * yq + matrices_1.srot[
	    9] * zq;
    zr = matrices_1.srot[2] * xq + matrices_1.srot[6] * yq + matrices_1.srot[
	    10] * zq;
    xs = xr * matrices_1.scale + matrices_1.sxcent;
    ys = yr * matrices_1.scale + matrices_1.sycent;
    zs = zr * matrices_1.scale;
    rs = radlim * matrices_1.scale;
    sdtail[1] = xs;
    sdtail[2] = ys;
    sdtail[3] = zs;
    sdtail[4] = rs;
/* 	tally shadow tiles the object might impinge on */
    ixlo = (xs - rs) / raster_1.npx + 1;
    ixhi = (xs + rs) / raster_1.npx + 1;
    iylo = (ys - rs) / raster_1.npy + 1;
    iyhi = (ys + rs) / raster_1.npy + 1;
    if (ixlo < 1) {
	ixlo = 1;
    }
    if (ixlo > 360) {
	goto L209;
    }
    if (ixhi < 1) {
	goto L209;
    }
    if (ixhi > 360) {
	ixhi = 360;
    }
    if (iylo < 1) {
	iylo = 1;
    }
    if (iylo > 360) {
	goto L209;
    }
    if (iyhi < 1) {
	goto L209;
    }
    if (iyhi > 360) {
	iyhi = 360;
    }
    i__1 = iyhi;
    for (iy = iylo; iy <= i__1; ++iy) {
	i__2 = ixhi;
	for (ix = ixlo; ix <= i__2; ++ix) {
	    ++lists_1.mount[ix + iy * 360 - 361];
	}
    }
/* 	transform QQ into shadow space as well */
    tmul4_(qt, qp, matrices_1.srtinv);
    tmul4_(qp, matrices_1.srtinvt, qt);
    sdtail[5] = qp[0];
    sdtail[6] = qp[5];
    sdtail[7] = qp[10];
    sdtail[8] = qp[4];
    sdtail[9] = qp[9];
    sdtail[10] = qp[8];
    sdtail[11] = qp[12];
    sdtail[12] = qp[13];
    sdtail[13] = qp[14];
    sdtail[14] = qp[15];

L209:

/* Done with input */

    return ret_val;
} /* qinp_ */

/* *********************************************************************** */
/* This routine does the exact test for pixel impingement of a quadric  * */
/* in both pixel and shadow space.                                      * */
/* Find Z coordinate of point on quadric surface with given X, Y        * */
/* Also return surface normal at that point                             * */
/* *********************************************************************** */
/* CC */
/* C */

logical qtest_(real *center, real *coeffs, real *xp, real *yp, real *zp, real 
	*qnorm, logical *shadowspace, logical *backside)
{
    /* System generated locals */
    logical ret_val;

    /* Builtin functions */
    double sqrt(doublereal);

    /* Local variables */
    static real x, y, z__, aa, bb, cc, dd, ee, qa, qb, qc, qd, qe, qf, qg, qh,
	     qi, qj;




    /* Parameter adjustments */
    --qnorm;
    --coeffs;
    --center;

    /* Function Body */
    qa = coeffs[1];
    qb = coeffs[2];
    qc = coeffs[3];
    qd = coeffs[4];
    qe = coeffs[5];
    qf = coeffs[6];
    qg = coeffs[7];
    qh = coeffs[8];
    qi = coeffs[9];
    qj = coeffs[10];

/* 	xp and yp are in pixel coordinates, rather than in normalized ones */
/* 	First displace center of quadric surface from the implicit origin */
/* 	to the point specified by center(3), then convert to pixel coords */

    x = (*xp - center[1]) / matrices_1.scale;
    y = (*yp - center[2]) / matrices_1.scale;

    aa = qc;
    bb = (qf * x + qe * y + qi) * 2.f;
    cc = qa * x * x + qb * y * y + (qd * x * y + qg * x + qh * y) * 2.f + qj;
    dd = bb * bb - aa * 4 * cc;
    if (dd < 0.f) {
	ret_val = FALSE_;
	return ret_val;
    } else {
	ret_val = TRUE_;
	ee = sqrt(dd);
    }
    if (aa == 0.f) {
/* DEBUG	    Can this happen???? */
	z__ = 9999.f;
    } else if (*backside) {
	if (aa <= 0.f) {
	    z__ = (-bb + ee) / (aa * 2);
	} else {
	    z__ = (-bb - ee) / (aa * 2);
	}
    } else {
	if (aa > 0.f) {
	    z__ = (-bb + ee) / (aa * 2);
	} else {
	    z__ = (-bb - ee) / (aa * 2);
	}
    }
    *zp = z__ * matrices_1.scale;
    *zp += center[3];

/* 	Surface normal comes from partial derivatives at this point */
    if (! (*shadowspace)) {
	qnorm[1] = qa * x + qd * y + qf * z__ + qg;
	qnorm[2] = qb * y + qd * x + qe * z__ + qh;
	qnorm[3] = qc * z__ + qf * x + qe * y + qi;
    }

    return ret_val;
} /* qtest_ */

/* CC	Convert ANISOU description of anisotropic displacement parameters */
/* C	into quadric surface enclosing a given probability volume */
/* 	Returns -1 if the Uij are non-positive definite, 1 otherwise */

integer anitoquad_(real *anisou, real *prob, real *quadric, real *eigens, 
	real *evecs)
{
    /* System generated locals */
    integer ret_val;
    real r__1;

    /* Builtin functions */
    double sqrt(doublereal);

    /* Local variables */
    static integer i__;
    static real uu[16]	/* was [4][4] */, det, uinv[16]	/* was [4][4] */;
    extern doublereal tinv4_(real *, real *);
    extern /* Subroutine */ int jacobi_(real *, integer *, integer *, real *, 
	    real *);





/* 	Save FPU traps later by checking that ANISOU is non-zero */
    /* Parameter adjustments */
    evecs -= 5;
    --eigens;
    --quadric;
    --anisou;

    /* Function Body */
    for (i__ = 1; i__ <= 3; ++i__) {
	eigens[i__] = 0.f;
    }
    if (anisou[1] == 0.f && anisou[2] == 0.f && anisou[3] == 0.f && anisou[4] 
	    == 0.f && anisou[5] == 0.f && anisou[6] == 0.f) {
	ret_val = -2;
	return ret_val;
    }

/* 	Build matrix from Uij coefficients and invert it */
    uu[0] = anisou[1];
    uu[5] = anisou[2];
    uu[10] = anisou[3];
/* Computing 2nd power */
    r__1 = *prob;
    uu[15] = -(1 / (r__1 * r__1));
    uu[4] = anisou[4];
    uu[1] = anisou[4];
    uu[9] = anisou[6];
    uu[6] = anisou[6];
    uu[8] = anisou[5];
    uu[2] = anisou[5];
    uu[12] = 0.f;
    uu[3] = 0.f;
    uu[13] = 0.f;
    uu[7] = 0.f;
    uu[14] = 0.f;
    uu[11] = 0.f;
    det = tinv4_(uinv, uu);

/* 	and from that we extract the coefficients of the surface */
    quadric[1] = uinv[0];
    quadric[2] = uinv[5];
    quadric[3] = uinv[10];
    quadric[4] = uinv[4];
    quadric[5] = uinv[9];
    quadric[6] = uinv[8];
    quadric[7] = uinv[12];
    quadric[8] = uinv[13];
    quadric[9] = uinv[14];
    quadric[10] = uinv[15];

/* 	Find eigenvalues of the ellipsoid */
    jacobi_(uu, &c__3, &c__4, &eigens[1], &evecs[5]);

/* 	Units of this matrix are A**2; we want values in A */
    ret_val = 1;
    for (i__ = 1; i__ <= 3; ++i__) {
	if (eigens[i__] > 0.f) {
	    eigens[i__] = sqrt(eigens[i__]);
	} else {
/* 		write (noise,*) 'Non-positive definite ellipsoid!' */
	    eigens[i__] = 0.f;
	    ret_val = -1;
	}
    }

    return ret_val;
} /* anitoquad_ */

/* CC	Matrix manipulation routines for 4x4 homogeneous transformations */
/* C */

/* Subroutine */ int tmul4_(real *c__, real *a, real *b)
{
    static integer i__, j;

    /* Parameter adjustments */
    b -= 5;
    a -= 5;
    c__ -= 5;

    /* Function Body */
    for (i__ = 1; i__ <= 4; ++i__) {
	for (j = 1; j <= 4; ++j) {
	    c__[i__ + (j << 2)] = a[i__ + 4] * b[(j << 2) + 1] + a[i__ + 8] * 
		    b[(j << 2) + 2] + a[i__ + 12] * b[(j << 2) + 3] + a[i__ + 
		    16] * b[(j << 2) + 4];
	}
    }
    return 0;
} /* tmul4_ */

/* Subroutine */ int trnsp4_(real *b, real *a)
{
    static integer i__, j;

    /* Parameter adjustments */
    a -= 5;
    b -= 5;

    /* Function Body */
    for (i__ = 1; i__ <= 4; ++i__) {
	for (j = 1; j <= 4; ++j) {
	    b[i__ + (j << 2)] = a[j + (i__ << 2)];
	}
    }
    return 0;
} /* trnsp4_ */

/* CC	Matrix inversion for a 4x4 matrix A */
/* C */

doublereal tinv4_(real *ai, real *a)
{
    /* System generated locals */
    real ret_val;

    /* Local variables */
    static real d__;
    static integer i__, j;
    static real tmp[16]	/* was [4][4] */;
    static integer index[4];
    extern /* Subroutine */ int lubksb_(real *, integer *, integer *, real *),
	     ludcmp_(real *, integer *, integer *, real *);



    /* Parameter adjustments */
    a -= 5;
    ai -= 5;

    /* Function Body */
    for (i__ = 1; i__ <= 4; ++i__) {
	for (j = 1; j <= 4; ++j) {
	    tmp[i__ + (j << 2) - 5] = a[i__ + (j << 2)];
	    ai[i__ + (j << 2)] = 0.f;
	}
	ai[i__ + (i__ << 2)] = 1.f;
    }
    ludcmp_(tmp, &c__4, index, &d__);
    ret_val = d__ * tmp[0] * tmp[5] * tmp[10] * tmp[15];
    for (j = 1; j <= 4; ++j) {
	lubksb_(tmp, &c__4, index, &ai[(j << 2) + 1]);
    }
    return ret_val;
} /* tinv4_ */

/* *********************************************************************** */
/*              Matrix inversion via LU decomposition                   * */
/* 	adapted from Numerical Recipes in Fortran (1986)               * */
/* *********************************************************************** */

/* CC	input  NxN matrix A is replaced by its LU decomposition */
/* C	output index(N) records row permutation due to pivoting */
/* 	output D is parity of row permutations */

/* Subroutine */ int ludcmp_(real *a, integer *n, integer *index, real *d__)
{
    /* System generated locals */
    integer a_dim1, a_offset, i__1, i__2, i__3;
    real r__1, r__2;
    logical L__1;

    /* Local variables */
    static integer i__, j, k;
    static real vv[10], dum, sum;
    static integer imax;
    static real aamax;
    extern /* Subroutine */ int assert_(logical *, char *, ftnlen);



    /* Parameter adjustments */
    --index;
    a_dim1 = *n;
    a_offset = 1 + a_dim1;
    a -= a_offset;

    /* Function Body */
    *d__ = 1.f;
    i__1 = *n;
    for (i__ = 1; i__ <= i__1; ++i__) {
	aamax = 0.f;
	i__2 = *n;
	for (j = 1; j <= i__2; ++j) {
	    if ((r__1 = a[i__ + j * a_dim1], dabs(r__1)) > aamax) {
		aamax = (r__2 = a[i__ + j * a_dim1], dabs(r__2));
	    }
	}
	L__1 = aamax != 0.f;
	assert_(&L__1, "Singular matrix", (ftnlen)15);
	vv[i__ - 1] = 1.f / aamax;
    }

    i__1 = *n;
    for (j = 1; j <= i__1; ++j) {
	if (j > 1) {
	    i__2 = j - 1;
	    for (i__ = 1; i__ <= i__2; ++i__) {
		sum = a[i__ + j * a_dim1];
		if (i__ > 1) {
		    i__3 = i__ - 1;
		    for (k = 1; k <= i__3; ++k) {
			sum -= a[i__ + k * a_dim1] * a[k + j * a_dim1];
		    }
		    a[i__ + j * a_dim1] = sum;
		}
	    }
	}
	aamax = 0.f;
	i__2 = *n;
	for (i__ = j; i__ <= i__2; ++i__) {
	    sum = a[i__ + j * a_dim1];
	    if (j > 1) {
		i__3 = j - 1;
		for (k = 1; k <= i__3; ++k) {
		    sum -= a[i__ + k * a_dim1] * a[k + j * a_dim1];
		}
		a[i__ + j * a_dim1] = sum;
	    }
	    dum = vv[i__ - 1] * dabs(sum);
	    if (dum >= aamax) {
		imax = i__;
		aamax = dum;
	    }
	}
	if (j != imax) {
	    i__2 = *n;
	    for (k = 1; k <= i__2; ++k) {
		dum = a[imax + k * a_dim1];
		a[imax + k * a_dim1] = a[j + k * a_dim1];
		a[j + k * a_dim1] = dum;
	    }
	    *d__ = -(*d__);
	    vv[imax - 1] = vv[j - 1];
	}
	index[j] = imax;
	L__1 = a[j + j * a_dim1] != 0.f;
	assert_(&L__1, "Singular matrix", (ftnlen)15);
	if (j != *n) {
	    dum = 1.f / a[j + j * a_dim1];
	    i__2 = *n;
	    for (i__ = j + 1; i__ <= i__2; ++i__) {
		a[i__ + j * a_dim1] *= dum;
	    }
	}
    }
    return 0;
} /* ludcmp_ */

/* CC	corresponding back-substitution routine */
/* C */

/* Subroutine */ int lubksb_(real *a, integer *n, integer *index, real *b)
{
    /* System generated locals */
    integer a_dim1, a_offset, i__1, i__2;

    /* Local variables */
    static integer i__, j, ii, ll;
    static real sum;



    /* Parameter adjustments */
    --b;
    --index;
    a_dim1 = *n;
    a_offset = 1 + a_dim1;
    a -= a_offset;

    /* Function Body */
    ii = 0;
    i__1 = *n;
    for (i__ = 1; i__ <= i__1; ++i__) {
	ll = index[i__];
	sum = b[ll];
	b[ll] = b[i__];
	if (ii != 0) {
	    i__2 = i__ - 1;
	    for (j = ii; j <= i__2; ++j) {
		sum -= a[i__ + j * a_dim1] * b[j];
	    }
	} else if (sum != 0.f) {
	    ii = i__;
	}
	b[i__] = sum;
    }

    for (i__ = *n; i__ >= 1; --i__) {
	sum = b[i__];
	if (i__ < *n) {
	    i__1 = *n;
	    for (j = i__ + 1; j <= i__1; ++j) {
		sum -= a[i__ + j * a_dim1] * b[j];
	    }
	}
	b[i__] = sum / a[i__ + i__ * a_dim1];
    }
    return 0;
} /* lubksb_ */

/* *********************************************************************** */
/* Find eigenvalues and eigenvectors of nxn symmetric matrix using      * */
/* cyclic Jacobi method. NP is (Fortran) physical storage dimension.    * */
/* On return A is destroyed, D contains eigenvalues,and each column of  * */
/* V is a normalized eigenvector                                        * */
/* Adapted from Numerical Recipes in Fortran (1986)                     * */
/* We only need it for 3x3 symmetric matrices, so it's overkill.        * */
/* *********************************************************************** */
/* CC */
/* C */

/* Subroutine */ int jacobi_(real *a, integer *n, integer *np, real *d__, 
	real *v)
{
    /* System generated locals */
    integer a_dim1, a_offset, v_dim1, v_offset, i__1, i__2, i__3;
    real r__1, r__2, r__3, r__4;
    logical L__1;

    /* Builtin functions */
    double sqrt(doublereal);

    /* Local variables */
    static real b[4], c__, g, h__;
    static integer i__, j;
    static real s, t, z__[4];
    static integer ip, iq;
    static real sm, tau;
    static integer nrot;
    static real theta, thresh;
    extern /* Subroutine */ int assert_(logical *, char *, ftnlen);

/* 	Machine dependent! (converge when off-diagonal sum is less than this) */

    /* Parameter adjustments */
    --d__;
    v_dim1 = *np;
    v_offset = 1 + v_dim1;
    v -= v_offset;
    a_dim1 = *np;
    a_offset = 1 + a_dim1;
    a -= a_offset;

    /* Function Body */
    L__1 = *n <= 4;
    assert_(&L__1, "Matrix too big for eigenvector routine", (ftnlen)38);

/* 	Initialize V to identity matrix, B and D to diagonal of A */
    i__1 = *n;
    for (ip = 1; ip <= i__1; ++ip) {
	i__2 = *n;
	for (iq = 1; iq <= i__2; ++iq) {
	    v[ip + iq * v_dim1] = 0.f;
	}
	v[ip + ip * v_dim1] = 1.f;
	b[ip - 1] = a[ip + ip * a_dim1];
	d__[ip] = b[ip - 1];
	z__[ip - 1] = 0.f;
    }
    nrot = 0;

/* 	50 sweeps is never expected to happen; Press et al (1986) claim */
/* 	6 to 10 are typical for moderate matrix sizes */
/* 	Empirical trials of rastep show 6 sweeps, 8-10 rotations */
    for (i__ = 1; i__ <= 50; ++i__) {
	sm = 0.f;
	i__1 = *n - 1;
	for (ip = 1; ip <= i__1; ++ip) {
	    i__2 = *n;
	    for (iq = ip + 1; iq <= i__2; ++iq) {
		sm += (r__1 = a[ip + iq * a_dim1], dabs(r__1));
	    }
	}
	if (sm < 1e-37f) {
	    return 0;
	}
/* 	    After 4 sweeps skip rotation if the off-diagonal is small */
	if (i__ < 4) {
	    thresh = sm * .2f / (*n * *n);
	} else {
	    thresh = 0.f;
	}
	i__1 = *n - 1;
	for (ip = 1; ip <= i__1; ++ip) {
	    i__2 = *n;
	    for (iq = ip + 1; iq <= i__2; ++iq) {
		g = (r__1 = a[ip + iq * a_dim1], dabs(r__1)) * 100.f;
		if (i__ > 4 && (r__1 = d__[ip], dabs(r__1)) + g == (r__2 = 
			d__[ip], dabs(r__2)) && (r__3 = d__[iq], dabs(r__3)) 
			+ g == (r__4 = d__[iq], dabs(r__4))) {
		    a[ip + iq * a_dim1] = 0.f;
		} else if ((r__1 = a[ip + iq * a_dim1], dabs(r__1)) > thresh) 
			{
		    h__ = d__[iq] - d__[ip];
		    if (dabs(h__) + g == dabs(h__)) {
			t = a[ip + iq * a_dim1] / h__;
		    } else {
			theta = h__ * .5f / a[ip + iq * a_dim1];
			t = 1.f / (dabs(theta) + sqrt(theta * theta + 1.f));
			if (theta < 0.f) {
			    t = -t;
			}
		    }
		    c__ = 1.f / sqrt(t * t + 1.f);
		    s = t * c__;
		    tau = s / (c__ + 1.f);
		    h__ = t * a[ip + iq * a_dim1];
		    z__[ip - 1] -= h__;
		    z__[iq - 1] += h__;
		    d__[ip] -= h__;
		    d__[iq] += h__;
		    a[ip + iq * a_dim1] = 0.f;
		    i__3 = ip - 1;
		    for (j = 1; j <= i__3; ++j) {
			g = a[j + ip * a_dim1];
			h__ = a[j + iq * a_dim1];
			a[j + ip * a_dim1] = g - s * (h__ + g * tau);
			a[j + iq * a_dim1] = h__ + s * (g - h__ * tau);
		    }
		    i__3 = iq - 1;
		    for (j = ip + 1; j <= i__3; ++j) {
			g = a[ip + j * a_dim1];
			h__ = a[j + iq * a_dim1];
			a[ip + j * a_dim1] = g - s * (h__ + g * tau);
			a[j + iq * a_dim1] = h__ + s * (g - h__ * tau);
		    }
		    i__3 = *n;
		    for (j = iq + 1; j <= i__3; ++j) {
			g = a[ip + j * a_dim1];
			h__ = a[iq + j * a_dim1];
			a[ip + j * a_dim1] = g - s * (h__ + g * tau);
			a[iq + j * a_dim1] = h__ + s * (g - h__ * tau);
		    }
		    i__3 = *n;
		    for (j = 1; j <= i__3; ++j) {
			g = v[j + ip * v_dim1];
			h__ = v[j + iq * v_dim1];
			v[j + ip * v_dim1] = g - s * (h__ + g * tau);
			v[j + iq * v_dim1] = h__ + s * (g - h__ * tau);
		    }
		    ++nrot;
		}
	    }
	}

/* 	    Update D with sum and reinitialize Z */
	i__1 = *n;
	for (ip = 1; ip <= i__1; ++ip) {
	    b[ip - 1] += z__[ip - 1];
	    d__[ip] = b[ip - 1];
	    z__[ip - 1] = 0.f;
	}
    }

    assert_(&c_false, "Failed to find eigenvectors", (ftnlen)27);
    return 0;
} /* jacobi_ */


/* 	This one really has nothing to do with quadrics per se, */
/* 	but since it's invoked by qinp, both render and rastep */
/* 	need to be able to see it. */
/* 	Should really be in separate file of support routines */

doublereal persp_(real *z__)
{
    /* System generated locals */
    real ret_val;

    if (*z__ / matrices_1.eyepos > .999f) {
	ret_val = 1e3f;
    } else {
	ret_val = 1.f / (1.f - *z__ / matrices_1.eyepos);
    }
    return ret_val;
} /* persp_ */