/*****************************************************************************
* Adaptive Iso Shader program to render freeform models.		     *
*									     *
* Written by:  Gershon Elber			     Ver 1.0, January 1994   *
*****************************************************************************/

#include <stdio.h>
#include <math.h>
#include <string.h>
#include "irit_sm.h"
#include "config.h"
#include "iritgrap.h"
#include "attribut.h"
#include "allocate.h"
#include "iritgrap.h"
#include "poly_cln.h"
#include "geomat3d.h"
#include "getarg.h"
#include "ip_cnvrt.h"

typedef enum {		           /* Type of distance function computation. */
    AISO_DIST_REGULAR = 0,
    AISO_DIST_PHONG,
    AISO_DIST_PHONG_2L,
    AISO_DIST_PHONG_SPEC,
    AISO_DIST_PHONG_2L_SPEC,
    AISO_DIST_ZNORMAL,
    AISO_DIST_POINT_E3
} AIsoDistCompType;

#define NORMAL_ANGLE_SCALE	1e5
#define TRUNCATED_DOMAIN	1e-3
#define MIN_DOMAIN_VALID	1e-2
#define MAX_LEVEL_VALID		100

#ifdef NO_CONCAT_STR
static char *VersionStr =
    "AIsoShad		Version 4.0,	Gershon Elber,\n\
	 (C) Copyright 1989/90/91/92/93 Gershon Elber, Non commercial use only.";
#else
static char *VersionStr =
    "AIsoShad		" VERSION ",	Gershon Elber,	"
	__DATE__ ",   " __TIME__ "\n" COPYRIGHT ", Non commercial use only.";
#endif /* NO_CONCAT_STR */

static char *CtrlStr =
#ifdef DOUBLE
    "AIsoShad o%-OutName!s m%- i%- r%-RndrMdl!d c%-CosPwr!F s%-SdrPwr!F l%-Lx|Ly|Lz!F!F!F R%-Random!F d%-AdapDir!d t%-SrfZTrans!F M%-MinSubdiv!d D%-AdapDist!F w%-AdapIsoWidth!F W%- z%- DFiles!*s";
#else
    "AIsoShad o%-OutName!s m%- i%- r%-RndrMdl!d c%-CosPwr!f s%-SdrPwr!f l%-Lx|Ly|Lz!f!f!f R%-Random!F d%-AdapDir!d t%-SrfZTrans!f M%-MinSubdiv!d D%-AdapDist!f w%-AdapIsoWidth!f W%- z%- DFiles!*s";
#endif

static char
    *GlblLightSourceStr = "",
    *GlblOutFileName = "aisoshad.dat";

static AIsoDistCompType
    GlblDistRndrModel = AISO_DIST_PHONG;

static int
    GlblSymbInterp = TRUE,
    GlblAdapIsoColor = IG_IRIT_WHITE,
    GlblAdapIsoDir[2] = { CAGD_CONST_U_DIR, CAGD_NO_DIR },
    GlblTalkative = FALSE,
    GlblMinSubdivLevel = 1,
    GlblVariableWidth = FALSE;

static RealType
    GlblCosinePower = 1,
    GlblShaderPower = 2,
    GlblLightSource[3] = { 1.0, 0.0, 0.0 },
    GlblSrfTranslate[3] = { 0.0, 0.0, 0.01 },
    GlblRandomDist = 1.0,
    GlblAdapDistance = 1.0,
    GlblAdapIsoWidth = 0.001,
    GlblVarWidthSegLen = 0.1;

static ConfigStruct SetUp[] =
{
  { "AdapDir",        "-d", (VoidPtr) &GlblAdapIsoDir[0],   SU_INTEGER_TYPE },
  { "AdapIsoColor",   "",   (VoidPtr) &GlblAdapIsoColor,    SU_INTEGER_TYPE },
  { "RenderModel",    "-r", (VoidPtr) &GlblDistRndrModel,   SU_INTEGER_TYPE },
  { "MinSubdiv",      "-M", (VoidPtr) &GlblMinSubdivLevel,  SU_INTEGER_TYPE },
  { "SymbInterp",     "-i", (VoidPtr) &GlblSymbInterp,	    SU_BOOLEAN_TYPE },
  { "MoreVerbose",    "-m", (VoidPtr) &GlblTalkative,	    SU_BOOLEAN_TYPE },
  { "LightSrcDir",    "-l", (VoidPtr) &GlblLightSourceStr,  SU_STRING_TYPE },
  { "AdapIsoWidth",   "-w", (VoidPtr) &GlblAdapIsoWidth,    SU_REAL_TYPE },
  { "CosinePower",    "-c", (VoidPtr) &GlblCosinePower,     SU_REAL_TYPE },
  { "ShaderPower",    "-s", (VoidPtr) &GlblShaderPower,     SU_REAL_TYPE },
  { "RandomDist",     "-R", (VoidPtr) &GlblRandomDist,      SU_REAL_TYPE },
  { "SrfZTrans",      "-t", (VoidPtr) &GlblSrfTranslate[2], SU_REAL_TYPE },
  { "AdapDistance",   "-D", (VoidPtr) &GlblAdapDistance,    SU_REAL_TYPE },
};
#define NUM_SET_UP	(sizeof(SetUp) / sizeof(ConfigStruct))

static CagdCrvStruct *AdapIsoDistSqr(int Level,
				     CagdCrvStruct *Crv1,
				     CagdCrvStruct *NCrv1,
				     CagdCrvStruct *Crv2,
				     CagdCrvStruct *NCrv2);
static RealType IntPower2(RealType t, int Power);
static void AIsoShaderExit(int ExitCode);

/*****************************************************************************
* DESCRIPTION:                                                               M
* Main module of AIsoShad - Read command line and do what is needed...	     M
*                                                                            *
* PARAMETERS:                                                                M
*   argc, argv:  Command line.                                               M
*                                                                            *
* RETURN VALUE:                                                              M
*   void                                                                     M
*                                                                            *
* KEYWORDS:                                                                  M
*   main                                                                     M
*****************************************************************************/
void main(int argc, char **argv)
{
    int Error,
	DistCompFlag = FALSE,
	CosinePowerFlag = FALSE,
	ShaderPowerFlag = FALSE,
	AdapDirFlag = FALSE,
	SrfZTransFlag = FALSE,
	AdapDistanceFlag = FALSE,
	AdapIsoWidthFlag = FALSE,
        LightSrcFlag = FALSE,
        MinSubdivFlag = FALSE,
	RandomFlag = FALSE,
	OutFileFlag = FALSE,
	VerFlag = FALSE,
	NumFiles = 0;
    char
	**FileNames = NULL;
    FILE *f;
    IPObjectStruct *PObjects, *NoProcessObjs, *PObj,
	*PObjWaves = NULL;
    MatrixType CrntViewMat;

    Config("aisoshad", SetUp, NUM_SET_UP);   /* Read config. file if exists. */
    if (GlblLightSourceStr != NULL && strlen(GlblLightSourceStr) > 0) {
#ifdef DOUBLE
	if (sscanf(GlblLightSourceStr, "%lf,%lf,%lf",
#else
	if (sscanf(GlblLightSourceStr, "%f,%f,%f",
#endif /* DOUBLE */
		   &GlblLightSource[0],
		   &GlblLightSource[1],
		   &GlblLightSource[2]) != 3) {
	    fprintf(stderr,
		    "Fail to parse LightSrcDir in configuration file.\n");
	    AIsoShaderExit(-1);
	}
    }

    if ((Error = GAGetArgs(argc, argv, CtrlStr,
			   &OutFileFlag, &GlblOutFileName, &GlblTalkative,
			   &GlblSymbInterp, &DistCompFlag, &GlblDistRndrModel,
			   &CosinePowerFlag, &GlblCosinePower,
			   &ShaderPowerFlag, &GlblShaderPower,
			   &LightSrcFlag, &GlblLightSource[0],
			   &GlblLightSource[1], &GlblLightSource[2],
			   &RandomFlag, &GlblRandomDist,
			   &AdapDirFlag, &GlblAdapIsoDir[0],
			   &SrfZTransFlag, &GlblSrfTranslate[2],
			   &MinSubdivFlag, &GlblMinSubdivLevel,
			   &AdapDistanceFlag, &GlblAdapDistance,
			   &AdapIsoWidthFlag, &GlblAdapIsoWidth,
			   &GlblVariableWidth,
			   &VerFlag, &NumFiles, &FileNames)) != 0) {
	GAPrintErrMsg(Error);
	GAPrintHowTo(CtrlStr);
	AIsoShaderExit(1);
    }

    if (VerFlag) {
	fprintf(stderr, "\n%s\n\n", VersionStr);
	GAPrintHowTo(CtrlStr);
	ConfigPrint(SetUp, NUM_SET_UP);
	AIsoShaderExit(0);
    }

    if (!NumFiles) {
	fprintf(stderr, "No data file names where given, exit.\n");
	GAPrintHowTo(CtrlStr);
	AIsoShaderExit(1);
    }

    switch (GlblAdapIsoDir[0]) {
        case 0:
	    GlblAdapIsoDir[0] = CAGD_CONST_U_DIR;
	    break;
	case 1:
	    GlblAdapIsoDir[0] = CAGD_CONST_V_DIR;
	    break;
	case 2:
	    GlblAdapIsoDir[0] = CAGD_CONST_U_DIR;
	    GlblAdapIsoDir[1] = CAGD_CONST_V_DIR;
	    break;
	default:
	    break;
    }

    BspMultInterpFlag(GlblSymbInterp);
    SymbSetAdapIsoExtractMinLevel(GlblMinSubdivLevel);

    if (GlblDistRndrModel != AISO_DIST_POINT_E3)
	PT_NORMALIZE(GlblLightSource);

    if (GlblRandomDist < 1) {
	fprintf(stderr,	"Randomization [-R] should not be less than one\n");
	AIsoShaderExit(0);
    }

    /* Get the data files: */
    if ((PObjects = IritPrsrGetDataFiles(FileNames, NumFiles, TRUE, FALSE)) ==
									NULL)
	AIsoShaderExit(0);

    if (IritPrsrWasPrspMat)
	MatMultTwo4by4(CrntViewMat, IritPrsrViewMat, IritPrsrPrspMat);
    else
	GEN_COPY(CrntViewMat, IritPrsrViewMat, sizeof(MatrixType));

    PObj = GMTransformObjectList(PObjects, CrntViewMat);
    IPFreeObjectList(PObjects);
    PObjects = PObj;

    /* Open output file, if necessary. */
    if (OutFileFlag) {
	if ((f = fopen(GlblOutFileName, "w")) == NULL) {
	    fprintf(stderr, "Failed to open \"%s\".\n", GlblOutFileName);
	    AIsoShaderExit(2);
	}
    }
    else
	f = stdout;

    /* Traverse all the objects and generate the coverage for them. */
    for (PObj = PObjects; PObj != NULL; PObj = PObj -> Pnext) {
	if (GlblTalkative)
	    fprintf(stderr, "Processing object \"%s\"\n", PObj -> Name);

	if (IP_IS_SRF_OBJ(PObj)) {
	    CagdSrfStruct *Srf;
	    int SrfAdapIsoDir[2],
		SrfAdapIsoMinSubdiv = AttrGetObjectIntAttrib(PObj,
							  "AdapIsoMinSubdiv");
	    RealType
		RelativeAdapIsoDist = AttrGetObjectRealAttrib(PObj,
							      "AdapIsoDist");

	    if (RelativeAdapIsoDist > IP_ATTR_BAD_REAL / 10.0)
		RelativeAdapIsoDist = 1.0;

	    SrfAdapIsoDir[0] = AttrGetObjectIntAttrib(PObj, "AdapIsoDir");
	    switch (SrfAdapIsoDir[0]) {
		case 0:
		    SrfAdapIsoDir[0] = CAGD_CONST_U_DIR;
		    break;
		case 1:
		    SrfAdapIsoDir[0] = CAGD_CONST_V_DIR;
		    break;
		case 2:
		    SrfAdapIsoDir[0] = CAGD_CONST_U_DIR;
		    SrfAdapIsoDir[1] = CAGD_CONST_V_DIR;
		    break;
		default:
		    SrfAdapIsoDir[0] = GlblAdapIsoDir[0];
		    SrfAdapIsoDir[1] = GlblAdapIsoDir[1];
		    break;
	    }

	    if (SrfAdapIsoMinSubdiv != IP_ATTR_BAD_INT)
	        SymbSetAdapIsoExtractMinLevel(SrfAdapIsoMinSubdiv);

	    /* Set the color of the surfaces to black. */
	    AttrSetObjectColor(PObj, IG_IRIT_BLACK);
	    AttrSetObjectRGBColor(PObj, 0, 0, 0);

	    IritPrsrPutObjectToFile(f, PObj);  /* Dump it as is. */

	    for (Srf = PObj -> U.Srfs; Srf != NULL; Srf = Srf -> Pnext) {
		RealType UMin, UMax, VMin, VMax;
		IPObjectStruct *PAdapIsoObj;
		CagdSrfStruct *TSrf, *TSrf1,
		    *NSrf = SymbSrfNormalSrf(Srf);
		CagdCrvStruct *Crv, *AdapIso;

		/* Remove a little from all four boundaries to prevent from */
		/* dealing with degenerated cases so common on boundary.    */
		CagdSrfDomain(Srf, &UMin, &UMax, &VMin, &VMax);
		TSrf1 = CagdSrfRegionFromSrf(Srf,
					     UMin + TRUNCATED_DOMAIN,
					     UMax - TRUNCATED_DOMAIN,
					     CAGD_CONST_U_DIR);
		TSrf = CagdSrfRegionFromSrf(TSrf1,
					    VMin + TRUNCATED_DOMAIN,
					    VMax - TRUNCATED_DOMAIN,
					    CAGD_CONST_V_DIR);
					     
		CagdSrfFree(TSrf1);

		CagdSrfTransform(TSrf, GlblSrfTranslate, 1.0);
		AdapIso = SymbAdapIsoExtract(TSrf, NSrf, AdapIsoDistSqr,
					     SrfAdapIsoDir[0],
					     GlblAdapDistance *
					         RelativeAdapIsoDist,
					     FALSE, FALSE);
		if (SrfAdapIsoDir[1] != CAGD_NO_DIR) {
		    CagdCrvStruct *AdapIsoLast,
		        *AdapIso2 = SymbAdapIsoExtract(TSrf, NSrf,
						       AdapIsoDistSqr,
						       SrfAdapIsoDir[1],
						       GlblAdapDistance *
						          RelativeAdapIsoDist,
						       FALSE, FALSE);

		    if (AdapIso != NULL) {
		        AdapIsoLast = CagdListLast(AdapIso);

			AdapIsoLast -> Pnext = AdapIso2;
		    }
		    else
		        AdapIso = AdapIso2;
		}

		CagdSrfFree(NSrf);
		CagdSrfFree(TSrf);

		/* Remove the normal curves from the list. */
		for (Crv = AdapIso; Crv != NULL; Crv = Crv -> Pnext) {
		    CagdCrvStruct
			*NextCrv = Crv -> Pnext -> Pnext;

		    Crv -> Pnext -> Pnext = NULL;
		    CagdCrvFree(Crv -> Pnext);
		    Crv -> Pnext = NextCrv;
		}

		if (GlblVariableWidth) {
		    /* Save every curve with its variable width function. */
		    for (Crv = AdapIso; Crv != NULL; ) {
			CagdCrvStruct
			    *NextCrv = Crv -> Pnext;

			Crv -> Pnext = NULL;

			PAdapIsoObj = GenCRVObject(Crv);
			PAdapIsoObj -> Attrs = Crv -> Attr;
			Crv -> Attr = NULL;
			AttrSetObjectColor(PAdapIsoObj, GlblAdapIsoColor);

			IritPrsrPutObjectToFile(f, PAdapIsoObj);
			IPFreeObject(PAdapIsoObj);

			Crv = NextCrv;
		    }
		}
		else {
		    PAdapIsoObj = GenCRVObject(AdapIso);
		    AttrSetObjectColor(PAdapIsoObj, GlblAdapIsoColor);
		    AttrSetObjectRealAttrib(PAdapIsoObj, "width",
					    GlblAdapIsoWidth);

		    IritPrsrPutObjectToFile(f, PAdapIsoObj);
		    IPFreeObject(PAdapIsoObj);
		}
	    }

	    SymbSetAdapIsoExtractMinLevel(GlblMinSubdivLevel);
	}
    }

    AIsoShaderExit(0);
}

/*****************************************************************************
* DESCRIPTION:                                                               *
*   Computes a distance square function between Crv1 and Crv2 as a scalar    *
* field. The normal vector fields (not unit size) of the two curves are also *
* provided as NCrv1, NCrv2.						     *
*                                                                            *
* PARAMETERS:                                                                *
*   Level:     Maximum depth of recursion of adaptive iso algorithm.         *
*   Crv1:      First curve to compute distance field from it to Crv2.        *
*   NCrv1:     Normal vector field of Crv1.                                  *
*   Crv2:      Second curve to compute distance field from it to Crv1.       *
*   NCrv2:     Normal vector field of Crv2.                                  *
*                                                                            *
* RETURN VALUE:                                                              *
*   CagdCrvStruct *:  A scalar distance function field.                      *
*****************************************************************************/
static CagdCrvStruct *AdapIsoDistSqr(int Level,
				     CagdCrvStruct *Crv1,
				     CagdCrvStruct *NCrv1,
				     CagdCrvStruct *Crv2,
				     CagdCrvStruct *NCrv2)
{
    int i;
    CagdRType TMin, TMax, *Points1, *Points2;
    CagdCrvStruct
	*DistCrv2Aux1 = NULL,
	*DistCrv2Aux2 = NULL,
        *DiffCrv = SymbCrvSub(Crv1, Crv2),
        *DiffCrv2D = CagdCoerceCrvTo(DiffCrv, CAGD_PT_E2_TYPE),
        *DistCrv2 = SymbCrvDotProd(DiffCrv2D, DiffCrv);

    CagdCrvFree(DiffCrv);
    CagdCrvFree(DiffCrv2D);

    CagdCrvDomain(Crv1, &TMin, &TMax);

    if (GlblVariableWidth) {
	if (DistCrv2 -> PType != CAGD_PT_E1_TYPE)
	    DistCrv2Aux1 = CagdCoerceCrvTo(DistCrv2, CAGD_PT_E1_TYPE);
	else
	    DistCrv2Aux1 = CagdCrvCopy(DistCrv2);
	DistCrv2Aux2 = CagdCrvCopy(DistCrv2Aux1);

	Points1 = DistCrv2Aux1 -> Points[1];
	Points2 = DistCrv2Aux2 -> Points[1];
    }

    if (TMax - TMin > MIN_DOMAIN_VALID && Level < MAX_LEVEL_VALID) {
	CagdRType TMin1, TMax1, TMin2, TMax2,
	    *Points = DistCrv2 -> Points[1],
	    *Nodes = CagdCrvNodes(DistCrv2);

	CagdCrvDomain(NCrv1, &TMin1, &TMax1);
	CagdCrvDomain(NCrv2, &TMin2, &TMax2);

	switch (GlblDistRndrModel) {
	    case AISO_DIST_PHONG:
	    case AISO_DIST_PHONG_2L:
	    case AISO_DIST_PHONG_SPEC:
	    case AISO_DIST_PHONG_2L_SPEC:
	        for (i = 0; i < DistCrv2 -> Length; i++) {
		    CagdVType N1, N2, RefDir;
		    CagdRType *R, Z, Z1, Z2,
		        t1 = BOUND(Nodes[i], TMin1, TMax1),
		        t2 = BOUND(Nodes[i], TMin2, TMax2);

		    R = CagdCrvEval(NCrv1, t1);
		    CagdCoerceToE3(N1, &R, -1, NCrv1 -> PType);
		    R = CagdCrvEval(NCrv2, t2);
		    CagdCoerceToE3(N2, &R, -1, NCrv2 -> PType);

		    PT_NORMALIZE(N1);
		    PT_NORMALIZE(N2);

		    if (GlblVariableWidth) {
			if (GlblDistRndrModel == AISO_DIST_PHONG_2L ||
			    GlblDistRndrModel == AISO_DIST_PHONG_2L_SPEC) {
			    Z1 = fabs(DOT_PROD(N1, GlblLightSource));
			    Z2 = fabs(DOT_PROD(N2, GlblLightSource));

			    if (GlblDistRndrModel == AISO_DIST_PHONG_2L_SPEC) {
				/* Add specular term, if have one. */
				if (DOT_PROD(N1, GlblLightSource) > 0.0) {
				    PT_COPY(RefDir, N1);
				    PT_SCALE(RefDir, 2 * Z);
				    PT_SUB(RefDir, RefDir, GlblLightSource);
				    Z1 += pow(fabs(RefDir[2]),
					      GlblCosinePower);
				}
				else {
				    PT_COPY(RefDir, N1);
				    PT_SCALE(RefDir, -2 * Z);
				    PT_SUB(RefDir, RefDir, GlblLightSource);
				    Z1 += pow(fabs(RefDir[2]),
					      GlblCosinePower);
				}

				if (DOT_PROD(N2, GlblLightSource) > 0.0) {
				    PT_COPY(RefDir, N2);
				    PT_SCALE(RefDir, 2 * Z);
				    PT_SUB(RefDir, RefDir, GlblLightSource);
				    Z2 += pow(fabs(RefDir[2]),
					      GlblCosinePower);
				}
				else {
				    PT_COPY(RefDir, N2);
				    PT_SCALE(RefDir, -2 * Z);
				    PT_SUB(RefDir, RefDir, GlblLightSource);
				    Z2 += pow(fabs(RefDir[2]),
					      GlblCosinePower);
				}
			    }
			}
			else { /* AISO_DIST_PHONG || AISO_DIST_PHONG_SPEC */
			    Z1 = (DOT_PROD(N1, GlblLightSource));
			    Z1 = (Z1 + 1.0) / 2.0;
			    Z2 = (DOT_PROD(N2, GlblLightSource));
			    Z2 = (Z2 + 1.0) / 2.0;

			    if (GlblDistRndrModel == AISO_DIST_PHONG_SPEC) {
				/* Add specular term, if have one. */
				if (DOT_PROD(N1, GlblLightSource) > 0.0) {
				    PT_COPY(RefDir, N1);
				    PT_SCALE(RefDir, 2 * Z);
				    PT_SUB(RefDir, RefDir, GlblLightSource);
				    Z1 += pow(fabs(RefDir[2]),
					      GlblCosinePower);
				}

				if (DOT_PROD(N2, GlblLightSource) > 0.0) {
				    PT_COPY(RefDir, N2);
				    PT_SCALE(RefDir, 2 * Z);
				    PT_SUB(RefDir, RefDir, GlblLightSource);
				    Z2 += pow(fabs(RefDir[2]),
					      GlblCosinePower);
				}
			    }
			}

			Z1 = pow(Z1, GlblShaderPower);
			Z2 = pow(Z2, GlblShaderPower);
	
			Points1[i] = Z1 + 1.0 / NORMAL_ANGLE_SCALE;
			Points2[i] = Z2 + 1.0 / NORMAL_ANGLE_SCALE;
		    }
		    else {
			if (GlblDistRndrModel == AISO_DIST_PHONG_2L ||
			    GlblDistRndrModel == AISO_DIST_PHONG_2L_SPEC) {
			    if (DOT_PROD(N1, N2) > 0.0) {
				Z = (DOT_PROD(N1, GlblLightSource) +
				     DOT_PROD(N2, GlblLightSource)) / 2.0;
				Z = fabs(Z);

				if (GlblDistRndrModel == AISO_DIST_PHONG_2L_SPEC) {
				    PT_BLEND(RefDir, N1, N2, 0.5);
				    PT_NORMALIZE(RefDir);
			
				    /* Add specular term, if have one. */
				    if (DOT_PROD(RefDir,
						 GlblLightSource) > 0.0) {
					PT_SCALE(RefDir, 2 * Z);
					PT_SUB(RefDir, RefDir,
					       GlblLightSource);
					Z += pow(fabs(RefDir[2]),
						 GlblCosinePower);
				    }
				    else {
					PT_SCALE(RefDir, -2 * Z);
					PT_SUB(RefDir, RefDir,
					       GlblLightSource);
					Z += pow(fabs(RefDir[2]),
						 GlblCosinePower);
				    }
				}
			    }
			    else
			      Z = 1.0;
			}
			else { /* AISO_DIST_PHONG || AISO_DIST_PHONG_SPEC */
			    if (DOT_PROD(N1, N2) > 0.0) {
				Z = (DOT_PROD(N1, GlblLightSource) +
				     DOT_PROD(N2, GlblLightSource)) / 2.0;
				Z = (Z + 1.0) / 2.0;

				if (GlblDistRndrModel == AISO_DIST_PHONG_SPEC) {
				    PT_BLEND(RefDir, N1, N2, 0.5);
				    PT_NORMALIZE(RefDir);
			
				    /* Add specular term, if have one. */
				    if (DOT_PROD(RefDir,
						 GlblLightSource) > 0.0) {
					PT_SCALE(RefDir, 2 * Z);
					PT_SUB(RefDir, RefDir,
					       GlblLightSource);
					Z += pow(fabs(RefDir[2]),
						 GlblCosinePower);
				    }
				}
			    }
			    else
			      Z = 1.0;
			}

			Z = pow(Z, GlblShaderPower);

			Points[i] *= Z + 1.0 / NORMAL_ANGLE_SCALE;
		    }
		}
	        break;
	    case AISO_DIST_ZNORMAL:
		for (i = 0; i < DistCrv2 -> Length; i++) {
		    CagdVType N1, N2, RefDir;
		    CagdRType *R, Z, Z1, Z2,
		        t1 = BOUND(Nodes[i], TMin1, TMax1),
		        t2 = BOUND(Nodes[i], TMin2, TMax2);

		    R = CagdCrvEval(NCrv1, t1);
		    CagdCoerceToE3(N1, &R, -1, NCrv1 -> PType);
		    R = CagdCrvEval(NCrv2, t2);
		    CagdCoerceToE3(N2, &R, -1, NCrv2 -> PType);

		    PT_NORMALIZE(N1);
		    PT_NORMALIZE(N2);

		    if (GlblVariableWidth) {
			Z1 = 1.0 + GlblCosinePower - fabs(N1[2]),
			Z2 = 1.0 + GlblCosinePower - fabs(N2[2]);

			Z1 = pow(Z1, GlblShaderPower);
			Z2 = pow(Z2, GlblShaderPower);

			Points1[i] = Z1;
			Points2[i] = Z2;
		    }
		    else {
			Z = 1.0 + GlblCosinePower -
			                   (fabs(N1[2]) + fabs(N2[2])) / 2.0;

			Z = pow(Z, GlblShaderPower);

			Points[i] *= Z;
		    }
		}
	        break;
	    case AISO_DIST_POINT_E3:
		for (i = 0; i < DistCrv2 -> Length; i++) {
		    CagdVType P1, P2;
		    CagdRType *R, Z, Z1, Z2,
		        t1 = BOUND(Nodes[i], TMin1, TMax1),
		        t2 = BOUND(Nodes[i], TMin2, TMax2);

		    R = CagdCrvEval(Crv1, t1);
		    CagdCoerceToE3(P1, &R, -1, NCrv1 -> PType);
		    R = CagdCrvEval(Crv2, t2);
		    CagdCoerceToE3(P2, &R, -1, NCrv2 -> PType);

		    if (GlblVariableWidth) {
			Z1 = CGDistPointPoint(P1, GlblLightSource);
			Z2 = CGDistPointPoint(P2, GlblLightSource);

			Z1 = pow(Z1, GlblShaderPower);
			Z2 = pow(Z2, GlblShaderPower);

			Points1[i] = 1 / (Z1 + GlblCosinePower);
			Points2[i] = 1 / (Z2 + GlblCosinePower);
		    }
		    else {
			PT_BLEND(P1, P1, P2, 0.5);
			Z = CGDistPointPoint(P1, GlblLightSource);

			Z = pow(Z, GlblShaderPower);

			Points[i] *= 1 / (Z + GlblCosinePower);
		    }
		}
	        break;
	    default:
		break;
	}

	IritFree((VoidPtr) Nodes);

	if (GlblRandomDist > 0.0) {
	    CagdRType
		*Points = DistCrv2 -> Points[1];

	    for (i = 0; i < DistCrv2 -> Length; i++)
		Points[i] *= IritRandom(1 / GlblRandomDist,
					GlblRandomDist);

	    if (GlblVariableWidth) {
		CagdRType
		    *Points1 = DistCrv2Aux1 -> Points[1],
		    *Points2 = DistCrv2Aux2 -> Points[1];

		for (i = 0; i < DistCrv2 -> Length; i++) {
		    Points1[i] *= IritRandom(1 / GlblRandomDist,
					     GlblRandomDist);
		    Points2[i] *= IritRandom(1 / GlblRandomDist,
					     GlblRandomDist);
		}
	    }
	}
    }

    if (GlblVariableWidth) {
	AttrSetObjAttrib(&Crv1 -> Attr, "VarWidth",
			 GenCRVObject(DistCrv2Aux1), FALSE);
	AttrSetObjAttrib(&Crv2 -> Attr, "VarWidth",
			 GenCRVObject(DistCrv2Aux2), FALSE);
    }

    return DistCrv2;
}

/*****************************************************************************
* DESCRIPTION:                                                               *
*   Approximates t to the power of Power, as powers of twos only.            *
*                                                                            *
* PARAMETERS:                                                                *
*   t:         To raise to the power of Power.                               *
*   Power:     If Power is of the form 2^n, n integer, result is exact.      *
*                                                                            *
* RETURN VALUE:                                                              *
*   RealType:  t^Power                                                       *
*****************************************************************************/
static RealType IntPower2(RealType t, int Power)
{
    int i;

    if (Power == 0)
        return 1.0;

    for (i = 1; i < Power; i <<= 1)
	t *= t;

    return t;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* AIsoShad Exit routine.						     M
*                                                                            *
* PARAMETERS:                                                                M
*   ExitCode:    To notify O.S. with result of program.                      M
*                                                                            *
* RETURN VALUE:                                                              M
*   void                                                                     M
*                                                                            *
* KEYWORDS:                                                                  M
*   AIsoShaderExit                                                           M
*****************************************************************************/
static void AIsoShaderExit(int ExitCode)
{
    exit(ExitCode);
}

#ifdef DEBUG

/*****************************************************************************
* DESCRIPTION:                                                               *
*    Dummy function to link at debugging time.                               *
*                                                                            *
* PARAMETERS:                                                                *
*                                                                            *
* RETURN VALUE:                                                              *
*   void                                                                     *
*                                                                            *
* KEYWORDS:                                                                  *
*****************************************************************************/
void DummyLinkCagdDebug(void)
{
    IritPrsrDbg();
}

#endif /* DEBUG */
