partibrains.c

  }

  if(by < 0 || by >= MAXVAL)
    by = 0;

  vd = &st->vdesc[curdata][by];

  lmin = vd->lmin, lmax = vd->lmax;
  if(lmin == lmax && lmin == 0 || vd->lall) {
	vd->lmin = lmin = vd->min;
	vd->lmax = lmax = vd->max;
  }

  if(lmax == lmin)
    normal = 1;
  else
    normal = 1 / (lmax - lmin);

  if(st->useemph && st->emphsel.use != SEL_NONE) {
    for(i = 0; i < nel; i++, sp = NextSpeck(sp, sl, 1)) {
	sp->size = (sp->val[by] - lmin) * normal;
	if(SELECTED(sl->sel[i], &st->emphsel)) {
	    sp->size *= st->emphfactor;
	}
    }
  } else {
    /* normal case */
    for(i = 0; i < nel; i++, sp = NextSpeck(sp, sl, 1))
	sp->size = (sp->val[by] - lmin) * normal;
  }
}

void specks_datawait(struct stuff *st) {
#if USE_IEEEIO
    while(st->fetching && st->fetchpid > 0
	    && st->fetchtime == st->curtime
	    && st->fetchdata == st->curdata)
	usleep(50000);
#endif
}

#ifdef USE_IEEEIO

void specks_ieee_server( void *vst ) {
  struct stuff *st = (struct stuff *)vst;
  struct specklist *sl;

  prctl(PR_TERMCHILD);	/* Die when parent dies */
  /* Await a request */
  for(;;) {
    time_t then;
    while(st->fetching <= 0)
	sginap(5);
    then = time(NULL);
    sl = specks_ieee_read_timestep( st, st->subsample,
				st->fetchdata, st->fetchtime );
    if(then + 5 < time(NULL))
	msg("... got %x (%d particles) from %s", sl, specks_count(sl),
		st->fname[st->fetchdata][st->fetchtime]);
    if(st->curtime == st->fetchtime && st->curdata == st->fetchdata
			&& sl != NULL) {
	st->sl = sl;
	parti_redraw();
    }
    st->fetching = 0;
  }
}

#endif /*USE_IEEEIO*/

void specks_set_timebase( struct stuff *st, double timebase ) {
  clock_set_timebase( st->clk, timebase );
  /* parti_set_timebase( st, timebase ); */
}

void specks_set_speed( struct stuff *st, double newspeed ) {
  clock_set_speed( st->clk, newspeed );
  parti_set_speed( st, newspeed );
}

void specks_set_fspeed( struct stuff *st, double newspeed ) {
  /* ignore it! */
  /* or maybe use fspeed to specify finite clock resolution */
}

void specks_set_timestep( struct stuff *st )
{
  struct specklist *sl = st->sl;
  struct specklist **slp;
  double realtime, tmin, tmax;
  int timestep;

  st->used++;

  tmin = st->clk->tmin;
  tmax = st->clk->tmax;
  specks_timerange( st, &tmin, &tmax );
  if(st->usertrange) {
    /*if(tmin < st->utmin)*/ tmin = st->utmin;
    /*if(tmax > st->utmax)*/ tmax = st->utmax;
  }
  /* Only expand clock's existing time range */
  if(tmin > st->clk->tmin) tmin = st->clk->tmin;
  if(tmax < st->clk->tmax) tmax = st->clk->tmax;
  clock_set_range( st->clk, tmin, tmax, st->utwrap );
  realtime = clock_time( st->clk );
  timestep = realtime;

  if(st->dyn.enabled && st->dyn.getspecks) {
    struct specklist *sl;

    if(realtime == st->currealtime && st->sl != NULL && st->sl->used >= 0)
	return;
    sl = (*st->dyn.getspecks)(&st->dyn, st, realtime);
    if(sl == NULL)
	return;
    sl->used = 1;
    st->currealtime = realtime;
    st->curtime = timestep;

    slp = specks_timespecksptr( st, 0, 0 );
    st->sl = sl;
    /* Do we need this? */ /* *slp = sl; */ /* Guess not... */
    parti_set_timestep( st, realtime );
    return;
  }

  if(timestep >= st->ntimes) timestep = st->ntimes - 1;
  if(timestep < 0) timestep = 0;

  slp = specks_timespecksptr( st, st->curdata, timestep );
  if(*slp != NULL && (*slp)->subsampled > st->subsample) {
    specks_discard( st, slp );
  }

  sl = *slp;

#ifdef USE_IEEEIO

  if(sl == NULL && (st->fetching == 0 || st->datasync)
		&& st->curdata>=0 && st->curdata<st->ndata
		&& st->datafile[st->curdata] != NULL
		&& timestep>=0 && timestep<st->ntimes
		&& st->datafile[st->curdata][timestep] != NULL) {

    st->fetchdata = st->curdata;
    st->fetchtime = timestep;

    if(getenv("NO_SPROC") || st->datasync) {
	sl = specks_ieee_read_timestep( st, st->subsample,
                                st->fetchdata, st->fetchtime );
        msg("Got %x (%d particles) from %s  (d%d t%d)\n",
		sl, specks_count(sl), st->fname[st->fetchdata][st->fetchtime],
		st->fetchdata, st->fetchtime);

    } else {
	st->fetching = 1;

	if(st->fetchpid <= 0) {
	    st->fetchpid = sproc( specks_ieee_server, PR_SADDR|PR_SFDS, st );
	    if(st->fetchpid < 0)
		perror("sproc");
	}
    }
  }

  if(st->fetchpid > 0) {
    if(kill(st->fetchpid, 0) < 0) {
	perror("specks server vanished: kill -0");
	st->fetching = 0;
	st->fetchpid = 0;
    }
  }
#endif /*USE_IEEEIO*/

  if(st->skipblanktimes
	&& sl == NULL
	&& !(st->boxtimes>timestep && st->boxes[timestep]!=NULL)
	&& st->meshes[st->curdata][timestep] == NULL)
    return;

  st->sl = sl;   /* st->sl <= anima[][] */
  st->curtime = timestep;
  st->currealtime = timestep;
#if !CAVEMENU
  st->frame_time = st->curtime;	/* if non-CAVE, we have no frame-function */
#endif

  { struct specklist *asl = CURDATATIME(annot);
    specks_set_current_annotation( st, asl ? asl->text : NULL );
  }
  parti_set_timestep( st, timestep );

}

void specks_set_current_annotation( struct stuff *st, char *annotation )
{
  st->annotation = annotation;
#ifdef CAVEMENU
  if(ment.annot)
    menu_settitle( ment.annot, st->annotation ? st->annotation : "" );
#endif
}

void specks_add_annotation( struct stuff *st, char *annotation, int timestep )
{
  struct specklist *sl, **slp;
  int curtime = (timestep < 0) ? st->curtime : timestep;

  if(annotation == NULL) annotation = "";

  specks_timespecksptr( st, st->curdata, curtime );
  if(curtime >= st->ntimes) curtime = 0;
  slp = &st->annot[st->curdata][curtime];	/* safe -- no need for CURDATATIME() check */
  if((sl = *slp) == NULL) {
    sl = *slp = NewN( struct specklist, 1 );
    memset(*slp, 0, sizeof(**slp));
    sl->speckseq = ++st->speckseq;
  } else if(sl->text) {
    Free(sl->text);
  }
  sl->text = NewN( char, strlen(annotation)+1 );
  strcpy( sl->text, annotation );
}


#if !WORDS_BIGENDIAN
void starswap(db_star *st) {
  int i, *wp;
  /* byte-swap x,y,z, dx,dy,dz, magnitude,radius,opacity fields (32-bit float),
   * 		num (32-bit int),
   *		color (16-bit short).
   * group and type fields shouldn't need swapping,
   * assuming the compiler packs bytes into a word in increasing
   * address order.  Seems safe.
   */
  for(i = 0, wp = (int *)st; i < 10; i++)
    wp[i] = htonl(wp[i]);
  st->color = htons(st->color);
}
#endif /*!WORDS_BIGENDIAN*/

void specks_read_sdb( struct stuff *st, char *sdbfname, int timestep )
{
  FILE *inf = fopen(sdbfname, "rb");
  long flen;
  int nspecks, i;
  float min[MAXVAL], max[MAXVAL], sum[MAXVAL];
  struct specklist *sl, **slp;
  register struct speck *sp;
  int dfltvars = (strcmp(st->sdbvars, "mcr") == 0);
  int nvars = strlen(st->sdbvars);

  if(inf == NULL) {
    msg("sdb: %s: cannot open: %s", sdbfname, strerror(errno));
    return;
  }
  /* Just measure file size */
  errno = 0;
  fseek(inf, 0, SEEK_END);
  flen = ftell(inf);
  if(flen == -1 || flen == 0) {
    msg("sdb: %s: can't measure length of file: %s", sdbfname, strerror(errno));

    return;
  }
  nspecks = (flen / sizeof(db_star));

  if(nspecks <= 0) {
    msg("sdb: %s: ignoring empty sdb file", sdbfname);
    return;
  }
  
  sl = NewN(struct specklist, 1);
  memset(sl, 0, sizeof(*sl));
  sl->speckseq = ++st->speckseq;
  
  sl->bytesperspeck = SMALLSPECKSIZE( nvars );
  if(nvars > MAXVAL) nvars = MAXVAL;

  sl->scaledby = st->spacescale;
  sp = NewNSpeck(sl, nspecks);
  sl->specks = sp;
  sl->sel = NewN( SelMask, nspecks );
  sl->nsel = nspecks;
  memset(sl->sel, 0, nspecks*sizeof(SelMask));

  fseek(inf, 0, SEEK_SET);
  for(i = 0; i < nspecks; i++, sp = NextSpeck(sp, sl, 1)) {
    db_star star;
    char *cp;
    int k;
    float *vp;
    if(fread(&star, sizeof(star), 1, inf) <= 0)
	break;
#if !WORDS_BIGENDIAN
    starswap(&star);
#endif
    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->dyn.enabled && st->dyn.trange &&
	(*st->dyn.trange)(&st->dyn, st, tminp, tmaxp)) {
    /* OK */
  } else {
    *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;
  struct mesh **nmesh;

  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 );
	nmesh = NewN( struct mesh *, 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));
	memset(nmesh, 0, needroom * sizeof(*nmesh));
	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) );

	if(d < st->ndata && st->meshes[d])
	    memcpy( nmesh, st->meshes[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->meshes[d] = nmesh;
    }
    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 );
  }
}

static int specks_nonempty_timestep( struct stuff *st, int timestep ) {
  
  if(specks_timespecks( st, st->curdata, timestep ) != NULL)
    return 1;

  if((st->datafile[st->curdata] != NULL &&
		st->datafile[st->curdata][timestep] != NULL))
    return 1;

  if(timestep >= 0 && timestep < st->ntimes &&
		st->meshes[st->curdata][timestep] != NULL)
    return 1;

  if(timestep >= 0 && timestep < st->boxtimes &&
		st->boxes[timestep] != NULL)
    return 1;

  return 0;
}


void specks_set_time( struct stuff *st, double newtime )
{

  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;
    if(st->skipblanktimes) {
	for(nudges = 0; nudges < st->ntimes; nudges++) {
	    ts = (st->curtime + nudges) % st->ntimes;
	    if( specks_nonempty_timestep( st, ts ) )
		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 */
    static Point lastinterest;
    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 );

extern void specks_draw_mesh( struct stuff *st, struct mesh *m, int *texturing );
extern void specks_draw_ellipsoid( struct stuff *st, struct ellipsoid *e );

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;
  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 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];
  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 & RGBBITS;
    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 fast = st->fast;
  int inpick = st->inpick;
  /*int usethresh = st->usethresh&P_USETHRESH ? THRESHBIT : 0;*/
#if USE_PTRACK
  int useptrack = (getenv("PTRACK") != NULL);
#endif
  SelOp seesel = st->seesel;
  float polyminrad, polymaxrad;
  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];