kira_parti.cc

#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"

static char local_id[] = "$Id$";

extern struct specklist *get_parti( struct stuff *st, double realtime );

typedef worldbundle *worldbundleptr;

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

    int maxstars, maxmarks;	// room allocated for each
    struct specklist *sl;
    struct specklist *marksl;
    int slvalid;
    struct speck *marksp;	// current pointer, updated by add_speck
    int leafcount;		// temp, used in recursion only

    worldstuff() {
	nh = 0;
	wh = NULL;
	ih = 0;
	tmin = 0, tmax = 1;
	tcur = treq = 0;

	treenodes = KIRA_ON;
	treerings = KIRA_OFF;
	treearcs = KIRA_ON;
	ringsizer = KIRA_RINGA;
	ringscale = 1.5;

	ringmin = 2;
	ringmax = 50;

	tickscale = 0.25;
	tracking = wastracking = 0;

	sl = marksl = NULL;
	maxstars = maxmarks = 0;
	slvalid = 0;
	marksp = NULL;
	leafcount = 0;
    }

};


    // speck fields:
    //  [0] worldline index (= kira index for singles, - and unique for others)
    //  [1] mass
    //  [2] number in clump
    //  [3] top-level name
    //  [4] address in clump
    //  [5] ring size
    //  [6] mass ratio (for nonleaves) (mu = child1/(child1+child2))
    //  [7..9] child separation vector (for nonleaves)

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_ROOTID = 3,	// worldline index of root of clump
	SPECK_TREEADDR = 4,	// bit-encoded tree address within our clump (0 for isolated stars)
	SPECK_RINGSIZE = 5,	// size of ring marker
	SPECK_MU = 6,		// mass ratio
	SPECK_SEPVEC = 7,	// separation vector[3]
	SPECK_NDATAFIELDS = 10
};

static char *fieldnames[] = {
	"id",
	"mass",
	"nclump",
	"rootid",
	"treeaddr",
	"ringsize",
	"mu",
	"sepvec3",
	NULL
};

void kira_invalidate( struct stuff *st ) {
    if(st->sl) st->sl->used = -1000000;
    worldstuff *ww = (worldstuff *)st->dyndatadata;
    if(ww) ww->slvalid = 0;
}

int kira_set( struct stuff *st, int what, double val )
{
    struct worldstuff *ww = (struct worldstuff *)st->dyndatadata;
    if(ww == NULL) return 0;

    if(kira_get(st, what) == val)
	return 2;

    switch(what) {
    case KIRA_NODES: ww->treenodes = (int)val; break;
    case KIRA_RINGS: ww->treerings = (int)val; break;
    case KIRA_TREE: ww->treearcs = (int)val; break;
    case KIRA_TICKSCALE: ww->tickscale = val; break;
    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;
    case KIRA_TRACK: if(ww->tracking != (int)val) ww->wastracking = 0;
		     ww->tracking = (int)val; break;
    default: return 0;
    }
    kira_invalidate( st );
    specks_set_timestep( st );
    parti_redraw();
    return 1;
}

double kira_get( struct stuff *st, int what )
{
    struct worldstuff *ww = (struct worldstuff *)st->dyndatadata;
    if(ww == NULL) return 0;

    switch(what) {
    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;
    }
    return 0;
}

int open_parti(struct stuff *st, char *filename, int verbose)
{
    ifstream s(filename);
    if (!s) {
	msg("Data file %s not found.", filename);
	return 0;
    }

    worldbundle *wb;
    struct worldstuff *ww = new worldstuff;
    ww->wh = new worldbundleptr[1024];

    ww->nh = 0;
    while (ww->nh < 1024 && (wb = read_bundle(s, verbose)))
	ww->wh[ww->nh++] = wb;

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

    ww->ih = 0;
    ww->tmin = ww->wh[0]->get_t_min();
    ww->tmax = ww->wh[ww->nh-1]->get_t_max();
    ww->sl = NULL;
    ww->tcur = ww->treq = -1.0;		// invalid

    st->dyndata = 1;
    st->dyndatafunc = get_parti;
    st->dyndatadata = ww;

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

    s.close();
    return 1;
}

double get_parti_time( struct stuff *st ) {
    if(st->dyndatadata == NULL)
	return 0.0;
    return ((struct worldstuff *)(st->dyndatadata))->tcur;
}

void set_parti_time( struct stuff *st, double reqtime ) {
    if(st->dyndatadata == NULL)
	return;
    // otherwise please seek to requested time ...
}

int get_parti_time_range( struct stuff *st, double *tmin, double *tmax )
{
    struct worldstuff *ww = (struct worldstuff *)st->dyndatadata;
    if(ww == NULL) {
	*tmin = 0, *tmax = 1;
	return 0;
    }

    *tmin = ww->tmin;
    *tmax = ww->tmax;
    return 1;
}


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

    sp->val[SPECK_ID] = b->get_index();
    if (!b->is_leaf())
	sp->val[SPECK_ID] = -sp->val[SPECK_ID];	// distinguish CM

    sp->val[SPECK_MASS] = b->get_mass();

    sp->val[SPECK_NCLUMP] = 0;			// complete (add n_leaves) in to_parti

    // Offset unperturbed binaries by 100 in field 2.

    if (b != top) {
	dyn *bb = b;
	while (bb->get_elder_sister())
	    bb = bb->get_elder_sister();
	if (bb->get_kepler()) sp->val[SPECK_NCLUMP] += 100;
    }

    sp->val[SPECK_ROOTID] = top->get_index();

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

    sp->val[SPECK_TREEADDR] = addr;

    sp->val[SPECK_RINGSIZE] = 0;

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

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

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

    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() &&
		(ww->treerings == KIRA_ON || ww->treearcs != KIRA_OFF
		 || (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

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

	    real size = 0;
	    switch(ww->ringsizer) {
	    case KIRA_RINGSEP:
		size = 0.5 * dist;
		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;
	    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];
	}
    }

    // Recursion.

    addr *= 2;
    for_all_daughters(dyn, b, bb)
	sp = assign_specks(bb, top, addr++, sp, sl, ww);

    return sp;
}


struct specklist *to_parti(dyn *root, struct stuff *st, struct specklist *sl)
{
    worldstuff *ww = (worldstuff *)st->dyndatadata;
    struct specklist *marksl = ww->marksl;

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

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

        sl = NewN( struct specklist, 1 );
	memset(sl, 0, sizeof(*sl));
	sl->bytesperspeck = SMALLSPECKSIZE( SPECK_NDATAFIELDS );
	sl->specks = NewNSpeck(sl, ww->maxstars);
	sl->scaledby = 1;
	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;

    }

    struct speck *sp = sl->specks;

    sl->nspecks = 0;
    marksl->nspecks = 0;
    ww->marksp = marksl->specks;
    ww->leafcount = 0;

    if(ww->wastracking) ww->wastracking = -1;

    for_all_daughters(dyn, root, b) {
	int ns = sl->nspecks;
	int mns = marksl->nspecks;
	int nl = ww->leafcount;
	speck *tsp = sp;
	speck *msp = ww->marksp;

	sp = assign_specks(b, b, !b->is_leaf(), sp, sl, ww);

	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
	    tsp = NextSpeck(tsp, sl, 1);
	}
	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);
	}
    }

    if(ww->wastracking < 0) ww->wastracking = 0;	// Detach if tracked pcle not found now

    return sl;
}

struct specklist *get_parti( struct stuff *st, double realtime )
{
    struct worldstuff *ww = (struct worldstuff *)(st->dyndatadata);

    if (!ww)
	return NULL;


    ww->treq = realtime;

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

    if (ww->tcur == realtime && ww->slvalid)
	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])
	;

    dyn *root = create_interpolated_tree(wb, realtime);
    ww->sl = to_parti(root, st, ww->sl);
    rmtree(root);

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

    ww->slvalid = 1;
    return ww->sl;
}

void kira_draw( struct stuff *st, struct specklist *slhead, Matrix *Tc2w, float radperpix )
{
    struct worldstuff *ww = (struct worldstuff *)st->dyndatadata;

    if(ww == NULL)
	return;

    struct specklist *sl;
    float halftickscale = 0.5 * ww->tickscale;

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

    vtfmpoint( &eyepoint, &zero, Tc2w );
    vtfmvector( &fwd, &fwdvec, Tc2w );
    fwdd = -vdot( &eyepoint, &fwd );
    vunit( &unitfwd, &fwd );

#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;

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

#endif /*USE_KIRA*/