自学教程：C++ ELL_3V_SET函数代码示例

int_echoRayIntx_Cube(RAYINTX_ARGS(Cube)) {  echoPos_t t;  int ax, dir;  AIR_UNUSED(parm);  if (!_echoRayIntx_CubeSurf(&t, &ax, &dir,                             -1, 1,                             -1, 1,                             -1, 1, ray))     return AIR_FALSE;  intx->obj = (echoObject *)obj;  intx->t = t;  switch(ax) {  case 0: ELL_3V_SET(intx->norm, dir, 0, 0); break;  case 1: ELL_3V_SET(intx->norm, 0, dir, 0); break;  case 2: ELL_3V_SET(intx->norm, 0, 0, dir); break;  }  intx->face = ax + 3*(dir + 1)/2;  if (tstate->verbose) {    fprintf(stderr, "%s%s: ax = %d --> norm = (%g,%g,%g)/n",            _echoDot(tstate->depth), "_echoRayIntx_Cube", ax,            intx->norm[0], intx->norm[1], intx->norm[2]);  }  /* does NOT set u, v */  return AIR_TRUE;}

voidmakeSceneGlass2(limnCamera *cam, echoRTParm *parm, echoObject **sceneP) {  echoObject *cube, *rect;  echoObject *scene;  Nrrd *ntext;  echoPos_t matx[16];    *sceneP = scene = echoObjectNew(echoList);  ELL_3V_SET(cam->from, 0, 0, 100);  ELL_3V_SET(cam->at,   0, 0, 0);  ELL_3V_SET(cam->up,   0, 1, 0);  cam->uRange[0] = -1.0;  cam->uRange[1] = 1.0;  cam->vRange[0] = -1.0;  cam->vRange[1] = 1.0;  parm->jitterType = echoJitterNone;  parm->numSamples = 1;  parm->imgResU = 300;  parm->imgResV = 300;  parm->aperture = 0.0;  parm->renderLights = AIR_FALSE;  parm->shadow = 0.0;  parm->seedRand = AIR_FALSE;  parm->maxRecDepth = 10;  parm->mrR = 1.0;  parm->mrG = 0.0;  parm->mrB = 1.0;  ELL_4M_SET_SCALE(matx, 0.5, 0.5, 0.5);  cube = echoRoughSphere(80, 40, matx);  /*  cube = echoObjectNew(echoSphere);  echoSphereSet(cube, 0, 0, 0, 0.5);  */  echoMatterGlassSet(cube,                     1.0, 1.0, 1.0,                     1.33333, 0.0, 0.0);  echoObjectAdd(scene, cube);  nrrdLoad(ntext=nrrdNew(), "check.nrrd", NULL);    rect = echoObjectNew(echoRectangle);  printf("rect = %p/n", rect);  echoRectangleSet(rect,                         -1, -1, -0.51,                         2, 0, 0,                         0, 2, 0);  echoMatterPhongSet(rect, 1.0, 1.0, 1.0, 1.0,                     0.0, 1.0, 0.0, 40);  echoMatterTextureSet(rect, ntext);  echoObjectAdd(scene, rect);    /*  light = echoLightNew(echoLightDirectional);  echoLightDirectionalSet(light, 1, 1, 1, 0, 0, 1);  echoLightArrayAdd(lightArr, light);  */}

voidincTest(char *me, int num, baneRange *range) {  double *val, tmp, incParm[BANE_PARM_NUM], omin, omax, rmin, rmax;  baneInc *inc;  Nrrd *hist;  int i, j;  airSrand48();  val = (double*)malloc(num*sizeof(double));  /* from <http://www.itl.nist.gov/div898/handbook/index.htm>:     the standard dev of a uniform distribution between A and B is     sqrt((B-A)^2/12) */  for (j=0; j<num; j++) {    tmp = AIR_AFFINE(0.0, airDrand48(), 1.0, -1.0, 1.0);    /* val[j] = tmp*tmp*tmp; */    val[j] = tmp;  }  rmin = rmax = val[0];  for (j=0; j<num; j++) {    rmin = AIR_MIN(rmin, val[j]);    rmax = AIR_MAX(rmax, val[j]);  }  fprintf(stderr, "incTest: real min,max = %g,%g/n", rmin, rmax);    for (i=1; i<baneIncLast; i++) {    /* NOTE: THIS IS BROKEN !!! */    inc = baneIncNew(i, NULL, incParm);    printf("%s: inclusion %s ------/n", me, inc->name);    switch(i) {    case baneIncAbsolute:      ELL_3V_SET(incParm, -0.8, 1.5, AIR_NAN);      break;    case baneIncRangeRatio:      ELL_3V_SET(incParm, 0.99, AIR_NAN, AIR_NAN);      break;    case baneIncPercentile:      ELL_3V_SET(incParm, 1024, 10, AIR_NAN);      break;    case baneIncStdv:      ELL_3V_SET(incParm, 1.0, AIR_NAN, AIR_NAN);      break;    }    fprintf(stderr, "!%s: THIS IS BROKEN!!!/n", "incTest");    /*    if (inc->passA) {      for (j=0; j<num; j++)        inc->process[0](hist, val[j], incParm);    }    if (inc->passB) {      for (j=0; j<num; j++)        inc->process[1](hist, val[j], incParm);    }    inc->ans(&omin, &omax, hist, incParm, range);    */    printf(" --> (%g,%g)/n", omin, omax);  }  free(val);}

hooverContext *hooverContextNew() {    hooverContext *ctx;    ctx = (hooverContext *)calloc(1, sizeof(hooverContext));    if (ctx) {        ctx->cam = limnCameraNew();        ELL_3V_SET(ctx->volSize, 0, 0, 0);        ELL_3V_SET(ctx->volSpacing, AIR_NAN, AIR_NAN, AIR_NAN);        ctx->volCentering = hooverDefVolCentering;        ctx->shape = NULL;        ctx->imgSize[0] = ctx->imgSize[1] = 0;        ctx->imgCentering = hooverDefImgCentering;        ctx->user = NULL;        ctx->numThreads = 1;        ctx->workIdx = 0;        ctx->workMutex = NULL;        ctx->renderBegin = hooverStubRenderBegin;        ctx->threadBegin = hooverStubThreadBegin;        ctx->rayBegin = hooverStubRayBegin;        ctx->sample = hooverStubSample;        ctx->rayEnd = hooverStubRayEnd;        ctx->threadEnd = hooverStubThreadEnd;        ctx->renderEnd = hooverStubRenderEnd;    }    return(ctx);}

void_cap2xyz(double xyz[3], double ca, double cp, int version, int whole) {  double cl, cs, mean;  cs = 1 - ca;  cl = 1 - cs - cp;  mean = (cs + cp + cl)/3;  /*    xyz[0] = cs*0.333 + cl*1.0 + cp*0.5;    xyz[1] = cs*0.333 + cl*0.0 + cp*0.5;    xyz[2] = cs*0.333 + cl*0.0 + cp*0.0;    xyz[0] = AIR_AFFINE(0, ca, 1, 1.1*xyz[0], 0.86*xyz[0]);    xyz[1] = AIR_AFFINE(0, ca, 1, 1.1*xyz[1], 0.86*xyz[1]);    xyz[2] = AIR_AFFINE(0, ca, 1, 1.1*xyz[2], 0.86*xyz[2]);  */  if (whole) {    ELL_3V_SET(xyz,               AIR_AFFINE(0.0, 0.9, 1.0, mean, cl),               AIR_AFFINE(0.0, 0.9, 1.0, mean, cp),               AIR_AFFINE(0.0, 0.9, 1.0, mean, cs));    ELL_3V_SET(xyz, cl, cp, cs);  } else {    if (1 == version) {      ELL_3V_SET(xyz,                 (3 + 3*cl - cs)/6,                 (2 - 2*cl + cp)/6,                 2*cs/6);    } else {      ELL_3V_SET(xyz, 1, 1 - cl, cs);    }  }}

void_ra2t(Nrrd *nten, double rad, double angle,       double mRI[9], double mRF[9], double hack) {  double x, y, xyz[3], XX[3], YY[3], CC[3], EE[3], VV[3], tmp, mD[9], mT[9];  float *tdata;  int xi, yi, sx, sy;  sx = nten->axis[1].size;  sy = nten->axis[2].size;  x = rad*sin(AIR_PI*angle/180);  y = rad*cos(AIR_PI*angle/180);  xi = airIndexClamp(0.0, x, sqrt(3.0)/2.0, sx);  yi = airIndexClamp(0.0, y, 0.5, sy);  ELL_3V_SET(VV, 0, 3, 0);  ELL_3V_SET(EE, 1.5, 1.5, 0);  ELL_3V_SET(CC, 1, 1, 1);  ELL_3V_SUB(YY, EE, CC);  ELL_3V_SUB(XX, VV, EE);  ELL_3V_NORM(XX, XX, tmp);  ELL_3V_NORM(YY, YY, tmp);  ELL_3V_SCALE_ADD3(xyz, 1.0, CC, hack*x, XX, hack*y, YY);    ELL_3M_IDENTITY_SET(mD);  ELL_3M_DIAG_SET(mD, xyz[0], xyz[1], xyz[2]);  ELL_3M_IDENTITY_SET(mT);  ell_3m_post_mul_d(mT, mRI);  ell_3m_post_mul_d(mT, mD);  ell_3m_post_mul_d(mT, mRF);  tdata = (float*)(nten->data) + 7*(xi + sx*(yi + 1*sy));  tdata[0] = 1.0;  TEN_M2T(tdata, mT);}

/* XX != AB */static intlattABtoXX(int dstLatt, double *dstParm, const double *srcParm) {  double AA[2], BB[2], theta, phase, radi, area, len;  int ret = 0;  /* we have to reduce the DOF, which always starts with the same     loss of orientation information */  ELL_2V_COPY(AA, srcParm + 0);  ELL_2V_COPY(BB, srcParm + 2);  getToPosY(AA, BB);  switch(dstLatt) {  case rvaLattPRA:   /* AB -> PRA (loss off orientation) */    theta = atan2(BB[1], BB[0]);    phase = AIR_AFFINE(AIR_PI/2, theta, AIR_PI/3, 0.0, 1.0);    radi = _rvaLen2(BB)/_rvaLen2(AA);    area = _rvaLen2(AA)*BB[1];    ELL_3V_SET(dstParm, phase, radi, area);    break;  case rvaLattUVW:  /* AB -> UVW (loss of orientation) */    ELL_3V_SET(dstParm, BB[0], BB[1], AA[0]);    break;  case rvaLattXY:   /* AB -> XY (loss of orientation and scale) */    len = _rvaLen2(AA);    ELL_2V_SET(dstParm, BB[0]/len, BB[1]/len);    break;  default: ret = 1; break; /* unimplemented */  }  return ret;}

voidwashQtoM3(double m[9], double q[4]) {  double p[4], w, x, y, z, len;  ELL_4V_COPY(p, q);  len = ELL_4V_LEN(p);  ELL_4V_SCALE(p, 1.0/len, p);  w = p[0];  x = p[1];  y = p[2];  z = p[3];  /* mathematica work implies that we should be     setting ROW vectors here */  ELL_3V_SET(m+0,             1 - 2*(y*y + z*z),             2*(x*y - w*z),             2*(x*z + w*y));  ELL_3V_SET(m+3,             2*(x*y + w*z),             1 - 2*(x*x + z*z),             2*(y*z - w*x));  ELL_3V_SET(m+6,             2*(x*z - w*y),             2*(y*z + w*x),             1 - 2*(x*x + y*y));}

voidmakeSceneDOF(limnCamera *cam, echoRTParm *parm, echoScene *scene) {  echoObject *rect;  Nrrd *ntext;    ELL_3V_SET(cam->from, 6, 6, 20);  ELL_3V_SET(cam->at,   0, 0, 0);  ELL_3V_SET(cam->up,   0, 1, 0);  cam->uRange[0] = -3.3;  cam->uRange[1] = 3.3;  cam->vRange[0] = -3.3;  cam->vRange[1] = 3.3;  parm->jitterType = echoJitterJitter;  parm->numSamples = 4;  parm->imgResU = 300;  parm->imgResV = 300;  parm->aperture = 0.5;  parm->renderLights = AIR_FALSE;  parm->renderBoxes = AIR_FALSE;  parm->seedRand = AIR_FALSE;  parm->maxRecDepth = 10;  parm->shadow = 1.0;  nrrdLoad(ntext = nrrdNew(), "tmp.png", NULL);  rect = echoObjectNew(scene, echoTypeRectangle);  echoRectangleSet(rect,                   -0.5, 1.5, -3,                   2, 0, 0,                   0, -2, 0);  echoColorSet(rect, 1, 0.5, 0.5, 1);  echoMatterPhongSet(scene, rect, 1.0, 0.0, 0.0, 1);  echoMatterTextureSet(scene, rect, ntext);  echoObjectAdd(scene, rect);  rect = echoObjectNew(scene, echoTypeRectangle);  echoRectangleSet(rect,                   -1, 1, 0,                   2, 0, 0,                   0, -2, 0);  echoColorSet(rect, 0.5, 1, 0.5, 1);  echoMatterPhongSet(scene, rect, 1.0, 0.0, 0.0, 1);  echoMatterTextureSet(scene, rect, ntext);  echoObjectAdd(scene, rect);  rect = echoObjectNew(scene, echoTypeRectangle);  echoRectangleSet(rect,                   -1.5, 0.5, 3,                   2, 0, 0,                   0, -2, 0);  echoColorSet(rect, 0.5, 0.5, 1, 1);  echoMatterPhongSet(scene, rect, 1.0, 0.0, 0.0, 1);  echoMatterTextureSet(scene, rect, ntext);  echoObjectAdd(scene, rect);  return;}

void_coilKind7TensorFilterFinish(coil_t *delta, coil_t **iv3,                              double spacing[3],                             double parm[COIL_PARMS_NUM]) {  coil_t rspX, rspY, rspZ,    rspsqX, rspsqY, rspsqZ;  double eval[3], evec[9], tens[7], tengrad[21], grad[3], LL, KK,    cnd,    dmu1[7], dmu2[7], dskw[7], phi3[7];  rspX = AIR_CAST(coil_t, 1.0/spacing[0]); rspsqX = rspX*rspX;  rspY = AIR_CAST(coil_t, 1.0/spacing[1]); rspsqY = rspY*rspY;  rspZ = AIR_CAST(coil_t, 1.0/spacing[2]); rspsqZ = rspZ*rspZ;  TENS(tens, iv3);  TENGRAD(tengrad, iv3, rspX, rspY, rspZ);  tenEigensolve_d(eval, evec, tens);  tenInvariantGradientsK_d(dmu1, dmu2, dskw, tens, 0.000001);  tenRotationTangents_d(NULL, NULL, phi3, evec);  /* /midhat{/nabla} /mu_1 ----------------- */  ELL_3V_SET(grad,             TEN_T_DOT(dmu1, tengrad + 0*7),             TEN_T_DOT(dmu1, tengrad + 1*7),             TEN_T_DOT(dmu1, tengrad + 2*7));  LL = ELL_3V_DOT(grad,grad);  KK = parm[1]*parm[1];  cnd = _COIL_CONDUCT(LL, KK);  /* /midhat{/nabla} /mu_2 ----------------- */  ELL_3V_SET(grad,             TEN_T_DOT(dmu2, tengrad + 0*7),             TEN_T_DOT(dmu2, tengrad + 1*7),             TEN_T_DOT(dmu2, tengrad + 2*7));  LL = ELL_3V_DOT(grad,grad);  KK = parm[2]*parm[2];  cnd *= _COIL_CONDUCT(LL, KK);  /* /midhat{/nabla} /skw and twist! ----------------- */  ELL_3V_SET(grad,             TEN_T_DOT(dskw, tengrad + 0*7),             TEN_T_DOT(dskw, tengrad + 1*7),             TEN_T_DOT(dskw, tengrad + 2*7));  LL = ELL_3V_DOT(grad,grad);  ELL_3V_SET(grad,             TEN_T_DOT(phi3, tengrad + 0*7),             TEN_T_DOT(phi3, tengrad + 1*7),             TEN_T_DOT(phi3, tengrad + 2*7));  LL += ELL_3V_DOT(grad,grad);  KK = AIR_CAST(coil_t, parm[3]*parm[3]);  cnd *= _COIL_CONDUCT(LL, KK);   delta[0]= 0.0f;  delta[1]= AIR_CAST(coil_t, parm[0]*cnd*LAPL(iv3, 1, rspsqX, rspsqY, rspsqZ));  delta[2]= AIR_CAST(coil_t, parm[0]*cnd*LAPL(iv3, 2, rspsqX, rspsqY, rspsqZ));  delta[3]= AIR_CAST(coil_t, parm[0]*cnd*LAPL(iv3, 3, rspsqX, rspsqY, rspsqZ));  delta[4]= AIR_CAST(coil_t, parm[0]*cnd*LAPL(iv3, 4, rspsqX, rspsqY, rspsqZ));  delta[5]= AIR_CAST(coil_t, parm[0]*cnd*LAPL(iv3, 5, rspsqX, rspsqY, rspsqZ));  delta[6]= AIR_CAST(coil_t, parm[0]*cnd*LAPL(iv3, 6, rspsqX, rspsqY, rspsqZ));}

/* normalized gradients of k or r invariants, in XYZ space,   to help determine if path is really a loxodrome */voidkgrads(double grad[3][3], const double eval[3]) {  double rtz[3];  tenTripleConvertSingle_d(rtz, tenTripleTypeRThetaZ,                           eval, tenTripleTypeEigenvalue);  ELL_3V_SET(grad[0], cos(rtz[1]), sin(rtz[1]), 0);  ELL_3V_SET(grad[1], -sin(rtz[1]), cos(rtz[1]), 0);  ELL_3V_SET(grad[2], 0, 0, 1);}

voidwheelGeomToABC(double ABC[3], double center, double radius, double angle) {  double geom[3], x[3], yroot[3];  ELL_3V_SET(geom, center, radius, angle);  wheelGeomToRoot(x, yroot, geom);  ELL_3V_SET(ABC,             -x[0] - x[1] - x[2],             x[0]*x[1] + x[1]*x[2] + x[0]*x[2],             -x[0]*x[1]*x[2]);  return;}

voidechoRectangleSet(echoObject *rect,                 echoPos_t ogx, echoPos_t ogy, echoPos_t ogz,                 echoPos_t e0x, echoPos_t e0y, echoPos_t e0z,                 echoPos_t e1x, echoPos_t e1y, echoPos_t e1z) {  if (rect && echoTypeRectangle == rect->type) {    ELL_3V_SET(RECTANGLE(rect)->origin, ogx, ogy, ogz);    ELL_3V_SET(RECTANGLE(rect)->edge0, e0x, e0y, e0z);    ELL_3V_SET(RECTANGLE(rect)->edge1, e1x, e1y, e1z);  }  return;}

dyeColor *dyeColorInit(dyeColor *col) {  if (col) {    ELL_3V_SET(col->val[0], AIR_NAN, AIR_NAN, AIR_NAN);    ELL_3V_SET(col->val[1], AIR_NAN, AIR_NAN, AIR_NAN);    col->xWhite = col->yWhite = AIR_NAN;    col->spc[0] = dyeSpaceUnknown;    col->spc[1] = dyeSpaceUnknown;    col->ii = 0;  }  return col;}

voidechoTriangleSet(echoObject *tri,                echoPos_t xx0, echoPos_t yy0, echoPos_t zz0,                echoPos_t xx1, echoPos_t yy1, echoPos_t zz1,                echoPos_t xx2, echoPos_t yy2, echoPos_t zz2) {  if (tri && echoTypeTriangle == tri->type) {    ELL_3V_SET(TRIANGLE(tri)->vert[0], xx0, yy0, zz0);    ELL_3V_SET(TRIANGLE(tri)->vert[1], xx1, yy1, zz1);    ELL_3V_SET(TRIANGLE(tri)->vert[2], xx2, yy2, zz2);  }  return;}

voidgageShapeReset(gageShape *shape) {    if (shape) {    shape->defaultCenter = gageDefDefaultCenter;    ELL_3V_SET(shape->size, 0, 0, 0);    shape->center = nrrdCenterUnknown;    shape->fromOrientation = AIR_FALSE;    ELL_3V_SET(shape->spacing, AIR_NAN, AIR_NAN, AIR_NAN);    ELL_3V_SET(shape->volHalfLen, AIR_NAN, AIR_NAN, AIR_NAN);    ELL_3V_SET(shape->voxLen, AIR_NAN, AIR_NAN, AIR_NAN);  }  return;}

voidmakeSceneRainLights(limnCamera *cam, echoRTParm *parm, echoScene *scene) {  echoObject *sphere, *rect;  int i, N;  echoPos_t w;  float r, g, b;  ELL_3V_SET(cam->from, 2.5, 0, 5);  ELL_3V_SET(cam->at,   0, 0, 0);  ELL_3V_SET(cam->up,   0, 0, 1);  cam->uRange[0] = -1.7;  cam->uRange[1] = 1.7;  cam->vRange[0] = -1.7;  cam->vRange[1] = 1.7;  parm->jitterType = echoJitterJitter;  parm->numSamples = 36;  parm->imgResU = 1000;  parm->imgResV = 1000;  parm->numSamples = 16;  parm->imgResU = 200;  parm->imgResV = 200;  parm->aperture = 0.0;  parm->renderLights = AIR_TRUE;  parm->shadow = 0.0;  ELL_3V_SET(scene->bkgr, 0.1, 0.1, 0.1);  /* create scene */  sphere = echoObjectNew(scene, echoTypeSphere);  echoSphereSet(sphere, 0, 0, 0, 1.0);  echoColorSet(sphere, 1.0, 1.0, 1.0, 1.0);  echoMatterPhongSet(scene, sphere, 0.02, 0.2, 1.0, 400);  echoObjectAdd(scene, sphere);  N = 8;  w = 1.7/N;  for (i=0; i<N; i++) {    rect = echoObjectNew(scene, echoTypeRectangle);    echoRectangleSet(rect,                     w/2, AIR_AFFINE(0, i, N-1, -1-w/2, 1-w/2), 1.5,                     0, w, 0,                     w, 0, 0);    _dyeHSVtoRGB(&r, &g, &b, AIR_AFFINE(0, i, N, 0.0, 1.0), 1.0, 1.0);    echoColorSet(rect, r, g, b, 1);    echoMatterLightSet(scene, rect, 1, 0);    echoObjectAdd(scene, rect);  }}

miteUser *miteUserNew() {  miteUser *muu;  int i;  muu = (miteUser *)calloc(1, sizeof(miteUser));  if (!muu)    return NULL;  muu->umop = airMopNew();  muu->nsin = NULL;  muu->nvin = NULL;  muu->ntin = NULL;  muu->ntxf = NULL;              /* managed by user (with miter: hest) */  muu->nout = NULL;              /* managed by user (with miter: hest) */  muu->debug = NULL;  muu->debugArr = NULL;  muu->ndebug = NULL;            /* not allocated until the debug pixel                                    is rendered, see miteRayBegin */  muu->ntxfNum = 0;  muu->shadeStr[0] = 0;  muu->normalStr[0] = 0;  for (i=0; i<MITE_RANGE_NUM; i++) {    muu->rangeInit[i] = 1.0;  }  muu->normalSide = miteDefNormalSide;  muu->refStep = miteDefRefStep;  muu->rayStep = AIR_NAN;  muu->opacMatters = miteDefOpacMatters;  muu->opacNear1 = miteDefOpacNear1;  muu->hctx = hooverContextNew();  ELL_3V_SET(muu->fakeFrom, AIR_NAN, AIR_NAN, AIR_NAN);  ELL_3V_SET(muu->vectorD, 0, 0, 0);  airMopAdd(muu->umop, muu->hctx, (airMopper)hooverContextNix, airMopAlways);  for (i=0; i<GAGE_KERNEL_NUM; i++) {    muu->ksp[i] = NULL;  }  muu->gctx0 = gageContextNew();  airMopAdd(muu->umop, muu->gctx0, (airMopper)gageContextNix, airMopAlways);  gageParmSet(muu->gctx0, gageParmRequireAllSpacings, AIR_FALSE);  muu->lit = limnLightNew();  airMopAdd(muu->umop, muu->lit, (airMopper)limnLightNix, airMopAlways);  muu->normalSide = miteDefNormalSide;  muu->verbUi = muu->verbVi = -1;  muu->rendTime = 0;  muu->sampRate = 0;  return muu;}

echoRTParm *echoRTParmNew(void) {  echoRTParm *parm;    parm = (echoRTParm *)calloc(1, sizeof(echoRTParm));  if (parm) {    parm->jitterType = echoJitterNone;    parm->reuseJitter = AIR_FALSE;    parm->permuteJitter = AIR_TRUE;    parm->textureNN = AIR_TRUE;    parm->numSamples = 1;    parm->imgResU = parm->imgResV = 256;    parm->maxRecDepth = 5;    parm->renderLights = AIR_TRUE;    parm->renderBoxes = AIR_FALSE;    parm->seedRand = AIR_TRUE;    parm->sqNRI = 15;    parm->numThreads = 1;    parm->sqTol = 0.0001;    parm->aperture = 0.0;     /* pinhole camera by default */    parm->timeGamma = 6.0;    parm->boxOpac = 0.2f;    parm->shadow = 1.0;    parm->glassC = 3;    ELL_3V_SET(parm->maxRecCol, 1.0, 0.0, 1.0);  }  return parm;}

int_tenFiberIntegrate_RK4(tenFiberContext *tfx, double forwDir[3]) {  double loc[3], k1[3], k2[3], k3[3], k4[3], c1, c2, c3, c4, h;  h = tfx->stepSize;  c1 = h/6.0; c2 = h/3.0; c3 = h/3.0; c4 = h/6.0;  _tenFiberStep[tfx->fiberType](tfx, k1);  ELL_3V_SCALE_ADD2(loc, 1, tfx->wPos, 0.5*h, k1);  if (_tenFiberProbe(tfx, loc)) return 1;  _tenFiberStep[tfx->fiberType](tfx, k2);  ELL_3V_SCALE_ADD2(loc, 1, tfx->wPos, 0.5*h, k2);  if (_tenFiberProbe(tfx, loc)) return 1;  _tenFiberStep[tfx->fiberType](tfx, k3);  ELL_3V_SCALE_ADD2(loc, 1, tfx->wPos, h, k3);  if (_tenFiberProbe(tfx, loc)) return 1;  _tenFiberStep[tfx->fiberType](tfx, k4);  ELL_3V_SET(forwDir,             c1*k1[0] + c2*k2[0] + c3*k3[0] + c4*k4[0],             c1*k1[1] + c2*k2[1] + c3*k3[1] + c4*k4[1],             c1*k1[2] + c2*k2[2] + c3*k3[2] + c4*k4[2]);    return 0;}

static voidbaryFind(double bcoord[3], const double uvp[2],         const double uv0[2],         const double uv1[2],         const double uv2[2]) {  double mat[9], a, a01, a02, a12;  ELL_3M_SET(mat,             uv0[0], uv0[1], 1,             uv1[0], uv1[1], 1,             uvp[0], uvp[1], 1);  a01 = ELL_3M_DET(mat); a01 = AIR_ABS(a01);  ELL_3M_SET(mat,             uv0[0], uv0[1], 1,             uv2[0], uv2[1], 1,             uvp[0], uvp[1], 1);  a02 = ELL_3M_DET(mat); a02 = AIR_ABS(a02);  ELL_3M_SET(mat,             uv1[0], uv1[1], 1,             uv2[0], uv2[1], 1,             uvp[0], uvp[1], 1);  a12 = ELL_3M_DET(mat); a12 = AIR_ABS(a12);  a = a01 + a02 + a12;  ELL_3V_SET(bcoord, a12/a, a02/a, a01/a);  return;}

/*** for each pvl: pvl's query --> pvl's needD*/void_gagePvlNeedDUpdate(gageContext *ctx) {  char me[]="_gagePvlNeedDUpdate";  gagePerVolume *pvl;  int que, needD[3];  unsigned int pvlIdx;  if (ctx->verbose) fprintf(stderr, "%s: hello/n", me);  for (pvlIdx=0; pvlIdx<ctx->pvlNum; pvlIdx++) {    pvl = ctx->pvl[pvlIdx];    if (pvl->flag[gagePvlFlagQuery]) {      ELL_3V_SET(needD, 0, 0, 0);      que = pvl->kind->itemMax+1;      do {        que--;        if (GAGE_QUERY_ITEM_TEST(pvl->query, que)) {          needD[pvl->kind->table[que].needDeriv] = 1;        }      } while (que);      if (!ELL_3V_EQUAL(needD, pvl->needD)) {        if (ctx->verbose) {          fprintf(stderr, "%s: updating pvl[%d]'s needD to (%d,%d,%d)/n",                  me, pvlIdx, needD[0], needD[1], needD[2]);        }        ELL_3V_COPY(pvl->needD, needD);        pvl->flag[gagePvlFlagNeedD] = AIR_TRUE;      }    }  }  if (ctx->verbose) fprintf(stderr, "%s: bye/n", me);  return;}

echoPos_t_echo_SuperquadZ_lvg(echoPos_t grad[3],                     echoPos_t x, echoPos_t y, echoPos_t z,                     echoPos_t A, echoPos_t B) {  echoPos_t R, xxa, yya, zzb, larg;  echoPos_t ret;  xxa = pow(x*x, 1/A);  yya = pow(y*y, 1/A);  zzb = pow(z*z, 1/B);  R = pow(xxa + yya, 1-(A/B))*zzb;  ELL_3V_SET(grad,             2*xxa/(B*x*(xxa + yya + R)),             2*yya/(B*y*(xxa + yya + R)),             2/(B*z*(1 + pow(xxa + yya, A/B)/zzb)));  larg = pow(xxa + yya, A/B) + zzb;  ret = larg > 0 ? log(larg) : ECHO_POS_MIN;  /*  if (!AIR_EXISTS(ret)) {    fprintf(stderr, "_echo_SuperquadZ_lvg: PANIC/n");    fprintf(stderr, "x = %g, y = %g, z = %g, A = %g, B = %g/n",            x, y, z, A, B);    fprintf(stderr, "pow(%g*%g = %g, %g) = %g/n",             x, x, x*x, 1/A, xxa);    fprintf(stderr, "pow(%g*%g = %g, %g) = %g/n",             y, y, y*y, 1/A, yya);    fprintf(stderr, "log(pow(%g, %g) = %g + %g) = %g/n",            xxa + yya, A/B, pow(xxa + yya, A/B), zzb, ret);    exit(0);  }  */  return ret;}

voidtenTripleCalcSingle_d(double dst[3], int ttype, double ten[7]) {  double eval[3];  /* in time this can become more efficient ... */  switch (ttype) {  case tenTripleTypeEigenvalue:    tenEigensolve_d(dst, NULL, ten);    break;  case tenTripleTypeMoment:  case tenTripleTypeXYZ:  case tenTripleTypeRThetaZ:  case tenTripleTypeRThetaPhi:  case tenTripleTypeK:  case tenTripleTypeJ:  case tenTripleTypeWheelParm:    tenEigensolve_d(eval, NULL, ten);    tenTripleConvertSingle_d(dst, ttype, eval, tenTripleTypeEigenvalue);    break;  case tenTripleTypeR:    dst[0] = sqrt(_tenAnisoTen_d[tenAniso_S](ten));    dst[1] = _tenAnisoTen_d[tenAniso_FA](ten);    dst[2] = _tenAnisoTen_d[tenAniso_Mode](ten);    break;  default:    /* what on earth? */    ELL_3V_SET(dst, AIR_NAN, AIR_NAN, AIR_NAN);  }  return;}

/*** all pvls' needD --> ctx's needD*/void_gageNeedDUpdate(gageContext *ctx) {  char me[]="_gageNeedDUpdate";  gagePerVolume *pvl;  int needD[3];  unsigned int pvlIdx;  if (ctx->verbose) fprintf(stderr, "%s: hello/n", me);  ELL_3V_SET(needD, 0, 0, 0);  for (pvlIdx=0; pvlIdx<ctx->pvlNum; pvlIdx++) {    pvl = ctx->pvl[pvlIdx];    needD[0] |= pvl->needD[0];    needD[1] |= pvl->needD[1];    needD[2] |= pvl->needD[2];  }  if (!ELL_3V_EQUAL(needD, ctx->needD)) {    if (ctx->verbose) {      fprintf(stderr, "%s: updating ctx's needD to (%d,%d,%d)/n",              me, needD[0], needD[1], needD[2]);    }    ELL_3V_COPY(ctx->needD, needD);    ctx->flag[gageCtxFlagNeedD] = AIR_TRUE;  }  if (ctx->verbose) fprintf(stderr, "%s: bye/n", me);  return;}

voidrgrads(double grad[3][3], const double eval[3]) {  double rtp[3];  tenTripleConvertSingle_d(rtp, tenTripleTypeRThetaPhi,                           eval, tenTripleTypeEigenvalue);  ELL_3V_SET(grad[0],             cos(rtp[1])*sin(rtp[2]),             sin(rtp[1])*sin(rtp[2]),             cos(rtp[2]));  ELL_3V_SET(grad[1], -sin(rtp[1]), cos(rtp[1]), 0);  ELL_3V_SET(grad[2],             cos(rtp[1])*cos(rtp[2]),             sin(rtp[1])*cos(rtp[2]),             -sin(rtp[2]));}

voidmakeSceneBVH(limnCamera *cam, echoRTParm *parm, echoObject **sceneP) {  echoObject *sphere;  int i, N;  float r, g, b;  echoObject *scene;  double time0, time1;    *sceneP = scene = echoObjectNew(echoList);  ELL_3V_SET(cam->from, 9, 6, 0);  ELL_3V_SET(cam->at,   0, 0, 0);  ELL_3V_SET(cam->up,   0, 0, 1);  cam->uRange[0] = -3;  cam->uRange[1] = 3;  cam->vRange[0] = -3;  cam->vRange[1] = 3;  parm->jitterType = echoJitterNone;  parm->numSamples = 1;  parm->imgResU = 500;  parm->imgResV = 500;  parm->aperture = 0.0;  parm->renderLights = AIR_TRUE;  parm->renderBoxes = AIR_FALSE;  parm->seedRand = AIR_FALSE;  parm->maxRecDepth = 10;  parm->shadow = 0.0;  N = 1000000;  airArrayLenSet(LIST(scene)->objArr, N);  for (i=0; i<N; i++) {    sphere = echoObjectNew(echoSphere);    echoSphereSet(sphere,                  4*airDrandMT()-2, 4*airDrandMT()-2, 4*airDrandMT()-2, 0.005);    _dyeHSVtoRGB(&r, &g, &b, AIR_AFFINE(0, i, N, 0.0, 1.0), 1.0, 1.0);    echoMatterPhongSet(sphere, r, g, b, 1.0,                       1.0, 0.0, 0.0, 50);    LIST(scene)->obj[i] = sphere;  }  time0 = airTime();  *sceneP = scene = echoListSplit3(scene, 8);  time1 = airTime();  printf("BVH build time = %g seconds/n", time1 - time0);}

/********* ell_3m2sub_eigenvalues_d**** for doing eigensolve of the upper-left 2x2 matrix sub-matrix of a** 3x3 matrix.  The other entries are assumed to be zero.  A 0 root is** put last (in eval[2]), possibly in defiance of the usual eigenvalue** ordering.*/intell_3m2sub_eigenvalues_d(double eval[3], const double _m[9]) {  double A, B, m[4], D, Dsq, eps=1.0E-11;  int roots;  /* static const char me[]="ell_3m2sub_eigenvalues_d"; */  m[0] = _m[0];  m[1] = _m[1];  m[2] = _m[3];  m[3] = _m[4];  /* cubic characteristic equation is L^3 + A*L^2 + B*L = 0 */  A = -m[0] - m[3];  B = m[0]*m[3] - m[1]*m[2];  Dsq = A*A - 4*B;  /*  fprintf(stderr, "!%s: m = {{%f,%f},{%f,%f}} -> A=%f B=%f Dsq=%.17f %s 0 (%.17f)/n", me,          m[0], m[1], m[2], m[3], A, B, Dsq,          (Dsq > 0 ? ">" : (Dsq < 0 ? "<" : "==")), eps);  fprintf(stderr, "!%s: Dsq = /n", me);  airFPFprintf_d(stderr, Dsq);  fprintf(stderr, "!%s: eps = /n", me);  airFPFprintf_d(stderr, eps);  ell_3m_print_d(stderr, _m);  */  if (Dsq > eps) {    D = sqrt(Dsq);    eval[0] = (-A + D)/2;    eval[1] = (-A - D)/2;    eval[2] = 0;    roots = ell_cubic_root_three;  } else if (Dsq < -eps) {    /* no quadratic roots; only the implied zero */    ELL_3V_SET(eval, AIR_NAN, AIR_NAN, 0);    roots = ell_cubic_root_single;  } else {    /* a quadratic double root */    ELL_3V_SET(eval, -A/2, -A/2, 0);    roots = ell_cubic_root_single_double;  }  /*  fprintf(stderr, "!%s: Dsq=%f, roots=%d (%f %f %f)/n",          me, Dsq, roots, eval[0], eval[1], eval[2]);  */  return roots;}

void_limnOptsPSDefaults(limnOptsPS *ps) {  ps->lineWidth[limnEdgeTypeUnknown] = AIR_NAN;  ps->lineWidth[limnEdgeTypeBackFacet] = 0.0;  ps->lineWidth[limnEdgeTypeBackCrease] = 0.0;  ps->lineWidth[limnEdgeTypeContour] = 2.0;  ps->lineWidth[limnEdgeTypeFrontCrease] = 1.0;  ps->lineWidth[limnEdgeTypeFrontFacet] = 0.0;  ps->lineWidth[limnEdgeTypeBorder] = 1.0;  ps->lineWidth[limnEdgeTypeLone] = 1.0;  ps->creaseAngle = 46;  ps->showpage = AIR_FALSE;  ps->wireFrame = AIR_FALSE;  ps->noBackground = AIR_FALSE;  ELL_3V_SET(ps->bg, 1, 1, 1);  ELL_3V_SET(ps->edgeColor, 0, 0, 0);}

voiddyeXYZtoRGB(float *R, float *G, float *B,            float  X, float  Y, float  Z) {  float in[3], out[3];    ELL_3V_SET(in, X, Y, Z);  ELL_3MV_MUL(out, dyeXYZtoRGBMatx, in);  ELL_3V_GET(*R, *G, *B, out);  return;}