Skip to content
Snippets Groups Projects
kira_parti.cc 21.4 KiB
Newer Older
slevy's avatar
slevy committed
#ifdef USE_KIRA

#include "worldline.h"

// #define CAVE 1		// NewN: use potentially-shared-memory malloc
#include "shmem.h"
#include "specks.h"
#include "partiviewc.h"

#include "kira_parti.h"		// declare following as extern "C"
slevy's avatar
slevy committed
#include <sys/types.h>
#include <sys/time.h>
slevy's avatar
slevy committed

slevy's avatar
 
slevy committed
static char local_id[] = "$Id$";
slevy's avatar
slevy committed

slevy's avatar
slevy committed
extern struct specklist *kira_get_parti( struct stuff *st, double realtime );
extern int kira_parse_args( struct dyndata *dd, struct stuff *st, int argc, char **argv );
extern int kira_draw( struct dyndata *dd, struct stuff *st, struct specklist *slhead,
		Matrix *Tc2w, float radperpix );
slevy's avatar
slevy committed

typedef worldbundle *worldbundleptr;
slevy's avatar
 
slevy committed

enum speckfields {
	SPECK_ID = 0,		// worldline index (= kira index for single stars,
				//			unique small negative int for others)
	SPECK_MASS = 1,		// mass
	SPECK_NCLUMP = 2,	// number of stars in clump
	SPECK_TLOG = 3,		// log10(Teff)
	SPECK_LUM = 4,		// L/Lsun?
	SPECK_STYPE = 5,	// star type: 1ms 2cd 3gs 4bh 5sg 6hb
	SPECK_ROOTID = 6,	// worldline index of root of clump
	SPECK_TREEADDR = 7,	// bit-encoded tree address within our clump (0 for isolated stars)
	SPECK_RINGSIZE = 8,	// size of ring marker
	SPECK_MU = 9,		// mass ratio
	SPECK_SEPVEC = 10,	// separation vector[3]
	SPECK_NDATAFIELDS = 13
};

static char *fieldnames[] = {  // Must match ``enum speckfields'' !!
	"id",		// worldline index (=? kira index for singles), unique <0 for others
	"mass",		// mass (fraction of cluster mass)
	"nclump",	// number of stars in clump, = 1 for singles
	"Tlog",		// log10( Teff )
	"Lum",		// L/Lsun?
	"stype",	// stellar type index
	"rootid",	// id of root of clump, = our id for singles
	"treeaddr",	// binary-coded address in clump
	"ringsize",	// size of ring
	"mu",		// mass ratio for nonleaf nodes
	NULL
};

struct vald {
	float min, max, sum;
};

slevy's avatar
slevy committed
struct worldstuff {
slevy's avatar
slevy committed
    int nh, maxnh;
slevy's avatar
slevy committed
    worldbundleptr *wh;
slevy's avatar
slevy committed
    ifstream *ins;
slevy's avatar
slevy committed
    int ih;			// current worldbundle index
    real tmin, tmax;
    real tcur;			// current time
    real treq;			// requested time
slevy's avatar
slevy committed
    int readflags;		// KIRA_VERBOSE | KIRA_READLATER
slevy's avatar
 
slevy committed
    int treenodes;		// KIRA_{OFF|ON|ROOTS}
    int treerings;		// KIRA_{OFF|ON|ROOTS}
slevy's avatar
 
slevy committed
    int treearcs;		// KIRA_{OFF|ON|CROSS|TICK}
    float tickscale;		// size of treearc cross mark (frac of sep)
slevy's avatar
 
slevy committed
    int ringsizer;		// KIRA_RINGSEP, KIRA_RINGA
    float ringscale;		// multiplier for above
    float ringmin, ringmax;	// range of pixel sizes for ring markers
slevy's avatar
 
slevy committed
    int tracking, wastracking;	// id of particle we're tracking, or zero
    Point trackpos;		// last known position of tracked particle

slevy's avatar
 
slevy committed
    int maxstars, maxmarks;	// room allocated for each
slevy's avatar
 
slevy committed
    struct specklist *sl;
slevy's avatar
 
slevy committed
    struct specklist *marksl;
slevy's avatar
 
slevy committed
    int slvalid;
slevy's avatar
slevy committed
    int *starset;		// per-star bitmask of "sets" it belongs to

    struct specklist *trailsl;	// per-star specklist of recent history
    int maxtrail, ntrail;

slevy's avatar
 
slevy committed
    struct speck *marksp;	// current pointer, updated by add_speck
    int leafcount;		// temp, used in recursion only

    struct vald vd[SPECK_NDATAFIELDS];
slevy's avatar
slevy committed

slevy's avatar
 
slevy committed
    worldstuff() {
slevy's avatar
slevy committed
	this->init();
    }

    void init() {
slevy's avatar
 
slevy committed
	nh = 0;
	wh = NULL;
	ih = 0;
	tmin = 0, tmax = 1;
	tcur = treq = 0;

	treenodes = KIRA_ON;
	treerings = KIRA_OFF;
slevy's avatar
 
slevy committed
	treearcs = KIRA_ON;
slevy's avatar
 
slevy committed
	ringsizer = KIRA_RINGA;
	ringscale = 1.5;

	ringmin = 2;
	ringmax = 50;

slevy's avatar
 
slevy committed
	tickscale = 0.25;
slevy's avatar
 
slevy committed
	tracking = wastracking = 0;

slevy's avatar
 
slevy committed
	sl = marksl = NULL;
	maxstars = maxmarks = 0;
slevy's avatar
 
slevy committed
	slvalid = 0;
slevy's avatar
slevy committed

	trailsl = NULL;
	maxtrail = ntrail = 0;

slevy's avatar
 
slevy committed
	marksp = NULL;
	leafcount = 0;
    }

};


void kira_invalidate( struct dyndata *dd, struct stuff *st ) {
slevy's avatar
 
slevy committed
    if(st->sl) st->sl->used = -1000000;
    worldstuff *ww = (worldstuff *)dd->data;
slevy's avatar
 
slevy committed
    if(ww) ww->slvalid = 0;
slevy's avatar
 
slevy committed
}

int kira_set( struct dyndata *dd, struct stuff *st, int what, double val )
slevy's avatar
 
slevy committed
{
    struct worldstuff *ww = (struct worldstuff *)dd->data;
slevy's avatar
 
slevy committed
    if(ww == NULL) return 0;

    if(kira_get(dd, st, what) == val)
slevy's avatar
 
slevy committed
	return 2;

    switch(what) {
    case KIRA_NODES: ww->treenodes = (int)val; break;
    case KIRA_RINGS: ww->treerings = (int)val; break;
slevy's avatar
 
slevy committed
    case KIRA_TREE: ww->treearcs = (int)val; break;
    case KIRA_TICKSCALE: ww->tickscale = val; break;
slevy's avatar
 
slevy committed
    case KIRA_RINGSIZE: ww->ringsizer = (int)val; break;
    case KIRA_RINGSCALE: ww->ringscale = val; break;
    case KIRA_RINGMIN: ww->ringmin = val; break;
    case KIRA_RINGMAX: ww->ringmax = val; break;
slevy's avatar
 
slevy committed
    case KIRA_TRACK: if(ww->tracking != (int)val) ww->wastracking = 0;
		     ww->tracking = (int)val; break;
slevy's avatar
 
slevy committed
    default: return 0;
    }
    kira_invalidate( dd, st );
slevy's avatar
 
slevy committed
    specks_set_timestep( st );
    parti_redraw();
    return 1;
}

double kira_get( struct dyndata *dd, struct stuff *st, int what )
slevy's avatar
 
slevy committed
{
    struct worldstuff *ww = (struct worldstuff *)dd->data;
slevy's avatar
 
slevy committed
    if(ww == NULL) return 0;

    switch(what) {
slevy's avatar
 
slevy committed
    case KIRA_NODES: return ww->treenodes;
    case KIRA_RINGS: return ww->treerings;
    case KIRA_TREE:  return ww->treearcs;
    case KIRA_TICKSCALE: return ww->tickscale;
    case KIRA_RINGSIZE: return ww->ringsizer;
    case KIRA_RINGSCALE: return ww->ringscale;
    case KIRA_RINGMIN: return ww->ringmin;
    case KIRA_RINGMAX: return ww->ringmax;
    case KIRA_TRACK: return ww->tracking;
slevy's avatar
 
slevy committed
    }
    return 0;
}

slevy's avatar
slevy committed
int hasdata( ifstream *s, double maxwait ) {
    struct timeval tv;
    fd_set fds;
    if(s == NULL) return 0;
    if(s->rdbuf()->in_avail() || s->eof()) return 1;
    if(maxwait < 0) maxwait = 0;
    tv.tv_sec = (int)maxwait;
    tv.tv_usec = (int) (1000000 * (maxwait - tv.tv_sec));
    FD_ZERO(&fds);
    FD_SET(s->rdbuf()->fd(), &fds);
    return(select(s->rdbuf()->fd() + 1, &fds, NULL, NULL, &tv) > 0);
}

int kira_help( struct dyndata *dd, struct stuff *st, int verbose )
{
    msg(" kira {node|ring|size|scale|span|track}  starlab controls; try \"kira ?\"");
    return 1;
}

int kira_read_more( struct dyndata *dd, struct stuff *st, int nmax, double tmax, double realdtmax )
slevy's avatar
slevy committed
{
    struct worldstuff *ww = (struct worldstuff *)dd->data;
slevy's avatar
slevy committed
    if(ww == NULL) return -1;
    ifstream *ins = ww->ins;
    if(ins == NULL || !ins->rdbuf()->is_open()) return -1;
slevy's avatar
slevy committed

    int paras = 0;
    int timelimit = (realdtmax < 1e6);
    double wallstop = realdtmax + wallclock_time();
    while( paras < nmax && (ww->nh == 0 || ww->wh[ww->nh-1]->get_t_max() < tmax) ) {

	if(timelimit) {
	    double dt = wallstop - wallclock_time();
	    if(!hasdata(ins, dt))
		break;
	    if(dt <= 0)
		break;
	}
slevy's avatar
slevy committed

slevy's avatar
slevy committed
	if(ww->nh >= ww->maxnh) {
	    ww->maxnh = (ww->nh > ww->maxnh*2) ? ww->nh + 1 : ww->maxnh*2;
	    ww->wh = RenewN( ww->wh, worldbundleptr, ww->maxnh );
	}
slevy's avatar
slevy committed

slevy's avatar
slevy committed
	worldbundle *wb;

	wb = read_bundle(*ins, ww->readflags & KIRA_VERBOSE);
	if(wb == NULL) {
	    ins->close();
	    if(paras == 0) paras = -1;
	    break;
	}
slevy's avatar
slevy committed
	ww->wh[ww->nh++] = wb;
slevy's avatar
slevy committed
	paras++;
    }

    if(paras > 0) {
	ww->tmin = ww->wh[0]->get_t_min();
	ww->tmax = ww->wh[ww->nh-1]->get_t_max();
    }
    return paras;
}
slevy's avatar
slevy committed

int kira_get_trange( struct dyndata *dd, struct stuff *st, double *tminp, double *tmaxp )
    struct worldstuff *ww = (struct worldstuff *)dd->data;
    if(ww == NULL) {
	*tminp = *tmaxp = 0;
	return 0;
    }
    *tminp = ww->tmin;
    *tmaxp = ww->tmax;
    return 1;
}

int kira_open( struct dyndata *dd, struct stuff *st, char *filename, int readflags)
slevy's avatar
slevy committed
{
    ifstream *ins = new ifstream(filename);
    if (!ins || !*ins) {
	msg("Kira/starlab data file %s not found.", filename);
	delete ins;
slevy's avatar
slevy committed
	return 0;
    }

slevy's avatar
slevy committed
    struct worldstuff *ww = NewN( worldstuff, 1 );	// from shmem -- avoid "new"
    ww->init();
    ww->maxnh = 1024;
    ww->wh = NewN( worldbundleptr, ww->maxnh );
    ww->ins = ins;
    ww->tmin = ww->tmax = 0;
    ww->tcur = ww->treq = -1.0;		/* invalidate */
slevy's avatar
slevy committed
    ww->readflags = readflags;
slevy's avatar
slevy committed

    memset( dd, 0, sizeof(*dd) );
    dd->data = ww;
    dd->getspecks = kira_get_parti;
    dd->trange = kira_get_trange;
    dd->ctlcmd = kira_parse_args;
    dd->draw = kira_draw;
    dd->help = kira_help;
    dd->free = NULL;			/* XXX create a 'free' function someday */
    dd->enabled = 1;
slevy's avatar
slevy committed

slevy's avatar
slevy committed
    ww->nh = 0;
    if(! (readflags & KIRA_READLATER)) {
	kira_read_more( dd, st, ww->maxnh, 1e38, 1e38 );
slevy's avatar
slevy committed

	if (ww->nh == 0) {
	    msg("Data file %s not in worldbundle format.", filename);
	    return 0;
	}
    }

    ww->ih = 0;
    ww->sl = NULL;

    // If data field names aren't already assigned, set them now.
    for(int i = 0; fieldnames[i] != NULL; i++) {
	if(st->vdesc[st->curdata][i].name[0] == '\0')
	    strcpy(st->vdesc[st->curdata][i].name, fieldnames[i]);
    }
slevy's avatar
slevy committed
    return 1;
}

double get_parti_time( struct dyndata *dd, struct stuff *st ) {
    return dd->data ? ((struct worldstuff *)(dd->data))->tcur : 0.0;
slevy's avatar
slevy committed
}

void set_parti_time( struct dyndata *dd, struct stuff *st, double reqtime ) {
    if(dd->data == NULL)
slevy's avatar
slevy committed
	return;
    // otherwise please seek to requested time ...
}

static float log10_of( float v ) {
    return v <= 0 ? 0 : log10f( v );
}
slevy's avatar
 
slevy committed

speck *add_speck(pdyn *b, pdyn *top, int addr, speck *sp, specklist *sl, worldstuff *ww)
slevy's avatar
slevy committed
{
slevy's avatar
 
slevy committed
    sp->p.x[0] = b->get_pos()[0];
    sp->p.x[1] = b->get_pos()[1];
    sp->p.x[2] = b->get_pos()[2];

slevy's avatar
 
slevy committed
    sp->val[SPECK_ID] = b->get_index();
slevy's avatar
 
slevy committed
    if (!b->is_leaf())
slevy's avatar
 
slevy committed
	sp->val[SPECK_ID] = -sp->val[SPECK_ID];	// distinguish CM
slevy's avatar
 
slevy committed

slevy's avatar
 
slevy committed
    sp->val[SPECK_MASS] = b->get_mass();
    sp->val[SPECK_TLOG] = log10_of( b->get_temperature() );
    sp->val[SPECK_LUM] = b->get_luminosity();
    sp->val[SPECK_STYPE] = b->get_stellar_type();	/* see starlab/inc/star_support.h */
slevy's avatar
 
slevy committed

slevy's avatar
slevy committed
    sp->val[SPECK_NCLUMP] = 0;		// complete (add n_leaves) in kira_to_parti
slevy's avatar
 
slevy committed

    // Offset unperturbed binaries by 100 in field 2.

    if (b != top) {
	pdyn *bb = b;
slevy's avatar
 
slevy committed
	while (bb->get_elder_sister())
	    bb = bb->get_elder_sister();
	if (bb->get_kepler()) sp->val[SPECK_NCLUMP] += 100;
slevy's avatar
 
slevy committed
    }

slevy's avatar
 
slevy committed
    sp->val[SPECK_ROOTID] = top->get_index();
slevy's avatar
 
slevy committed

    // Addr is 0 for top-level leaf, 1 for top-level CM,
    // constructed from parent addr otherwise.

slevy's avatar
 
slevy committed
    sp->val[SPECK_TREEADDR] = addr;

    sp->val[SPECK_RINGSIZE] = 0;
slevy's avatar
 
slevy committed

    sl->nspecks++;
    sp = NextSpeck(sp, sl, 1);
slevy's avatar
 
slevy committed
    return sp;
}
    
speck *assign_specks(pdyn *b, pdyn *top, int addr, speck *sp, specklist *sl, worldstuff *ww)
slevy's avatar
 
slevy committed
{
    // Convert relative coordinates to absolute.

    if (b != top)
	b->inc_pos(b->get_parent()->get_pos());

slevy's avatar
 
slevy committed
    if(ww->tracking != 0 && abs(ww->tracking) == b->get_index()) {
	Point newpos;
	newpos.x[0] = b->get_pos()[0];
	newpos.x[1] = b->get_pos()[1];
	newpos.x[2] = b->get_pos()[2];
	if(ww->wastracking) {
	    Point delta;
	    vsub( &delta, &newpos,&ww->trackpos );
	    parti_nudge_camera( &delta );
	}
	ww->trackpos = newpos;
	ww->wastracking = 1;
    }

slevy's avatar
 
slevy committed
    if(b->is_leaf()
		|| ww->treenodes == KIRA_ON
		|| (ww->treenodes == KIRA_ROOTS && b == top)) {
	sp = add_speck(b, top, addr, sp, sl, ww);

	if(b->is_leaf())
	    ww->leafcount++;
    }

    if(!b->is_leaf() &&
slevy's avatar
 
slevy committed
		(ww->treerings == KIRA_ON || ww->treearcs != KIRA_OFF
slevy's avatar
 
slevy committed
		 || (ww->treerings == KIRA_ROOTS && b == top))) {

	struct speck *marksp = ww->marksp;
	ww->marksp = add_speck(b, top, addr, marksp, ww->marksl, ww);

	// find our two children

	pdyn *b1 = b->get_oldest_daughter();
	pdyn *b2 = b1 ? b1->get_younger_sister() : NULL;
slevy's avatar
 
slevy committed
	if(b2) {
	    vector sep = b1->get_pos() - b2->get_pos();
	    real dist = sqrt(sep * sep);
slevy's avatar
 
slevy committed
	    real mass = b->get_mass();	// = sum of b1&b2 masses.
slevy's avatar
 
slevy committed

	    real size = 0;
	    switch(ww->ringsizer) {
	    case KIRA_RINGSEP:
slevy's avatar
 
slevy committed
		size = 0.5 * dist;
slevy's avatar
 
slevy committed
		break;
	    case KIRA_RINGA:
	      {
		vector vel = b1->get_vel() - b2->get_vel();
		real speed2 = vel * vel;
		real E = 0.5 * speed2 - mass / dist;
		size = -mass / (2*E);	// semimajor axis
	      }
	    }
	    marksp->val[SPECK_RINGSIZE] = size * ww->ringscale;
slevy's avatar
 
slevy committed
	    marksp->val[SPECK_MU] = b1->get_mass() / mass;
	    marksp->val[SPECK_SEPVEC] = sep[0];
	    marksp->val[SPECK_SEPVEC+1] = sep[1];
	    marksp->val[SPECK_SEPVEC+2] = sep[2];
slevy's avatar
 
slevy committed
	}
    }
slevy's avatar
 
slevy committed

    // Recursion.

    addr *= 2;
    for_all_daughters(pdyn, b, bb)
slevy's avatar
 
slevy committed
	sp = assign_specks(bb, top, addr++, sp, sl, ww);
slevy's avatar
 
slevy committed

    return sp;
slevy's avatar
slevy committed
}

slevy's avatar
 
slevy committed

struct specklist *kira_to_parti(pdyn *root, struct dyndata *dd, struct stuff *st, struct specklist *sl)
slevy's avatar
slevy committed
{
    worldstuff *ww = (worldstuff *)dd->data;
slevy's avatar
slevy committed

    if(ww == NULL || ww->nh == 0) return NULL;

slevy's avatar
 
slevy committed
    struct specklist *marksl = ww->marksl;

    if(sl == NULL || sl != ww->sl) {
slevy's avatar
slevy committed
	// Use NewN to allocate these in (potentially) shared memory,
	// rather than new which will allocate in local malloc pool.
	// Needed for virdir/cave version.
slevy's avatar
 
slevy committed

	int smax = root->n_leaves();
	ww->maxstars = 2*smax;
	ww->maxmarks = 1+smax;

slevy's avatar
slevy committed
        sl = NewN( struct specklist, 1 );
	memset(sl, 0, sizeof(*sl));
slevy's avatar
 
slevy committed
	sl->bytesperspeck = SMALLSPECKSIZE( SPECK_NDATAFIELDS );
	sl->specks = NewNSpeck(sl, ww->maxstars);
slevy's avatar
slevy committed
	sl->scaledby = 1;
slevy's avatar
 
slevy committed
	marksl = NewN( struct specklist, 1 );
	*marksl = *sl;
	marksl->specks = NewNSpeck(marksl, ww->maxmarks);
	marksl->special = MARKERS;

	sl->next = marksl;
	marksl->next = NULL;

	ww->sl = sl;
	ww->marksl = marksl;

slevy's avatar
slevy committed
    }
slevy's avatar
 
slevy committed

slevy's avatar
slevy committed
    struct speck *sp = sl->specks;

slevy's avatar
 
slevy committed
    sl->nspecks = 0;
    sl->colorseq = sl->sizeseq = sl->threshseq = 0;
slevy's avatar
 
slevy committed
    marksl->nspecks = 0;
    ww->marksp = marksl->specks;
    ww->leafcount = 0;
slevy's avatar
 
slevy committed

    int i;
    struct vald *vd = ww->vd;
    for(i = 0, vd = ww->vd; i < SPECK_NDATAFIELDS; i++, vd++) {
	vd->min = 1e9;
	vd->max = -1e9;
	vd->sum = 0;
    }
    int ntotal = 0;

slevy's avatar
 
slevy committed
    if(ww->wastracking) ww->wastracking = -1;

    for_all_daughters(pdyn, root, b) {
slevy's avatar
 
slevy committed
	int ns = sl->nspecks;
slevy's avatar
 
slevy committed
	int mns = marksl->nspecks;
	int nl = ww->leafcount;
slevy's avatar
 
slevy committed
	speck *tsp = sp;
slevy's avatar
 
slevy committed
	speck *msp = ww->marksp;
slevy's avatar
 
slevy committed

slevy's avatar
 
slevy committed
	sp = assign_specks(b, b, !b->is_leaf(), sp, sl, ww);
slevy's avatar
 
slevy committed

slevy's avatar
 
slevy committed
	int k;
	int nleaves = ww->leafcount - nl;
	int count = sl->nspecks - ns;

	for(k = 0; k < count; k++) {
	    tsp->val[SPECK_NCLUMP] += nleaves;
	    if(tsp->val[0] < 0)
		tsp->val[SPECK_NCLUMP] = -tsp->val[SPECK_NCLUMP]; // negate nclump for CM nodes
	    for(i = 0, vd = ww->vd; i <= SPECK_STYPE; i++, vd++) {
		float v = tsp->val[i];
		if(vd->min > v) vd->min = v;
		if(vd->max < v) vd->max = v;
		vd->sum += v;
	    }
	    ntotal++;
slevy's avatar
 
slevy committed
	    tsp = NextSpeck(tsp, sl, 1);
slevy's avatar
slevy committed
	}
slevy's avatar
 
slevy committed
	count = marksl->nspecks - mns;
	for(k = 0; k < count; k++) {
	    msp->val[SPECK_NCLUMP] += nleaves;
	    if(msp->val[0] < 0)
		msp->val[SPECK_NCLUMP] = -msp->val[SPECK_NCLUMP]; // negate nclump for CM nodes
	    msp = NextSpeck(msp, sl, 1);
	}
slevy's avatar
slevy committed
    }
slevy's avatar
 
slevy committed

slevy's avatar
 
slevy committed
    if(ww->wastracking < 0) ww->wastracking = 0;	// Detach if tracked pcle not found now

    for(i = 0, vd = ww->vd; i < SPECK_NDATAFIELDS && i < MAXVAL; i++, vd++) {
	struct valdesc *tvd = &st->vdesc[ st->curdata ][ i ];

	tvd->nsamples = ntotal;
	if(vd->min > vd->max) {
	    vd->min = vd->max = 0;
	    tvd->nsamples = 0;
	}
	tvd->min = vd->min;
	tvd->max = vd->max;
	tvd->sum = vd->sum;
	tvd->mean = tvd->nsamples > 0 ? vd->sum / tvd->nsamples : 0;
    }

slevy's avatar
slevy committed
    return sl;
}

struct specklist *kira_get_parti( struct dyndata *dd, struct stuff *st, double realtime )
slevy's avatar
slevy committed
{
    struct worldstuff *ww = (struct worldstuff *)dd->data;
slevy's avatar
slevy committed

slevy's avatar
slevy committed
    if (ww == NULL || ww->nh == 0)
slevy's avatar
slevy committed
	return NULL;

    ww->treq = realtime;

    if (realtime < ww->tmin) realtime = ww->tmin;
    if (realtime > ww->tmax) realtime = ww->tmax;

slevy's avatar
 
slevy committed
    if (ww->tcur == realtime && ww->slvalid)
slevy's avatar
slevy committed
	return ww->sl;

    int ih = ww->ih;
    int nh = ww->nh;

    worldbundle *wb = ww->wh[ih];
    for(; realtime > wb->get_t_max() && ih < nh-1; ih++, wb = ww->wh[ih])
	;
    for(; realtime < wb->get_t_min() && ih > 0; ih--, wb = ww->wh[ih])
	;

#if OLDTREE
    pdyn *root = create_interpolated_tree(wb, realtime);
    ww->sl = kira_to_parti(root, dd, st, ww->sl);
    rmtree(root);
#else
    pdyn *root = create_interpolated_tree2(wb, realtime);
    ww->sl = kira_to_parti(root, dd, st, ww->sl);
#endif
slevy's avatar
slevy committed

    ww->ih = ih;
    ww->tcur = realtime;

slevy's avatar
 
slevy committed
    ww->slvalid = 1;
slevy's avatar
slevy committed
    return ww->sl;
}

int kira_draw( struct dyndata *dd, struct stuff *st, struct specklist *slhead, Matrix *Tc2w, float radperpix )
slevy's avatar
 
slevy committed
{
    struct worldstuff *ww = (struct worldstuff *)dd->data;
slevy's avatar
 
slevy committed

    if(ww == NULL || !dd->enabled)
	return 0;
slevy's avatar
 
slevy committed

slevy's avatar
 
slevy committed
    struct specklist *sl;
slevy's avatar
 
slevy committed
    float halftickscale = 0.5 * ww->tickscale;
slevy's avatar
 
slevy committed

    static Point zero = {0,0,0};
    Point fwdvec = {0,0,-radperpix};	// scaled by pixels-per-radian
slevy's avatar
 
slevy committed
    Point eyepoint, fwd, unitfwd;
slevy's avatar
 
slevy committed
    float fwdd;

    vtfmpoint( &eyepoint, &zero, Tc2w );
    vtfmvector( &fwd, &fwdvec, Tc2w );
    fwdd = -vdot( &eyepoint, &fwd );
slevy's avatar
 
slevy committed
    vunit( &unitfwd, &fwd );
slevy's avatar
 
slevy committed

#define MAXRING 32
    Point fan[MAXRING];
    int nring = 24;
    int i;
    for(i = 0; i < nring; i++) {
	float th = 2*M_PI*i / nring;
	vcomb( &fan[i],  cos(th),(Point *)&Tc2w->m[0*4+0],
			 sin(th),(Point *)&Tc2w->m[1*4+0] );
    }

    for(sl = slhead; sl != NULL; sl = sl->next) {
	if(sl->special == MARKERS) {
	    int ns = sl->nspecks;
	    struct speck *sp = sl->specks;
	    for(i = 0; i < ns; i++, sp = NextSpeck(sp, sl, 1)) {
		float unitperpix = vdot( &fwd, &sp->p ) + fwdd;
		if(unitperpix <= 0) continue;
slevy's avatar
 
slevy committed

		if(ww->treerings == KIRA_ON ||
			(ww->treerings == KIRA_ROOTS && sp->val[SPECK_TREEADDR] == 1)) {
		    float ringpixels = fabs(sp->val[SPECK_RINGSIZE]) / unitperpix;
		    if(ringpixels < ww->ringmin) ringpixels = ww->ringmin;
		    if(ringpixels > ww->ringmax) ringpixels = ww->ringmax;
		    float rring = ringpixels * unitperpix;
		    int step = ringpixels>20 ? 1 : ringpixels>8 ? 2 : 3;

		    glColor3ubv( (GLubyte *)&sp->rgba );
		    glBegin( GL_LINE_LOOP );
		    for(int k = 0; k < nring; k+=step)
			glVertex3f( sp->p.x[0] + rring*fan[k].x[0],
				    sp->p.x[1] + rring*fan[k].x[1],
				    sp->p.x[2] + rring*fan[k].x[2] );
		    glEnd();
		}
		if(ww->treearcs != KIRA_OFF) {
		    float mu = sp->val[SPECK_MU];
		    float *sep = &sp->val[SPECK_SEPVEC];
		    glColor3ubv( (GLubyte *)&sp->rgba );
		    glBegin( GL_LINES );

		    if(ww->treearcs != KIRA_TICK) {
			glVertex3f(
			    sp->p.x[0] - mu*sep[0],
			    sp->p.x[1] - mu*sep[1],
			    sp->p.x[2] - mu*sep[2] );
			glVertex3f(
			    sp->p.x[0] + (1-mu)*sep[0],
			    sp->p.x[1] + (1-mu)*sep[1],
			    sp->p.x[2] + (1-mu)*sep[2] );
		    }

		    if(ww->treearcs == KIRA_CROSS || ww->treearcs == KIRA_TICK) {
			Point cr;
			vcross( &cr,  &unitfwd, (Point *)sep );
			float sepsep = vdot((Point *)sep, (Point *)sep);
			float sepfwd = vdot((Point *)sep, &unitfwd);
			float scaleby = halftickscale / sqrt(1 - sepfwd*sepfwd/sepsep);
			glVertex3f(
			   sp->p.x[0] - scaleby*cr.x[0],
			   sp->p.x[1] - scaleby*cr.x[1],
			   sp->p.x[2] - scaleby*cr.x[2] );
			glVertex3f(
			   sp->p.x[0] + scaleby*cr.x[0],
			   sp->p.x[1] + scaleby*cr.x[1],
			   sp->p.x[2] + scaleby*cr.x[2] );
		    }
		    glEnd();
		}
slevy's avatar
 
slevy committed
	    }
	}
    }
slevy's avatar
slevy committed

    glPointSize( 1.0 );
    glBegin( GL_POINTS );
    for(i = 1; i < st->nghosts; i++) {
	sl = specks_timespecks( st, st->curdata, i );
	if(sl == NULL) continue;
	int k, rgba = 0;
	struct speck *sp;
	for(k = 0, sp = sl->specks; k < sl->nspecks; k++, sp = NextSpeck(sp, sl, 1)) {
	    if(rgba != sp->rgba) {
		glColor3ubv( (GLubyte *)&sp->rgba );
		rgba = sp->rgba;
	    }
	    glVertex3fv( &sp->p.x[0] );
	}
    }
    glEnd();
    return 1;
}

int kira_parse_args( struct dyndata *dd, struct stuff *st, int argc, char **argv )
{
    char *swhat = argv[1];
    char *sval = argc>2 ? argv[2] : NULL;
    int what;
    double val;

    if(0!=strncmp(argv[0], "kira", 4))	/* accept "kira" or "kiractl" */
	return 0;
slevy's avatar
slevy committed

    if(swhat == NULL) swhat = "?";

    if(!strncmp(swhat, "sep", 3) || !strncmp(swhat, "semi", 4)) {
	kira_set( dd, st, KIRA_RINGSIZE, swhat[2]=='p' ? KIRA_RINGSEP : KIRA_RINGA );
	msg("kiractl ringsize %s",
	    kira_get( dd, st, KIRA_RINGSIZE ) == KIRA_RINGSEP ? "separation" : "semimajor");

    } else if(!strcmp(swhat, "ringsize") || !strcmp(swhat, "ringscale") || !strcmp(swhat, "size")) {
	if(sval) {
	    val = !strncmp(sval,"sep",3) ? KIRA_RINGSEP
		: !strncmp(sval,"semi",4) ? KIRA_RINGA
		: !strcmp(sval,"a") ? KIRA_RINGA
		: kira_get( dd, st, KIRA_RINGSIZE );
	    kira_set( dd, st, KIRA_RINGSIZE, val );
	    if(argc > 3)
		kira_set( dd, st, KIRA_RINGSCALE,
		    getfloat( argv[3], kira_get( dd, st, KIRA_RINGSCALE ) ) );
	}
	msg("kiractl ringsize %s %g",
	    kira_get( dd, st, KIRA_RINGSIZE ) == KIRA_RINGSEP ? "separation" : "semimajor",
	    kira_get( dd, st, KIRA_RINGSCALE ));

    } else if(!strcmp(swhat, "ringscale") || !strcmp(swhat, "scale")) {
	if(sval)
	    kira_set( dd, st, KIRA_RINGSCALE, getfloat( sval, kira_get(dd,st,KIRA_RINGSCALE) ) );
	msg("kiractl ringscale %g", kira_get(dd,st,KIRA_RINGSCALE));

    } else if(!strcmp(swhat, "span") || !strcmp(swhat, "ringspan")) {
	if(argc>2 && (val = getfloat( argv[2], -1 )) >= 0)
		kira_set( dd, st, KIRA_RINGMIN, val );
	if(argc>3 && (val = getfloat( argv[3], -1 )) >= 0)
		kira_set( dd, st, KIRA_RINGMAX, val );
	msg("kiractl ringspan %.0f %.0f",
	    kira_get( dd, st, KIRA_RINGMIN ), kira_get( dd, st, KIRA_RINGMAX ));

    } else if(!strncmp(swhat, "nod", 3) || !strncmp(swhat, "ring", 4)) {
	what = swhat[0]=='n' ? KIRA_NODES : KIRA_RINGS;
	if(sval)
	    kira_set( dd, st, what,
		 (0==strncmp(sval, "root", 4)) ? KIRA_ROOTS : getbool(sval, KIRA_ON) );
	val = kira_get( dd, st, what );
	msg("kiractl %s %s %g",
	    what==KIRA_NODES ? "nodes" : "rings",
	    val==2 ? "root" : val==1 ? "on" : "off",
	    kira_get( dd, st, KIRA_TICKSCALE ));

    } else if(!strncmp(swhat, "tree", 3) || !strcmp(swhat, "arc")) {
	if(sval)
	    kira_set( dd, st, KIRA_TREE,
		sval[0]=='c' ? KIRA_CROSS : sval[0]=='t' ? KIRA_TICK
		: getbool(sval, KIRA_ON) );
	if(argc>3 && sscanf(argv[3], "%lf", &val)>0)
	    kira_set( dd, st, KIRA_TICKSCALE, val );
	val = kira_get( dd, st, KIRA_TREE );
	msg("kiractl tree %s",
	    val==KIRA_CROSS ? "cross" : val==KIRA_TICK ? "tick"
	    : val ? "on" : "off");

    } else if(!strncmp(swhat, "track", 2)) {
	if(sval)
	    kira_set( dd, st, KIRA_TRACK, getbool(sval, 0) );
	val = kira_get( dd, st, KIRA_TRACK );
	msg(val == 0 ? "kiractl track off" : "kiractl track %d", (int)val);

    } else {
	msg("kiractl {node|ring} {on|off|root} | tree {on|off|cross|tick} [<tickscale>] | size {sep|semimaj} | scale <fac> | span <minpix> <maxpix> | track <id>");
    }
    return 1;
slevy's avatar
 
slevy committed
}

#endif /* USE_KIRA */