partibrains.c

    sp->p.x[0] = star.x * sl->scaledby;
    sp->p.x[1] = star.y * sl->scaledby;
    sp->p.x[2] = star.z * sl->scaledby;
    if(dfltvars) {
	sp->val[0] = exp((-18-star.magnitude)*.921/*log(100)/5*/);
	sp->val[1] = star.color;
	sp->val[2] = star.radius;
    } else {
	for(vp = &sp->val[0], cp = st->sdbvars; *cp; cp++, vp++) {
	    switch(*cp) {
	    case 'm': *vp = exp((-18-star.magnitude)*.921/*log(100)/5*/); break;
	    case 'M': *vp = star.magnitude; break;
	    case 'c': *vp = star.color; break;
	    case 'r': *vp = star.radius; break;
	    case 'o': *vp = star.opacity; break;
	    case 'g': *vp = star.group; break;
	    case 't': *vp = star.type; break;
	    case 'x': *vp = star.dx; break;
	    case 'y': *vp = star.dy; break;
	    case 'z': *vp = star.dz; break;
	    case 'S': *vp = sqrt(star.dx*star.dx + star.dy*star.dy + star.dz*star.dz); break;
	    case 'n': *vp = star.num; break;
	    default: *vp = 1; break;
	    }
	}
    }

    if(i == 0) {
	for(k = 0; k < nvars; k++)
	    sum[k] = min[k] = max[k] = sp->val[k];
    } else {
	for(k = 0; k < nvars; k++) {
	    if(min[k] > sp->val[k]) min[k] = sp->val[k];
	    else if(max[k] < sp->val[k]) max[k] = sp->val[k];
	    sum[k] += sp->val[k];
	}
    }
  }
  sl->nspecks = i;
  sl->sizedby = 0;
  sl->coloredby = 1;

  /* Update statistics */
  if(sl->nspecks > 0) {
    struct valdesc *vdp = &st->vdesc[st->curdata][0];
    for(i = 0; i < nvars; i++, vdp++) {
	if(vdp->min > min[i]) vdp->min = min[i];
	if(vdp->max < max[i]) vdp->max = max[i];
	vdp->nsamples += sl->nspecks;
	vdp->sum += sum[i];
	vdp->mean = vdp->sum / vdp->nsamples;

	if(vdp->name[0] == '\0') {
	    char *name = "unk";
	    switch(st->sdbvars[i]) {
	    case 'm': name = "lumsdb"; break;
	    case 'M': name = "magsdb"; break;
	    case 'c': name = vdp->max > 16384 ? "rgb565" : "colorsdb";
		      vdp->cexact = 1;
		      break;
	    case 'r': name = "radius"; break;
	    case 'o': name = "opacity"; break;
	    case 'g': name = "group"; break;
	    case 't': name = "type"; break;
	    case 'x': name = "dx"; break;
	    case 'y': name = "dy"; break;
	    case 'z': name = "dz"; break;
	    case 'S': name = "speed"; break;
	    case 'n': name = "number"; break;
	    }
	    strcpy(vdp->name, name);
	}
    }
    specks_recolor( st, sl, st->coloredby );
    specks_resize( st, sl, st->sizedby );
  }

  /* Add to running list */
  slp = specks_timespecksptr( st, st->curdata, timestep );
  sl->next = *slp;
  *slp = sl;

  fclose(inf);
    
}

void specks_timerange( struct stuff *st, double *tminp, double *tmaxp )
{
  if(!st->dyndata
#if USE_KIRA
	 || !get_parti_time_range( st, tminp, tmaxp )
#endif
   ) {
    *tminp = 0;
    *tmaxp = st->ntimes == 0 ? 0 : st->ntimes - 1;
  }
}

int specks_get_datastep( struct stuff *st )
{
  return st->curtime;
}

double specks_get_realtime( struct stuff *st )
{
  return st->currealtime;
}

void set_interest_point( Point *p )
{
#if CAVEMENU
  char cmd[80];
  sprintf(cmd, "interest %g %g %g", p->x[0], p->x[1], p->x[2]);
  VIDI_queue_commandstr( cmd );
#else

  parti_center( p );
#endif
}

/* Only specks_timespecksptr() extends the spans of time or datasets */
struct specklist **
specks_timespecksptr( struct stuff *st, int dataset, int timestep )
{
  int d, needroom;
  struct specklist **na, **nan;
  void **ndf;
  char **nfn;

  if((timestep >= st->ntimes || dataset >= st->ndata)) {
    needroom = st->timeroom;
    if(timestep >= st->timeroom)
	needroom = 2*timestep + 15;

    for(d = 0; d < st->ndata || (dataset < MAXFILES && d <= dataset); d++) {

	if(needroom == st->timeroom && d < st->ndata)
	    continue;

	na = NewN( struct specklist *, needroom );
	nan = NewN( struct specklist *, needroom );
	ndf = NewN( void *, needroom );
	nfn = NewN( char *, needroom );
	memset(na, 0, needroom * sizeof(*na));
	memset(nan, 0, needroom * sizeof(*nan));
	memset(ndf, 0, needroom * sizeof(*ndf));
	memset(nfn, 0, needroom * sizeof(*nfn));
	if(d < st->ndata && st->anima[d])
	    memcpy( na, st->anima[d], st->ntimes * sizeof(*na) );

	if(d < st->ndata && st->annot[d])
	    memcpy( nan, st->annot[d], st->ntimes * sizeof(*nan) );

	if(d < st->ndata && st->datafile[d])
	    memcpy( ndf, st->datafile[d], st->ntimes * sizeof(*ndf) );

	if(d < st->ndata && st->fname[d])
	    memcpy( nfn, st->fname[d], st->ntimes * sizeof(*nfn) );

	/* Don't free old pointers, just in case they're in use. */
	st->anima[d] = na;
	st->annot[d] = nan;
	st->datafile[d] = ndf;
	st->fname[d] = nfn;
    }
    st->timeroom = needroom;

    if(timestep >= st->ntimes)
	st->ntimes = timestep + 1;
    if(dataset >= st->ndata && dataset < MAXFILES)
	st->ndata = dataset + 1;
  }

  return (timestep >= 0 && timestep < st->ntimes &&
			dataset >= 0 && dataset < st->ndata)
	? &st->anima[dataset][timestep] : NULL;
}

struct specklist *
specks_timespecks( struct stuff *st, int dataset, int timestep )
{
  return (timestep >= 0 && timestep < st->ntimes &&
			dataset >= 0 && dataset < st->ndata)
	? st->anima[dataset][timestep] : NULL;
}

void specks_reupdate( struct stuff *st, struct specklist *sl )
{
  struct specklist *tsl;

  if(sl != NULL && sl->threshseq != st->threshseq) {
    for(tsl = sl; tsl != NULL; tsl = tsl->next)
	specks_rethresh( st, tsl, st->threshvar );
  }

  if(sl != NULL && sl->colorseq != st->colorseq) {
    for(tsl = sl; tsl != NULL; tsl = tsl->next)
	specks_recolor( st, tsl, st->coloredby );
  }

  if(sl != NULL && sl->sizeseq != st->sizeseq) {
    for(tsl = sl; tsl != NULL; tsl = tsl->next)
	specks_resize( st, tsl, st->sizedby );
  }
}

void specks_set_time( struct stuff *st, double newtime )
{
  static Point lastinterest;

  clock_set_time( st->clk, newtime );
  specks_set_timestep( st );
  if(st->sl == NULL && st->ntimes > 1 && st->clk->parent == NULL) {
    /* Skip blank time-slots -- keep incrementing until either:
     *  - we find a time-slot that has (or could have) some data, or
     *  - we've run through all time-steps (avoid infinite loops!).
     */
    int nudges, ts;
    for(nudges = 0; nudges < st->ntimes; nudges++) {
	ts = (st->curtime + nudges) % st->ntimes;
	if( specks_timespecks( st, st->curdata, ts ) != NULL ||
	    (st->datafile[st->curdata] != NULL &&
		st->datafile[st->curdata][ts] != NULL) )
	    break;
    }
    clock_set_time( st->clk, ts );
    specks_set_timestep( st );
  }
  /* st->playnext = now + (st->fspeed != 0 ? 1/st->fspeed : 0); */

#if CAVE
  if(ment.tknob[0] == NULL) {	/* if not SC99DEMO */
    struct specklist *sl = st->sl;
    if(sl != NULL && memcmp(&sl->interest, &lastinterest, sizeof(Point))
	    && (sl->interest.x[0]!=0 || sl->interest.x[1]!=0
					    || sl->interest.x[2]!=0)) {
	char str[128];
	float scale = .004;
	lastinterest = sl->interest;
	set_interest_point( &lastinterest );
	sprintf(str, "\002setjump peak %g %g %g 0 0 0 %g",
	    lastinterest.x[0], lastinterest.x[1]-scale*5, lastinterest.x[2] + scale*5,
	    scale);
	VIDI_queue_commandstr( str );
    }
  }
#endif

  specks_reupdate( st, st->sl );

#if CAVE
  specks_refresh_menu( st );
#endif
}

static int specks_freenow( struct specklist **slp, int maxage )
{
  struct specklist *sl, **sprev;
  int any = 0;

  for(sprev = slp; (sl = *sprev) != NULL && sl->used <= maxage; ) {
    if(sl->used <= maxage) {
	*sprev = sl->freelink;
	if(sl->specks != NULL)
	    Free(sl->specks);
	Free(sl);
	any++;
    } else {
	/* too recent - might still be in use */
	sprev = &sl->freelink;
    }
  }
  return any;
}

void specks_discard( struct stuff *st, struct specklist **slp )
{
  struct specklist *sl, *slnext;

  for(sl = *slp; sl != NULL; sl = slnext) {
    slnext = sl->next;
    sl->freelink = st->scrap;
    st->scrap = sl;
  }
  *slp = NULL;
}

int specks_purge( void *vst, int nbytes, void *aarena )
{
  struct stuff *st = (struct stuff *)vst;
  int oldused = st->used;
  int oldtime = -1, oldds = -1;
  int t, ds;
  struct specklist *sl;

#ifdef sgi
  static int first = 1;
  struct mallinfo mi;
  mi = amallinfo( aarena );

  if(first) {
    first = 0;
    msg("Purging %dKbyte shmem arena (currently %dK used in %d blks, %dK free)\n",
	mi.arena>>10, mi.uordblks>>10, mi.ordblks, mi.fordblks>>10);
  }
#endif

  /* Free any known scrap first */

#define OLD_ENOUGH  4

  if(specks_freenow( &st->scrap, st->used - OLD_ENOUGH ) > 0)
    return 1;

  for(t = 0; t < st->ntimes; t++) {
    if(t == st->curtime) continue;
    for(ds = 0; ds < st->ndata; ds++) {
	sl = st->anima[ds][t];
	if(sl != NULL && sl != st->sl && sl->used < oldused) {
	    oldused = st->used;
	    oldtime = t;
	    oldds = ds;
	}
    }
  }
  if(oldtime >= 0) {
    specks_discard( st, &st->anima[oldds][oldtime] );
    specks_freenow( &st->scrap, oldtime );
    return 1;	/* We freed something, so try allocating again */
  } else {
    msg("Ran out of shmem, couldn't find anything more to purge\n");
#ifdef sgi
    msg("%dKbyte shmem arena (currently %dK used in %d blks, %dK free)\n",
	mi.arena>>10, mi.uordblks>>10, mi.ordblks, mi.fordblks>>10);

#endif
    return 0;	/* No progress made -- give up */
  }
}


/*
 * Just stash these values in our frame function so they won't change visibly
 * during a frame.  Each frame function will do this; we'll just hope that
 * they don't change as the various cave-wall processes start.
 */
void specks_current_frame( struct stuff *st, struct specklist *sl )
{
  st->frame_sl = sl;
  st->frame_time = st->curtime;
  st->frame_data = st->curdata;
  st->frame_annotation = st->annotation;
  specks_reupdate( st, sl );
}

#define	MAXXYFAN 16

extern void specks_draw_boxes( struct stuff *st, struct AMRbox *boxes, int levelmask, Matrix Ttext, int oriented );

struct cpoint {
    int rgba;
    Point p;
};

void dumpcpoints( struct cpoint *cp, int n )
{
    while(--n >= 0) {
	glColor4ubv( (GLubyte *)&cp->rgba );
	glVertex3fv( &cp->p.x[0] );
	cp++;
    }
}

static Point depth_fwd;
static float depth_d;

struct order {
  float z;
  struct speck *sp;
  struct specklist *sl;
};

static int depthcmp( const void *a, const void *b )
{
  return ((struct order *)a)->z < ((struct order *)b)->z ?    1
	: ((struct order *)a)->z > ((struct order *)b)->z ? -1 : 0;
}

static int additive_blend;

void sortedpolys( struct stuff *st, struct specklist *slhead, Matrix *Tc2wp, float radperpix, float polysize )
{
  struct speck *sp, *sbase;
  struct order *op, *obase;
  int i, k, total, skip;
  struct specklist *sl;
  int usethresh = st->usethresh&P_USETHRESH ? THRESHBIT : 0;
  int bps = 0;
  int prevrgba = -1;
  int usearea = st->polyarea;
  int sizevar = st->polysizevar;
  int polyorivar = st->polyorivar0;
  int texturevar = st->texturevar;
  int txno;
  int texturing = -1;
  float s;
  float polyminrad = st->polymin * radperpix;
  float polymaxrad = st->polymax * radperpix;
  float mins2d = polyminrad * polyminrad;
  int rgba;
  int alpha = st->alpha * 255;
  int nfan = st->npolygon<MAXXYFAN ? st->npolygon : MAXXYFAN;
  float xyfan[MAXXYFAN][2];
  Point sfan[MAXXYFAN], pfan[MAXXYFAN];
  Matrix Tc2w = *Tc2wp;
  float scl = vlength( (Point *)&Tc2w.m[0*4+0] );
  float fanscale;
  Texture *wanttx;
  int additive = additive_blend;
  int wantblend = additive;

  int useclip = (st->clipbox.level > 0);
  Point clipp0 = st->clipbox.p0;
  Point clipp1 = st->clipbox.p1;

  for(total = 0, sl = slhead; sl != NULL; sl = sl->next) {
    if(sl->text != NULL || sl->nspecks == 0 || sl->special != SPECKS)
	continue;
    skip = st->subsample;
    if(sl->subsampled != 0)	/* if already subsampled */
	skip /= sl->subsampled;
    if(skip <= 0) skip = 1;

    if(bps < sl->bytesperspeck) bps = sl->bytesperspeck;
    if(usethresh) {
	for(i=sl->nspecks, sp=sl->specks; i>0; i-=skip, sp=NextSpeck(sp,sl,skip)) {
	    if((usethresh & sp->rgba) == 0)
		total++;
	}
    } else {
	total += sl->nspecks / skip;
    }
  }

  obase = op = (struct order *)malloc( (total+1) * sizeof(struct order) );
  for(sl = slhead; sl != NULL; sl = sl->next) {
    if(sl->text != NULL || sl->nspecks == 0 || sl->special != SPECKS)
	continue;
    skip = st->subsample;
    if(sl->subsampled != 0)	/* if already subsampled */
	skip /= sl->subsampled;
    if(skip <= 0) skip = 1;
    for(i=sl->nspecks, sp=sl->specks; i > 0; i-=skip, sp=NextSpeck(sp,sl,skip)) {
	float dist;
	if(usethresh & sp->rgba)
	    continue;
	dist = VDOT( &sp->p, &depth_fwd ) + depth_d;
	if(dist < 0)
	    continue;
	if(useclip &&
	  (sp->p.x[0] < clipp0.x[0] ||
	   sp->p.x[0] > clipp1.x[0] ||
	   sp->p.x[1] < clipp0.x[1] ||
	   sp->p.x[1] > clipp1.x[1] ||
	   sp->p.x[2] < clipp0.x[2] ||
	   sp->p.x[2] > clipp1.x[2]))
	    continue;
	op->z = dist;
	op->sp = sp;
	op->sl = sl;
	op++;
    }
  }

  total = op - obase;
  qsort( obase, total, sizeof(*obase), depthcmp );

  prevrgba = 0;

  /* Build prototype fan -- unit disk in screen plane */
  /* Scale the fan big enough that the unit disk is inscribed in our polygon */
  fanscale = 1 / (scl * cos(M_PI/nfan));
  for(i = 0; i < nfan; i++) {
    float theta = 2*M_PI*(i+0.5f)/nfan;
    xyfan[i][0] = cos(theta);
    xyfan[i][1] = sin(theta);
    vcomb( &sfan[i], xyfan[i][0] * fanscale, (Point *)&Tc2w.m[0*4+0],
		     xyfan[i][1] * fanscale, (Point *)&Tc2w.m[1*4+0] );
  }

  if(st->usetextures == 0 || SMALLSPECKSIZE(texturevar) > bps)
    texturevar = -1;
  if(SMALLSPECKSIZE(polyorivar) > bps)
    polyorivar = -1;

  for(i = 0, op = obase; i < total; i++, op++) {
    float dist, size;

    sp = op->sp;
    dist = op->z;
    size = sp->val[sizevar] * polysize;
    if(usearea) {
	if(size < dist * dist * mins2d)
	    continue;
	size = sqrtf(size);
    } else {
	if(size < dist * polyminrad)
	    continue;
    } 
    if(size > dist * polymaxrad)
	size = dist * polymaxrad;

    rgba = sp->rgba & ~THRESHBIT;
    if(rgba != prevrgba) {
	prevrgba = rgba;
	rgba = RGBALPHA( prevrgba, alpha );
	glColor4ubv( (GLubyte *)&rgba );
    }

    if(texturevar >= 0 &&
	    (txno = sp->val[texturevar]) >= 0 &&
	    txno < st->ntextures &&
	    (wanttx = st->textures[txno]) != NULL) {

	txbind( wanttx, &texturing );
	wantblend = (wanttx->flags & TXF_ADD) ? 1 : additive_blend;
    } else if(texturing) {
	glDisable( GL_TEXTURE_2D );
	texturing = 0;
    }

    if(wantblend != additive) {
	additive = wantblend;
	glBlendFunc( GL_SRC_ALPHA, additive ? GL_ONE : GL_ONE_MINUS_SRC_ALPHA );
    }

    if(polyorivar >= 0 && sp->val[polyorivar] < 9) {
	float *xv = &sp->val[polyorivar];
	float *yv = &sp->val[polyorivar+3];

	glBegin( GL_TRIANGLE_FAN );
	if(texturing) {
	    for(k = 0; k < nfan; k++) {
		glTexCoord2fv( &xyfan[k][0] );
		glVertex3f(
		    sp->p.x[0] + size*(xyfan[k][0]*xv[0] + xyfan[k][1]*yv[0]),
		    sp->p.x[1] + size*(xyfan[k][0]*xv[1] + xyfan[k][1]*yv[1]),
		    sp->p.x[2] + size*(xyfan[k][0]*xv[2] + xyfan[k][1]*yv[2]));
	    }
	} else {
	    for(k = 0; k < nfan; k++) {
		glVertex3f(
		    sp->p.x[0] + size*(xyfan[k][0]*xv[0] + xyfan[k][1]*yv[0]),
		    sp->p.x[1] + size*(xyfan[k][0]*xv[1] + xyfan[k][1]*yv[1]),
		    sp->p.x[2] + size*(xyfan[k][0]*xv[2] + xyfan[k][1]*yv[2]));
	    }
	}
	glEnd();

    } else {
	glBegin( GL_TRIANGLE_FAN );
	if(texturing) {
	    for(k = 0; k < nfan; k++) {
		glTexCoord2fv( &xyfan[k][0] );
		glVertex3f(
		    sp->p.x[0] + size*sfan[k].x[0],
		    sp->p.x[1] + size*sfan[k].x[1],
		    sp->p.x[2] + size*sfan[k].x[2] );
	    }
	} else {
	    for(k = 0; k < nfan; k++) {
		glVertex3f(
		    sp->p.x[0] + size*sfan[k].x[0],
		    sp->p.x[1] + size*sfan[k].x[1],
		    sp->p.x[2] + size*sfan[k].x[2] );
	    }
	}
	glEnd();
    }
  }
  free(obase);
  if(texturing > 0)
    txbind( NULL, NULL );
}
  
  

void drawspecks( struct stuff *st )
{
  int i, slno, k;
  int rgba, alpha, prevrgba = 0;
  float prevsize = 0;
  struct specklist *sl, *slhead;
  register struct speck *p;
  Matrix Tw2c, Tc2w, Ttext, Tproj, Ttemp;
  static Point zero = {0,0,0};
  int xywh[4];
  float radperpix;
  Point tp, fan[MAXXYFAN];
  Point eyepoint;
  Point fwd;
  float fwdd;
  float tscale, scl, fanscale;
  int skip;
  static int nxyfan = 0;
  static float xyfan[MAXXYFAN][2];
  static unsigned char randskip[256];
  int randix = 0;
  int fast = st->fast;
  int inpick = st->inpick;
  int usethresh = st->usethresh&P_USETHRESH ? THRESHBIT : 0;
  int fixeddist;
  float polyminrad, polymaxrad;
  float threshmin = st->thresh[0];
  float threshmax = st->thresh[1];
  int useclip = (st->clipbox.level > 0);
  Point clipp0 = st->clipbox.p0;
  Point clipp1 = st->clipbox.p1;

  float plum = st->psize;

  float knee1dist2 = st->fadeknee1 * st->fadeknee1;
  float knee2dist2 = st->fadeknee2 * st->fadeknee2;
  float orthodist2 = st->fadeknee2 * st->fadeknee2;
  float steep2knee2 = st->knee2steep * st->knee2steep / knee2dist2;
  float faderball2 = 1 / (st->fadeknee2 * st->fadeknee2);
  Point fadecen = st->fadecen;
  enum FadeType fademodel = st->fade;

  if(!st->useme)
    return;

  switch(fademodel) {
  case F_CONSTANT:
	if(orthodist2 <= 0)
	    orthodist2 = 1;
	break;
  case F_KNEE12:
	if(st->fadeknee1 >= st->fadeknee2 || st->fadeknee1 <= 0)
	    fademodel = F_KNEE2;	/* and fall into... */
  case F_KNEE2:
	if(st->fadeknee2 <= 0)
	    fademodel = F_SPHERICAL;
	break;
  case F_LREGION:
	if(st->fadeknee2 <= 0) faderball2 = 1;
	break;
  }


  { float r=0,g=0,b=0;	/* Ugh. Allow background to be non-black */
    sscanf(parti_bgcolor(NULL), "%f%f%f", &r,&g,&b);
    additive_blend = (r+g+b == 0);
  }

  if(st->clipbox.level != 0) {
    GLdouble plane[4];

    if(st->clipbox.level == 1) {
	struct AMRbox b[2];
	b[0] = st->clipbox;
	b[0].level = 0;
	b[1].level = -1;
	specks_draw_boxes( st, b, ~0, Tidentity, 1 );
    }
    for(i = 0; i < 3; i++) {
	plane[0] = plane[1] = plane[2] = 0;
	plane[i] = 1;
	plane[3] = -st->clipbox.p0.x[i];
	glClipPlane( GL_CLIP_PLANE0 + i, plane );
	glEnable( GL_CLIP_PLANE0 + i );
	plane[i] = -1;
	plane[3] = st->clipbox.p1.x[i];
	glClipPlane( GL_CLIP_PLANE0 + 3 + i, plane );
	glEnable( GL_CLIP_PLANE0 + 3 + i );
    }
  }


  if(nxyfan != st->npolygon && st->npolygon > 0) {
    if(st->npolygon > MAXXYFAN) st->npolygon = MAXXYFAN;
    nxyfan = st->npolygon;
    fanscale = 1 / cos(M_PI/nxyfan);
    for(i = 0; i < nxyfan; i++) {
	float th = (i+.5f)*2*M_PI / nxyfan;
	xyfan[i][0] = cos(th) * fanscale;
	xyfan[i][1] = sin(th) * fanscale;
    }
    srandom(11);
    for(i = 0; i < 256; i++)
	randskip[i] = random() & 0xFF;
  }


  alpha = st->alpha * 255;
  rgba = RGBALPHA( RGBWHITE, alpha );	/* BIG-ENDIAN (1,1,1,alpha) */

  /* Find displacements which lie in the screen plane, for making
   * billboard-style polygonal patches, and for text.
   */
  glGetFloatv( GL_MODELVIEW_MATRIX, Tw2c.m );
  eucinv( &Tc2w, &Tw2c );
  scl = vlength( (Point *)&Tw2c.m[0] );
  for(i = 0; i < nxyfan; i++) {
    tp.x[0] = scl*st->polysize*xyfan[i][0];
    tp.x[1] = scl*st->polysize*xyfan[i][1];
    tp.x[2] = 0;
    vtfmvector( &fan[i], &tp, &Tc2w );
  }

  /* Find projection matrix and screen (well, viewport) size,
   * so we can convert angular sizes to screen (pixel) sizes,
   * in radians per pixel.
   */
  glGetFloatv( GL_PROJECTION_MATRIX, Tproj.m );
  glGetIntegerv( GL_VIEWPORT, xywh );
  radperpix = 1 / (.5*xywh[2] * Tproj.m[0*4+0]);

  /* Construct a "forward" vector in object coords too, for measuring
   * distance from camera plane.  Note camera looks toward its -Z axis not +Z!
   */
  tp.x[0] = 0, tp.x[1] = 0, tp.x[2] = -1;
  vtfmvector( &fwd, &tp, &Tc2w );
  vunit( &fwd, &fwd );
  /*
   * Actually we want a plane equation, whose value is zero in the
   * eye plane.  Camera-space distance from camera plane = 
   *		vdot( &objectpoint, &fwd ) + fwdd.
   */
  vtfmpoint( &eyepoint, &zero, &Tc2w );
  fwdd = -vdot( &eyepoint, &fwd );

  tscale = scl * st->textsize;
  mcopy( &Ttemp, &Tidentity );
  Ttemp.m[0*4+0] = Ttemp.m[1*4+1] = Ttemp.m[2*4+2] = tscale;
  mmmul( &Ttext, &Ttemp, &Tc2w );
  vsettranslation( &Ttext, &zero );

  glDisable( GL_LIGHTING );

  /* Draw any boxes (even if we have no specks) for this timestep */
  if(st->useboxes && st->boxlevelmask != 0
	&& st->frame_time >= 0 && st->frame_time < st->boxtimes
	&& st->boxes[st->frame_time] != NULL) {
    specks_draw_boxes( st, st->boxes[st->frame_time], st->boxlevelmask, Ttext, 0 );
  }
  if(st->useboxes && st->staticboxes != NULL) {
    specks_draw_boxes( st, st->staticboxes, st->boxlevelmask, Ttext, st->boxaxes );
  }

  slhead = st->frame_sl;	/* st->sl as snapped by specks_ffn() */
  if(slhead == NULL)
    slhead = st->sl;		/* maybe there is no specks_ffn() */
  if(slhead == NULL)
    return;

#if USE_KIRA
  kira_draw( st, slhead, &Tc2w, radperpix );
#endif

  skip = st->subsample;
  if(slhead->subsampled != 0)	/* if already subsampled */
    skip /= slhead->subsampled;
  if(skip == 0) skip = 1;

  for(sl = slhead; sl != NULL; sl = sl->next)
    sl->used = st->used;

  if((unsigned int)st->sizedby <= MAXVAL
			&& (unsigned int)st->curdata < MAXFILES
			&& st->vdesc[st->curdata][st->sizedby].lum != 0) {
	plum *= st->vdesc[st->curdata][st->sizedby].lum;
  }
  if(st->subsample > 0 && st->everycomp)
      plum *= st->subsample;	/* Compensate for "every" subsampling */


  if(st->alpha >= 1) {
    glDisable(GL_BLEND);
    glEnable(GL_DEPTH_TEST);
  } else {
    glEnable(GL_BLEND);
    glBlendFunc( GL_SRC_ALPHA, additive_blend ? GL_ONE : GL_ONE_MINUS_SRC_ALPHA );
    glEnable(GL_DEPTH_TEST);
    glDepthMask( GL_FALSE );
  }

  if(inpick) glLoadName(0);
  

  if(st->usepoly && st->polysize > 0) {
    int texturing = 0;
    int texturevar = st->usetextures && st->texturevar >= 0
			&& st->texturevar < MAXVAL
		   ? st->texturevar : -1;
    int usearea = st->polyarea;
    int sizevar = st->polysizevar;
    float polysize = st->polysize;
    float mins2d;

    int txno;

    glMatrixMode( GL_TEXTURE );
    glLoadIdentity();
    glTranslatef( .5, .5, 0 );
    glScalef( st->txscale, st->txscale, st->txscale );
    glMatrixMode( GL_MODELVIEW );


    polyminrad = st->polymin * radperpix;
    polymaxrad = st->polymax * radperpix;
    mins2d = polyminrad*polyminrad;

    if(sizevar == -1) {
	/* If polygon size is tied to point size,
	 * then include pointsize scale factors in polygon scaling.
	 */
      if((unsigned int)st->sizedby <= MAXVAL
		&& (unsigned int)st->curdata < MAXFILES
		&& st->vdesc[st->curdata][st->sizedby].lum != 0)
	polysize *= st->vdesc[st->curdata][st->sizedby].lum;
      if(st->subsample > 0 && st->everycomp)
        polysize *= st->subsample; /* Compensate for "every" subsampling */
    }

    if(st->depthsort && !inpick) {
	depth_fwd = fwd;
	depth_d = fwdd;
	sortedpolys( st, slhead, &Tc2w, radperpix, polysize );

    } else {
      for(sl = slhead, slno = 1; sl != NULL; sl = sl->next, slno++) {
	if(sl->text != NULL || sl->special != SPECKS) continue;
	if(inpick) {
	    glLoadName(slno);
	    glPushName(0);
	}
	for(i = 0, p = sl->specks; i < sl->nspecks; i+=skip, p = NextSpeck( p, sl, skip )) {
	    float dist = VDOT( &p->p, &fwd ) + fwdd;
	    float size;

	    if(dist <= 0) continue;

	    if(usethresh & p->rgba)
		continue;

	    size = p->val[sizevar] * polysize;
	    if(usearea) {
		if(size < dist * dist * mins2d)
		    continue;
		size = sqrtf(size);
	    } else {
		if(size < dist * polyminrad)
		    continue;
	    } 
	    if(size > dist * polymaxrad)
		size = dist * polymaxrad;

	    rgba = p->rgba & ~THRESHBIT;
	    if(rgba != prevrgba) {
		prevrgba = rgba;
		rgba = RGBALPHA( prevrgba, alpha );
		glColor4ubv( (GLubyte *)&rgba );
	    }
	    if(st->polyorivar0 >= 0 && p->val[st->polyorivar0] < 9) {
		for(k = 0; k < nxyfan; k++) {
		    vcomb( &fan[k],
			size*xyfan[k][0], (Point *)&p->val[st->polyorivar0],
			size*xyfan[k][1], (Point *)&p->val[st->polyorivar0+3] );
		}
	    } else if(p->size != prevsize) {
		float s = scl*size;
		for(k = 0; k < nxyfan; k++) {
		    vcomb( &fan[k], s*xyfan[k][0], (Point *)&Tc2w.m[0*4+0],
				    s*xyfan[k][1], (Point *)&Tc2w.m[1*4+0] );
		}
		prevsize = size;
	    }

#define PFAN(vno, comp)  p->p.x[comp] + fan[vno].x[comp]

	    if(inpick) {
		glLoadName( i );
		glBegin( GL_TRIANGLE_FAN );
		for(k = 0; k < nxyfan; k++) {
		    glVertex3f( PFAN(k,0), PFAN(k,1), PFAN(k,2) );
		}
		glEnd();

	    } else if(texturevar >= 0
		    && (txno = p->val[texturevar]) >= 0
		    && txno < st->ntextures &&
		    st->textures[txno] != NULL) {

		txbind( st->textures[txno], &texturing );

		glBegin( GL_TRIANGLE_FAN );
		for(k = 0; k < nxyfan; k++) {
		    glTexCoord2fv( &xyfan[k][0] );
		    glVertex3f( PFAN(k,0), PFAN(k,1), PFAN(k,2) );
		}
		glEnd();

	    } else {
		if(texturing) {
		    texturing = 0;
		    glDisable( GL_TEXTURE_2D );
		}
		glBegin(GL_TRIANGLE_FAN);
		for(k = 0; k < nxyfan; k++)
		    glVertex3f( PFAN(k,0), PFAN(k,1), PFAN(k,2) );
		glEnd();
	    }
#undef PFAN

	}
	if(inpick) glPopName();
      }
    }
    if(texturing) {
	txbind( NULL, NULL );
	glDisable( GL_TEXTURE_2D );
    }
    glMatrixMode( GL_TEXTURE );
    glLoadIdentity();
    glMatrixMode( GL_MODELVIEW );
  }

  if(st->usepoint && !(st->useboxes == 2)) {

#define MAXPTSIZE 16	/* in half-point units */
#define PERBUCKET 64	/* max points per bucket */

    struct cpoint sized[MAXPTSIZE*2][PERBUCKET];
    int nsized[MAXPTSIZE*2];
    unsigned char invgamma[256];
    float invgam = (st->gamma <= 0) ? 0 : 1/st->gamma;

    for(i = 0; i < 256; i++)
	invgamma[i] = (int) (255.99 * pow( i/255., invgam ));

    if(inpick) {
	for(sl = slhead, slno = 1; sl != NULL; sl = sl->next, slno++) {
	    if(sl->text != NULL || sl->special != SPECKS) continue;
	    glLoadName(slno);
	    glPushName(0);
	    for(i = 0, p = sl->specks; i < sl->nspecks; i+=skip, p = NextSpeck(p, sl, skip)) {
		if(usethresh & p->rgba)
		    continue;

		glLoadName(i);
		glBegin(GL_POINTS);
		glVertex3fv( &p->p.x[0] );
		glEnd();
	    }
	    glPopName();
	}

    } else if(fast) {
	static unsigned char apxsize[MAXPTSIZE*MAXPTSIZE];
	unsigned char faintrand[256];
	int pxsize, oldpxsize;
	int pxmin, pxmax;

	pxmin = 256 * st->pfaint;
	if(st->plarge > MAXPTSIZE) st->plarge = MAXPTSIZE;
	pxmax = 256 * st->plarge * st->plarge;

	if(apxsize[1] == 0) {
	    for(i=0; i<COUNT(apxsize); i++)
		apxsize[i] = (int)ceil(sqrtf(i+1));
	}
	for(i = 0; i < 256; i++)
	    faintrand[i] = randskip[i] * st->pfaint;

	/* Render using fast (non-antialiased) points */
	glDisable( GL_POINT_SMOOTH );
	glEnable( GL_BLEND );
	glBlendFunc( GL_SRC_ALPHA, additive_blend ? GL_ONE : GL_ONE_MINUS_SRC_ALPHA );

	prevrgba = 0;
	prevsize = 0;

	sl = slhead;

	pxsize = oldpxsize = 1;
	glPointSize( pxsize );
	glColor4ubv( (GLubyte *)&rgba );

	for(i = 0; i < MAXPTSIZE*2; i++)
	    nsized[i] = 0;

	glBegin( GL_POINTS );
	for(sl = slhead; sl != NULL; sl = sl->next) {
	    if(sl->text != NULL || sl->special != SPECKS) continue;
	    for(i = 0, p = sl->specks; i < sl->nspecks; i+=skip, p= NextSpeck(p, sl, skip)) {
		int lum, myalpha;
		float dist = VDOT( &p->p, &fwd ) + fwdd;
		if(dist <= 0)	/* Behind eye plane */
		    continue;