Newer
Older
#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"
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 );
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
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;
};
int ih; // current worldbundle index
real tmin, tmax;
real tcur; // current time
real treq; // requested time
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 *starset; // per-star bitmask of "sets" it belongs to
struct specklist *trailsl; // per-star specklist of recent history
int maxtrail, ntrail;
struct speck *marksp; // current pointer, updated by add_speck
int leafcount; // temp, used in recursion only
nh = 0;
wh = NULL;
ih = 0;
tmin = 0, tmax = 1;
tcur = treq = 0;
treenodes = KIRA_ON;
treerings = KIRA_OFF;
ringsizer = KIRA_RINGA;
ringscale = 1.5;
ringmin = 2;
ringmax = 50;
void kira_invalidate( struct dyndata *dd, struct stuff *st ) {
worldstuff *ww = (worldstuff *)dd->data;
int kira_set( struct dyndata *dd, struct stuff *st, int what, double val )
struct worldstuff *ww = (struct worldstuff *)dd->data;
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;
double kira_get( struct dyndata *dd, struct stuff *st, int what )
struct worldstuff *ww = (struct worldstuff *)dd->data;
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;
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 )
struct worldstuff *ww = (struct worldstuff *)dd->data;
if(ww == NULL) return -1;
ifstream *ins = ww->ins;
if(ins == NULL || !ins->rdbuf()->is_open()) return -1;
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;
}
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 );
}
worldbundle *wb;
wb = read_bundle(*ins, ww->readflags & KIRA_VERBOSE);
if(wb == NULL) {
ins->close();
if(paras == 0) paras = -1;
break;
}
paras++;
}
if(paras > 0) {
ww->tmin = ww->wh[0]->get_t_min();
ww->tmax = ww->wh[ww->nh-1]->get_t_max();
}
return paras;
}
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)
{
ifstream *ins = new ifstream(filename);
if (!ins || !*ins) {
msg("Kira/starlab data file %s not found.", filename);
delete ins;
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 */
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;
ww->nh = 0;
if(! (readflags & KIRA_READLATER)) {
kira_read_more( dd, st, ww->maxnh, 1e38, 1e38 );
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]);
}
double get_parti_time( struct dyndata *dd, struct stuff *st ) {
return dd->data ? ((struct worldstuff *)(dd->data))->tcur : 0.0;
void set_parti_time( struct dyndata *dd, struct stuff *st, double reqtime ) {
if(dd->data == NULL)
static float log10_of( float v ) {
return v <= 0 ? 0 : log10f( v );
}
speck *add_speck(pdyn *b, pdyn *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_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 */
sp->val[SPECK_NCLUMP] = 0; // complete (add n_leaves) in kira_to_parti
while (bb->get_elder_sister())
bb = bb->get_elder_sister();
if (bb->get_kepler()) sp->val[SPECK_NCLUMP] += 100;
// Addr is 0 for top-level leaf, 1 for top-level CM,
// constructed from parent addr otherwise.
speck *assign_specks(pdyn *b, pdyn *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_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;
if(b2) {
vector sep = b1->get_pos() - b2->get_pos();
real dist = sqrt(sep * sep);
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];
struct specklist *kira_to_parti(pdyn *root, struct dyndata *dd, struct stuff *st, struct specklist *sl)
worldstuff *ww = (worldstuff *)dd->data;
// 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->bytesperspeck = SMALLSPECKSIZE( SPECK_NDATAFIELDS );
sl->specks = NewNSpeck(sl, ww->maxstars);
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;
sl->colorseq = sl->sizeseq = sl->threshseq = 0;
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;
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++;
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
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;
}
struct specklist *kira_get_parti( struct dyndata *dd, struct stuff *st, double realtime )
struct worldstuff *ww = (struct worldstuff *)dd->data;
return NULL;
ww->treq = realtime;
if (realtime < ww->tmin) realtime = ww->tmin;
if (realtime > ww->tmax) realtime = ww->tmax;
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
int kira_draw( struct dyndata *dd, struct stuff *st, struct specklist *slhead, Matrix *Tc2w, float radperpix )
struct worldstuff *ww = (struct worldstuff *)dd->data;
if(ww == NULL || !dd->enabled)
return 0;
static Point zero = {0,0,0};
Point fwdvec = {0,0,-radperpix}; // scaled by pixels-per-radian
float fwdd;
vtfmpoint( &eyepoint, &zero, Tc2w );
vtfmvector( &fwd, &fwdvec, Tc2w );
fwdd = -vdot( &eyepoint, &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;
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
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();
}
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;
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
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;