diff options
Diffstat (limited to 'test/ccured_olden/bh/code.c')
| -rw-r--r-- | test/ccured_olden/bh/code.c | 200 |
1 files changed, 0 insertions, 200 deletions
diff --git a/test/ccured_olden/bh/code.c b/test/ccured_olden/bh/code.c deleted file mode 100644 index ee3c0894..00000000 --- a/test/ccured_olden/bh/code.c +++ /dev/null @@ -1,200 +0,0 @@ -/****************************************************************************/ -/* CODE.C: hierarchical N-body code. */ -/* */ -/* Copyright (c) 1993 by Joshua E. Barnes, Honolulu, HI. */ -/* It's free because it's yours. */ -/****************************************************************************/ - -#define global /* nada */ - -#include <stdlib.h> -#include "code.h" - -/* - * Default values for input parameters. - */ - -string defv[] = { - - /* file names for input/output */ - "in=", /* Input file with initial conditions */ - "out=", /* Output file of N-body frames */ - - /* params to control N-body integration */ - "dtime=0.03125", /* Integration time-step */ - "eps=0.025", /* Potential softening parameter */ - "theta=1.0", /* Cell subdivision tolerence */ - "usequad=false", /* If true, use quad moments */ - "rsize=4.0", /* Size of initial root cell */ - "options=", /* Various control options */ - "tstop=2.0", /* Time to stop integration */ - "dtout=0.25", /* Data output interval */ - - /* params used if no input specified to make a Plummer model */ - "nbody=1024", /* Number of particles for test run */ - "seed=123", /* Random number seed for test run */ - - NULL, -}; - -string headline = "Hierarchical N-body Code"; /* default id for run */ - -local void startrun(void); /* initialize system state */ -local void stepsystem(void); /* advance by one time-step */ -local void testdata(void); /* generate test data */ -local void pickshell(vector, real); /* pick point on shell */ - -/* - * MAIN: toplevel routine for hierarchical N-body code. - */ - -int main(int argc, string argv[]) -{ - initparam(argv, defv); /* setup parameter access */ - startrun(); /* set params, input data */ - initoutput(); /* begin system output */ - while (tnow < tstop - 1.0/(1024*freq)) /* while not past tstop */ - stepsystem(); /* advance N-body system */ - stopoutput(); /* finish up output */ - return 0; // Change this to zero when the running time is 3.8s again -} - -/* - * STARTRUN: startup hierarchical N-body code. - */ - -local void startrun(void) -{ - infile = getparam("in"); /* set I/O file names */ - outfile = getparam("out"); - freq = 1.0 / getrparam("dtime"); /* inverse timestep */ - eps = getrparam("eps"); /* softening length */ - theta = getrparam("theta"); /* force accuracy parameter */ - usequad = getbparam("usequad"); /* flag for quad moments */ - rsize = getrparam("rsize"); /* initial size of root */ - options = getparam("options"); /* and control options */ - tstop = getrparam("tstop"); /* time to stop integration */ - freqout = 1.0 / getrparam("dtout"); /* inverse output interval */ - if (! streq(infile, "")) { /* was data file given? */ - - inputdata(); /* read inital data */ - } else { /* make initial conds? */ - nbody = getiparam("nbody"); /* get nbody parameter */ - if (nbody < 1) /* check input value */ - error("startrun: nbody = %d is absurd\n", nbody); - srand((long) getiparam("seed")); /* set random generator */ - testdata(); /* make test model */ - } - nstep = 0; /* start counting steps */ - tout = tnow; /* schedule first output */ -} - -/* - * TESTDATA: generate Plummer model initial conditions for test runs, - * scaled to units such that M = -4E = G = 1 (Henon, Hegge, etc). - * See Aarseth, SJ, Henon, M, & Wielen, R (1974) Astr & Ap, 37, 183. - */ - -#define MFRAC 0.999 /* cut off 1-MFRAC of mass */ - -local void testdata(void) -{ - real rsc, vsc, r, v, x, y; - vector cmr, cmv; - bodyptr p; - - headline = "Hierarchical code: Plummer model"; - /* supply default headline */ - tnow = 0.0; /* reset elapsed model time */ - bodytab = (bodyptr) allocate(nbody * sizeof(body)); - /* alloc space for bodies */ - rsc = 3 * PI / 16; /* set length scale factor */ - vsc = rsqrt(1.0 / rsc); /* and recip. speed scale */ - CLRV(cmr); /* init cm pos, vel */ - CLRV(cmv); - for (p = bodytab; p < bodytab+nbody; p++) { /* loop over particles */ - Type(p) = BODY; /* tag as a body */ - Mass(p) = 1.0 / nbody; /* set masses equal */ - r = 1 / rsqrt(rpow(xrandom(0.0, MFRAC), /* pick r in struct units */ - -2.0/3.0) - 1); - pickshell(Pos(p), rsc * r); /* pick scaled position */ - ADDV(cmr, cmr, Pos(p)); /* add to running sum */ - do { /* select from fn g(x) */ - x = xrandom(0.0, 1.0); /* for x in range 0:1 */ - y = xrandom(0.0, 0.1); /* max of g(x) is 0.092 */ - } while (y > x*x * rpow(1 - x*x, 3.5)); /* using von Neumann tech */ - v = rsqrt(2.0) * x/rpow(1 + r*r, 0.25); /* find v in struct units */ - pickshell(Vel(p), vsc * v); /* pick scaled velocity */ - ADDV(cmv, cmv, Vel(p)); /* add to running sum */ - } - DIVVS(cmr, cmr, (real) nbody); /* normalize cm coords */ - DIVVS(cmv, cmv, (real) nbody); - for (p = bodytab; p < bodytab+nbody; p++) { /* loop over particles */ - SUBV(Pos(p), Pos(p), cmr); /* offset by cm coords */ - SUBV(Vel(p), Vel(p), cmv); - } -} - -/* - * PICKSHELL: pick a random point on a sphere of specified radius. - */ - -local void pickshell(vector vec, real rad) -{ - int k; - real rsq, rsc; - - do { /* pick point in NDIM-space */ - for (k = 0; k < NDIM; k++) /* loop over dimensions */ - vec[k] = xrandom(-1.0, 1.0); /* pick from unit cube */ - DOTVP(rsq, vec, vec); /* compute radius squared */ - } while (rsq > 1.0); /* reject if outside sphere */ - rsc = rad / rsqrt(rsq); /* compute scaling factor */ - MULVS(vec, vec, rsc); /* rescale to radius given */ -} - -/* - * STEPSYSTEM: advance N-body system one time-step. - */ - -local void stepsystem(void) -{ - bodyptr p; - real dt; - vector dvel, dpos; - - if (nstep == 0) { /* about to take 1st step? */ - maketree(bodytab, nbody); /* build tree structure */ - nfcalc = n2bcalc = nbccalc = 0; /* zero counters */ - for (p = bodytab; p < bodytab+nbody; p++) { - /* loop over all bodies */ - hackgrav(p, Mass(p) > 0.0); /* get force on each */ - nfcalc++; /* count force calcs */ - n2bcalc += n2bterm; /* and 2-body terms */ - nbccalc += nbcterm; /* and body-cell terms */ - } - output(); /* do initial output */ - } - dt = 1.0 / freq; /* set basic time-step */ - for (p = bodytab; p < bodytab+nbody; p++) { /* loop over all bodies */ - MULVS(dvel, Acc(p), 0.5 * dt); /* get velocity increment */ - ADDV(Vel(p), Vel(p), dvel); /* advance v by 1/2 step */ - MULVS(dpos, Vel(p), dt); /* get positon increment */ - ADDV(Pos(p), Pos(p), dpos); /* advance r by 1 step */ - } - maketree(bodytab, nbody); /* build tree structure */ - nfcalc = n2bcalc = nbccalc = 0; /* zero counters */ - for (p = bodytab; p < bodytab+nbody; p++) { /* loop over bodies */ - hackgrav(p, Mass(p) > 0.0); /* get force on each */ - nfcalc++; /* count force calcs */ - n2bcalc += n2bterm; /* and 2-body terms */ - nbccalc += nbcterm; /* and body-cell terms */ - } - for (p = bodytab; p < bodytab+nbody; p++) { /* loop over all bodies */ - MULVS(dvel, Acc(p), 0.5 * dt); /* get velocity increment */ - ADDV(Vel(p), Vel(p), dvel); /* advance v by 1/2 step */ - } - nstep++; /* count another time step */ - tnow = tnow + dt; /* finally, advance time */ - output(); /* do major or minor output */ -} |