void MakeCylModelStereo(PicWidth, FocalLength, CamSeparation, Cellsize,
			MaxRange, EvLo, EvHi, smodels)
MainInt PicWidth;
MainFloat FocalLength, CamSeparation, Cellsize, MaxRange;
MainInt  EvLo, EvHi;
CylSensorModelArray *smodels;
/*** 3D ***************************************************************\
*	MakeCylModelStereo(PicWidth, FocalLength,
*                 CamSeparation, GridResolution,ProbLo, ProbHi)
*	sensor model for stereoscopic ranging
*	PicWidth  image width, in pixels   FocalLength  in pixels
*       CamSeparation  stereo pair camera separation, in meters
*       Cellsize  cell to cell grid spacing, in meters
*       MaxRange is the length of the longest ray to paint
*       EvLo  lowest "empty" evidence (log odds) to use
*       EvHi  highest "occupied" evidence to use
*	smodels holds the resulting model
\**********************************************************************/
    {
     BestFloat lr, ld, hr, hd, r, d, dd, r0, rng, rngmax, rngmin, dr, t;
     BestFloat  rb;
     MainInt rsz, dsz, ir, id, k, gi, lrsz, nmods, displacement;
     Int8 *grid; Int16 *prof;

#define broad 8.0 /* beam broadening fudge */

     nmods=PicWidth/2;
     (*smodels).nmodels = nmods;
     (*smodels).models = (CylSensorModel *) calloc(nmods,sizeof(CylSensorModel));

     for (k=0; k<nmods; k++)
       { 
	   ld = lr = 0;
	   displacement = PicWidth/2-k;
	   rng = CamSeparation*FocalLength/displacement;
	   rngmax = CamSeparation*FocalLength/(MAX(displacement,2)-1);
	   rngmin = CamSeparation*FocalLength/(displacement+1);
	   if ((t=rngmax-rngmin) < Cellsize)
	     {
		 t = (Cellsize-t)/2;
		 rngmin -= t;  rngmax += t;
	     };

	   hd = MIN(rngmax,MaxRange);
	   hr = broad*hd/FocalLength;

	   dsz = (hd-ld)/Cellsize+0.5;
	   rsz = (hr-lr)/Cellsize;  rsz++;
	   lrsz = ILOG2C(rsz);  /* for r, lr is center, not edge, of row 0 */

	   (*smodels).models[k].rangem = rng;
	   (*smodels).models[k].rsize = rsz;  (*smodels).models[k].dsize = dsz;
	   (*smodels).models[k].lorv = lr;  (*smodels).models[k].lodv = ld;
	   (*smodels).models[k].hirv = hr;  (*smodels).models[k].hidv = hd;
	   (*smodels).models[k].modelm = grid =
	     (Int8 *) calloc(dsz<<lrsz,sizeof(Int8));   /* The grid */
	   (*smodels).models[k].profile = prof =
	     (Int16 *) calloc(dsz,sizeof(Int16));   /* radius silhouette */

	   dd = (hd-ld)/dsz;    d = 0.5*dd + ld;
	   dr = (hr-lr)/rsz;    r0 = lr - dr;  /* lr is midline of row 0 */
	   if (rsz==1 && lr==0.0)  /* thin ray! zoom zoom */
	     for (id=0; id<dsz; ++id)
	       {
		   prof[id]=0;
		   grid[id] = d<rngmin?EvLo:(d<=rngmax?EvHi:0);
		   d += dd;
	       }
	   else /* fat ray: more work */
	     for (id=0; id<dsz; ++id)
	       { 
		   rb = broad*d/FocalLength; /* beam radius at d */
		   r = r0;  gi = (id << lrsz) - 1;  prof[id]=0;
		   for (ir=0; ir<rsz; ++ir)
		     { r += dr;  gi++;
		       if (r>rb) grid[gi]=0; else
			 grid[id] = d<rngmin?EvLo:(d<=rngmax?EvHi:0);
		       if ((grid[gi])) prof[id]=ir+1; /* if nonzero! */
		   }
		   d += dd;
	       }
       }
   }