void
plotDots(plotPS *pps, plotParm *pparm, Nrrd **ndata, int nidx) {
  int ii, npts;
  double xx, yy, orad, irad, *data, val;

  if (!( pparm->dotDiameter[nidx] > 0 )) {
    return;
  }

  fprintf(pps->file, "gsave\n");
  fprintf(pps->file, "newpath\n");
  plotWidth(pps, pparm, 0);
  data = (double *)(ndata[nidx]->data);
  npts = ndata[nidx]->axis[1].size;
  orad = pparm->dotDiameter[nidx]/2;
  irad = pparm->dotInnerDiameterFraction*orad;
  for (ii=0; ii<npts; ii++) {
    val = data[ii];
    xx = AIR_AFFINE(0, ii, npts-1,
                    ndata[nidx]->axis[1].min, ndata[nidx]->axis[1].max);
    xx = AIR_AFFINE(pparm->dbox[0], xx, pparm->dbox[2],
                    pps->bbox[0], pps->bbox[2]);
    yy = AIR_AFFINE(pparm->dbox[1], val, pparm->dbox[3],
                    pps->bbox[1], pps->bbox[3]);
    plotGray(pps, pparm, pparm->dotGray[nidx]);
    fprintf(pps->file, "%g %g %g 0 360 arc closepath fill\n", 
            PPS_X(xx), PPS_Y(yy), PPS_S(orad));
    if (irad) {
      plotGray(pps, pparm, 1.0);
      fprintf(pps->file, "%g %g %g 0 360 arc closepath fill\n", 
              PPS_X(xx), PPS_Y(yy), PPS_S(irad));
    }
  }
  fprintf(pps->file, "grestore\n");
}

/*
******** tenGradientRandom
**
** generates num random unit vectors of type double
*/
int
tenGradientRandom(Nrrd *ngrad, unsigned int num, unsigned int seed) {
  static const char me[]="tenGradientRandom";
  double *grad, len;
  unsigned int gi;

  if (nrrdMaybeAlloc_va(ngrad, nrrdTypeDouble, 2,
                        AIR_CAST(size_t, 3), AIR_CAST(size_t, num))) {
    biffMovef(TEN, NRRD, "%s: couldn't allocate output", me);
    return 1;
  }
  airSrandMT(seed);
  grad = AIR_CAST(double*, ngrad->data);
  for (gi=0; gi<num; gi++) {
    do {
      grad[0] = AIR_AFFINE(0, airDrandMT(), 1, -1, 1);
      grad[1] = AIR_AFFINE(0, airDrandMT(), 1, -1, 1);
      grad[2] = AIR_AFFINE(0, airDrandMT(), 1, -1, 1);
      len = ELL_3V_LEN(grad);
    } while (len > 1 || !len);
    ELL_3V_SCALE(grad, 1.0/len, grad);
    grad += 3;
  }
  return 0;
}

void
plotGraph(plotPS *pps, plotParm *pparm, Nrrd **ndata, int nidx) {
  int ii, npts;
  double xx, yy, *data, val;

  if (!( pparm->graphThick[nidx] > 0 )) {
    return;
  }

  data = (double *)(ndata[nidx]->data);
  npts = ndata[nidx]->axis[1].size;
  plotGray(pps, pparm, pparm->graphGray[nidx]);
  fprintf(pps->file, "%g W\n", pps->psc*pparm->graphThick[nidx]);
  for (ii=0; ii<npts; ii++) {
    val = data[ii];
    xx = AIR_AFFINE(0, ii, npts-1,
                    ndata[nidx]->axis[1].min, ndata[nidx]->axis[1].max);
    xx = AIR_AFFINE(pparm->dbox[0], xx, pparm->dbox[2],
                    pps->bbox[0], pps->bbox[2]);
    yy = AIR_AFFINE(pparm->dbox[1], val, pparm->dbox[3],
                    pps->bbox[1], pps->bbox[3]);
    fprintf(pps->file, "%g %g %s\n", PPS_X(xx), PPS_Y(yy),
            ii ? "L" : "M");
  }
  fprintf(pps->file, "S\n");
}

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);
    }
  }
}

/*
******** nrrdAxisInfoPosRange()
**
** given a nrrd, an axis, and two (floating point) index space positions,
** return the range of positions implied the axis's min, max, and center
** The opposite of nrrdAxisIdxRange()
*/
void
nrrdAxisInfoPosRange(double *loP, double *hiP,
                     const Nrrd *nrrd, unsigned int ax, 
                     double loIdx, double hiIdx) {
  int center, flip = 0;
  size_t size;
  double min, max, tmp;

  if (!( loP && hiP && nrrd && ax <= nrrd->dim-1 )) {
    *loP = *hiP = AIR_NAN;
    return;
  }
  center = _nrrdCenter(nrrd->axis[ax].center);
  min = nrrd->axis[ax].min;
  max = nrrd->axis[ax].max;
  size = nrrd->axis[ax].size;

  if (loIdx > hiIdx) {
    flip = 1;
    tmp = loIdx; loIdx = hiIdx; hiIdx = tmp;
  }
  if (nrrdCenterCell == center) {
    *loP = AIR_AFFINE(0, loIdx, size, min, max);
    *hiP = AIR_AFFINE(0, hiIdx+1, size, min, max);
  } else {
    *loP = AIR_AFFINE(0, loIdx, size-1, min, max);
    *hiP = AIR_AFFINE(0, hiIdx, size-1, min, max);
  }
  if (flip) {
    tmp = *loP; *loP = *hiP; *hiP = tmp;
  }

  return;
}

void
makeSceneRainLights(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);
  }

}

echoObject *
echoRoughSphereNew(echoScene *scene, int theRes, int phiRes, echoPos_t *matx) {
  echoObject *trim;
  echoPos_t *_pos, *pos, tmp[3];
  int *_vert, *vert, thidx, phidx, n;
  echoPos_t th, ph;

  trim = echoObjectNew(scene, echoTypeTriMesh);
  TRIMESH(trim)->numV = 2 + (phiRes-1)*theRes;
  TRIMESH(trim)->numF = (2 + 2*(phiRes-2))*theRes;

  _pos = pos = (echoPos_t *)calloc(3*TRIMESH(trim)->numV, sizeof(echoPos_t));
  _vert = vert = (int *)calloc(3*TRIMESH(trim)->numF, sizeof(int));

  ELL_3V_SET(tmp, 0, 0, 1); _echoPosSet(pos, matx, tmp); pos += 3;
  for (phidx=1; phidx<phiRes; phidx++) {
    ph = AIR_AFFINE(0, phidx, phiRes, 0.0, AIR_PI);
    for (thidx=0; thidx<theRes; thidx++) {
      th = AIR_AFFINE(0, thidx, theRes, 0.0, 2*AIR_PI);
      ELL_3V_SET(tmp, cos(th)*sin(ph), sin(th)*sin(ph), cos(ph));
      _echoPosSet(pos, matx, tmp); pos += 3;
    }
  }
  ELL_3V_SET(tmp, 0, 0, -1); _echoPosSet(pos, matx, tmp);

  for (thidx=0; thidx<theRes; thidx++) {
    n = AIR_MOD(thidx+1, theRes);
    ELL_3V_SET(vert, 0, 1+thidx, 1+n); vert += 3;
  }
  for (phidx=0; phidx<phiRes-2; phidx++) {
    for (thidx=0; thidx<theRes; thidx++) {
      n = AIR_MOD(thidx+1, theRes);
      ELL_3V_SET(vert, 1+phidx*theRes+thidx, 1+(1+phidx)*theRes+thidx,
                 1+phidx*theRes+n); vert += 3;
      ELL_3V_SET(vert, 1+(1+phidx)*theRes+thidx, 1+(1+phidx)*theRes+n, 
                 1+phidx*theRes+n); vert += 3;
    }
  }
  for (thidx=0; thidx<theRes; thidx++) {
    n = AIR_MOD(thidx+1, theRes);
    ELL_3V_SET(vert, 1+(phiRes-2)*theRes+thidx, TRIMESH(trim)->numV-1,
               1+(phiRes-2)*theRes+n); 
    vert += 3;
  }

  echoTriMeshSet(trim, TRIMESH(trim)->numV, _pos, TRIMESH(trim)->numF, _vert);
  return(trim);
}

/* Calculate the Q-ball profile from DWIs */
void
_tenQball(const double b, const int gradcount, const double svals[],
          const double grads[], double qvals[] ) {
  /* Not an optimal Q-ball implementation! (Taken from submission to
     MICCAI 2006) Should be solved analytically in the future,
     implemented from recent papers. */
  int i,j;
  double d, dist, weight, min, max;

  AIR_UNUSED(b);
  min = max = svals[1] / svals[0];
  for( i = 0; i < gradcount; i++ ) {
    d = svals[i+1] / svals[0];
    if( d > max )
      max = d;
    else if( d < min )
      min = d;
  }

  for( i = 0; i < gradcount; i++ ) {
    qvals[i] = 0;
    for( j = 0; j < gradcount; j++ ) {
      d = AIR_AFFINE( min, svals[j+1] / svals[0], max, 0,1 );
      dist = ELL_3V_DOT(grads + 3*i, grads + 3*j);
      dist = AIR_ABS(dist);
      weight = cos( 0.5 * AIR_PI * dist );
      qvals[i] += d * weight*weight*weight*weight;
    }
  }
  return;
}

/* XX != AB */
static int
lattABtoXX(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;
}

/* XX != AB */
static int
lattXXtoAB(double *dstParm, int srcLatt, const double *srcParm) {
  double AA[2], BB[2], area, theta, radi, scl;
  int ret = 0;

  switch(srcLatt) {
  case rvaLattPRA:  /* PRA -> AB */
    area = AIR_ABS(srcParm[2]);
    theta = AIR_AFFINE(0, srcParm[0], 1, AIR_PI/2, AIR_PI/3);
    radi = srcParm[1];
    ELL_2V_SET(AA, 1.0, 0.0);
    ELL_2V_SET(BB, radi*cos(theta), radi*sin(theta));
    /* area from AA and BB is BB[1], but need to scale
       these to get to requested area */
    scl = sqrt(area/BB[1]);
    ELL_4V_SET(dstParm, scl*AA[0], scl*AA[1], scl*BB[0], scl*BB[1]);
    break;
  case rvaLattAB:  /* UVW -> AB */
    ELL_4V_SET(dstParm, srcParm[2], 0.0, srcParm[0], srcParm[1]);
    break;
  case rvaLattXY:  /* XY -> AB */
    ELL_4V_SET(dstParm, 1.0, 0.0, srcParm[0], srcParm[1]);
    break;
  default: ret = 1; break; /* unimplemented */
  }
  return ret;
}

/* ----------------------------------------------------------------
** ------------------------------ SPRING --------------------------
** ----------------------------------------------------------------
** 1 parms:
** parm[0]: width of pull region (beyond 1.0)
**
** learned: "1/2" is not 0.5 !!!!!
*/
double
_pullEnergySpringEval(double *frc, double dist, const double *parm) {
  /* char me[]="_pullEnergySpringEval"; */
  double enr, xx, pull;

  pull = parm[0];
  /* support used to be [0,1 + pull], but now is scrunched to [0,1],
     so hack "dist" to match old parameterization */
  dist = AIR_AFFINE(0, dist, 1, 0, 1+pull);
  xx = dist - 1.0;
  if (xx > pull) {
    enr = 0;
    *frc = 0;
  } else if (xx > 0) {
    enr = xx*xx*(xx*xx/(4*pull*pull) - 2*xx/(3*pull) + 1.0/2.0);
    *frc = xx*(xx*xx/(pull*pull) - 2*xx/pull + 1);
  } else {
    enr = xx*xx/2;
    *frc = xx;
  }
  /*
  if (!AIR_EXISTS(ret)) {
    fprintf(stderr, "!%s: dist=%g, pull=%g, blah=%d --> ret=%g\n",
            me, dist, pull, blah, ret);
  }
  */
  return enr;
}

int
limnSplineSample(Nrrd *nout, limnSpline *spline,
                 double minT, size_t M, double maxT) {
  char me[]="limnSplineSample", err[BIFF_STRLEN];
  airArray *mop;
  Nrrd *ntt;
  double *tt;
  size_t I;

  if (!(nout && spline)) {
    sprintf(err, "%s: got NULL pointer", me);
    biffAdd(LIMN, err); return 1;
  }
  mop = airMopNew();
  airMopAdd(mop, ntt=nrrdNew(), (airMopper)nrrdNuke, airMopAlways);
  if (nrrdMaybeAlloc_va(ntt, nrrdTypeDouble, 1,
                        M)) {
    sprintf(err, "%s: trouble allocating tmp nrrd", me);
    biffMove(LIMN, err, NRRD); airMopError(mop); return 1;
  }
  tt = (double*)(ntt->data);
  for (I=0; I<M; I++) {
    tt[I] = AIR_AFFINE(0, I, M-1, minT, maxT);
  }
  if (limnSplineNrrdEvaluate(nout, spline, ntt)) {
    sprintf(err, "%s: trouble", me);
    biffAdd(LIMN, err); airMopError(mop); return 1;
  }
  airMopOkay(mop);
  return 0;
}

void
incTest(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);
}

/*
******** tenGradientJitter
**
** moves all gradients by amount dist on tangent plane, in a random
** direction, and then renormalizes. The distance is a fraction
** of the ideal edge length (via tenGradientIdealEdge)
*/
int
tenGradientJitter(Nrrd *nout, const Nrrd *nin, double dist) {
  static const char me[]="tenGradientJitter";
  double *grad, perp0[3], perp1[3], len, theta, cc, ss, edge;
  unsigned int gi, num;

  if (nrrdConvert(nout, nin, nrrdTypeDouble)) {
    biffMovef(TEN, NRRD, "%s: trouble converting input to double", me);
    return 1;
  }
  if (tenGradientCheck(nout, nrrdTypeDouble, 3)) {
    biffAddf(TEN, "%s: didn't get valid gradients", me);
    return 1;
  }
  grad = AIR_CAST(double*, nout->data);
  num = AIR_UINT(nout->axis[1].size);
  /* HEY: possible confusion between single and not */
  edge = tenGradientIdealEdge(num, AIR_FALSE);
  for (gi=0; gi<num; gi++) {
    ELL_3V_NORM(grad, grad, len);
    ell_3v_perp_d(perp0, grad);
    ELL_3V_CROSS(perp1, perp0, grad);
    theta = AIR_AFFINE(0, airDrandMT(), 1, 0, 2*AIR_PI);
    cc = dist*edge*cos(theta);
    ss = dist*edge*sin(theta);
    ELL_3V_SCALE_ADD3(grad, 1.0, grad, cc, perp0, ss, perp1);
    ELL_3V_NORM(grad, grad, len);
    grad += 3;
  }

  return 0;
}

/*
** Eval from UV
*/
void
tenGlyphBqdEvalUv(double eval[3], const double uv[2]) {
  double xx, yy, zz, ll;

  yy = AIR_AFFINE(0, uv[0], 1, -1, 1);
  if (uv[0] + uv[1] > 1) {
    zz = AIR_AFFINE(0, uv[1], 1, -1, 1) - 1 + yy;
    xx = 1;
  } else {
    xx = AIR_AFFINE(0, uv[1], 1, -1, 1) + yy + 1;
    zz = -1;
  }
  ELL_3V_SET(eval, xx, yy, zz);
  ELL_3V_NORM(eval, eval, ll);
  return;
}

void
_echoRayIntxUV_Sphere(echoIntx *intx) {
  echoPos_t u, v;

  if (intx->norm[0] || intx->norm[1]) {
    u = atan2(intx->norm[1], intx->norm[0]);
    intx->u = AIR_AFFINE(-AIR_PI, u, AIR_PI, 0.0, 1.0);
    v = -asin(intx->norm[2]);
    intx->v = AIR_AFFINE(-AIR_PI/2, v, AIR_PI/2, 0.0, 1.0);
  }
  else {
    intx->u = 0;
    /* this is valid because if we're here, then intx->norm[2]
       is either 1.0 or -1.0 */
    intx->v = AIR_AFFINE(1.0, intx->norm[2], -1.0, 0.0, 1.0);
  }
}

/*
** hist[0]: hue vs sat
** hist[1]: hue vs val
*/
void
imageProc(Nrrd *nhproj[3], Nrrd *nhist[2], unsigned int sH,
          float *rgb, unsigned int size0, unsigned int sXY,
          unsigned int overSampleNum, float overSampleScale) {
  unsigned int xyi, hi, si, vi, oi;
  float rr, gg, bb, hh, ss, vv, *hist[2];
  double rndA, rndB;

  nrrdZeroSet(nhist[0]);
  nrrdZeroSet(nhist[1]);
  hist[0] = AIR_CAST(float *, nhist[0]->data);
  hist[1] = AIR_CAST(float *, nhist[1]->data);
  for (xyi=0; xyi<sXY; xyi++) {
    rr = AIR_CLAMP(0, rgb[0], 255);
    gg = AIR_CLAMP(0, rgb[1], 255);
    bb = AIR_CLAMP(0, rgb[2], 255);
    rr = AIR_AFFINE(-1, rr, 256, 0, 1);
    gg = AIR_AFFINE(-1, gg, 256, 0, 1);
    bb = AIR_AFFINE(-1, bb, 256, 0, 1);
    dyeRGBtoHSV(&hh, &ss, &vv, rr, gg, bb);
    si = airIndexClamp(0, ss, 1, sH);
    vi = airIndexClamp(0, vv, 1, sH);

#define UPDATE_HIST(rnd)                                                \
    hi = airIndexClamp(0, hh + overSampleScale*(1-ss)*(rnd), 1, sH);    \
    hist[0][hi + sH*si] += 1.0/overSampleNum;                           \
    hist[1][hi + sH*vi] += 1.0/overSampleNum

    if (overSampleNum % 2 == 1) {
      airNormalRand(&rndA, NULL);
      UPDATE_HIST(rndA);
      overSampleNum -= 1;
    }
    for (oi=0; oi<overSampleNum; oi+=2) {
      airNormalRand(&rndA, &rndB);
      UPDATE_HIST(rndA);
      UPDATE_HIST(rndB);
    }
    rgb += size0;
  }
  nrrdProject(nhproj[0], nhist[0], 1, nrrdMeasureHistoMean, nrrdTypeFloat);
  nrrdProject(nhproj[1], nhist[1], 1, nrrdMeasureHistoMean, nrrdTypeFloat);
  nrrdProject(nhproj[2], nhist[1], 1, nrrdMeasureSum, nrrdTypeFloat);
}

int
main(int argc, char *argv[]) {
  char *me;
  hestOpt *hopt=NULL;
  airArray *mop;
  
  pushEnergySpec *ensp;
  unsigned int pi, xi, nn;
  double xx, supp, del;

  mop = airMopNew();
  me = argv[0];
  hestOptAdd(&hopt, "energy", "spec", airTypeOther, 1, 1, &ensp, NULL,
             "specification of force function to use",
             NULL, NULL, pushHestEnergySpec);
  hestParseOrDie(hopt, argc-1, argv+1, NULL,
                 me, info, AIR_TRUE, AIR_TRUE, AIR_TRUE);
  airMopAdd(mop, hopt, (airMopper)hestOptFree, airMopAlways);
  airMopAdd(mop, hopt, (airMopper)hestParseFree, airMopAlways);

  fprintf(stderr, "%s: parsed energy \"%s\", with %u parms.\n", me,
          ensp->energy->name, ensp->energy->parmNum);
  for (pi=0; pi<ensp->energy->parmNum; pi++) {
    fprintf(stderr, "%u: %g\n", pi, ensp->parm[pi]);
  }
  fprintf(stderr, "\n");

  nn = 600;
  supp = ensp->energy->support(ensp->parm);
  del = AIR_DELTA(0, 2, nn, 0, supp);
  for (xi=1; xi<nn; xi++) {
    double x0, x1, ee, ff, e0, e1, dummy;
    xx = AIR_AFFINE(0, xi,   nn, 0, supp);
    x1 = AIR_AFFINE(0, xi+1, nn, 0, supp);
    x0 = AIR_AFFINE(0, xi-1, nn, 0, supp);
    ensp->energy->eval(&e1, &dummy, x1, ensp->parm);
    ensp->energy->eval(&e0, &dummy, x0, ensp->parm);
    ensp->energy->eval(&ee, &ff, xx, ensp->parm);
    printf("%g %g %g %g\n", xx, ee, ff, (e1 - e0)/del);
  }

  airMopOkay(mop);
  return 0;
}

void
_echoRayIntxUV_TriMesh(echoIntx *intx) {
  echoPos_t u, v, norm[3], len;
  echoTriMesh *trim;

  trim = TRIMESH(intx->obj);
  ELL_3V_SUB(norm, intx->pos, trim->meanvert);
  ELL_3V_NORM(norm, norm, len);
  if (norm[0] || norm[1]) {
    u = atan2(norm[1], norm[0]);
    intx->u = AIR_AFFINE(-AIR_PI, u, AIR_PI, 0.0, 1.0);
    v = -asin(norm[2]);
    intx->v = AIR_AFFINE(-AIR_PI/2, v, AIR_PI/2, 0.0, 1.0);
  }
  else {
    intx->u = 0;
    intx->v = AIR_AFFINE(1.0, norm[2], -1.0, 0.0, 1.0);
  }
}

int
main(int argc, char *argv[]) {
  char *me;
  unsigned int ii, NN;
  double min, max, *out;
  Nrrd *nout;

  me = argv[0];
  if (5 != argc) {
    usage(me);
  }
  if (3 != (sscanf(argv[2], "%u", &NN)
            + sscanf(argv[1], "%lf", &min)
            + sscanf(argv[3], "%lf", &max))) {
    fprintf(stderr, "%s: couldn't parse %s %s %s double uint double",
            me, argv[1], argv[2], argv[3]);
    usage(me);
  }
  nout = nrrdNew();
  if (nrrdAlloc_va(nout, nrrdTypeDouble, 2,
                   AIR_CAST(size_t, 5),
                   AIR_CAST(size_t, NN))) {
    fprintf(stderr, "%s: trouble allocating:\n%s", me,
            biffGetDone(NRRD));
    exit(1);
  }
  out = AIR_CAST(double *, nout->data);
  for (ii=0; ii<NN; ii++) {
    double xx, rr, ff, gg;
    xx = AIR_AFFINE(0, ii, NN-1, min, max);
    rr = exp(xx);
    ff = airFastExp(xx);
    gg = airExp(xx);
    if (rr < 0 || ff < 0 || gg < 0
        || !AIR_EXISTS(rr) || !AIR_EXISTS(ff) || !AIR_EXISTS(gg)) {
      fprintf(stderr, "%s: problem: %f -> real %f, approx %f %f\n",
              me, xx, rr, ff, gg);
      exit(1);
    }
    out[0] = rr;
    out[1] = ff;
    out[2] = (ff-rr)/rr;
    out[3] = gg;
    out[4] = (gg-rr)/rr;
    out += 5;
  }

  if (nrrdSave(argv[4], nout, NULL)) {
    fprintf(stderr, "%s: trouble saving:\n%s", me,
            biffGetDone(NRRD));
    exit(1);
  }
  exit(0);
}