```/****************************************************************************/
/* TREEIO.C: I/O routines for hierarchical N-body code. Public routines:    */
/* inputdata(), startoutput(), output(), savestate(), restorestate().       */
/* Copyright (c) 2001 by Joshua E. Barnes, Honolulu, Hawai`i.               */
/****************************************************************************/

#include "stdinc.h"
#include "mathfns.h"
#include "vectmath.h"
#include "getparam.h"
#include "treecode.h"

#include <sys/types.h>
#include <sys/stat.h>
#include <strings.h>

/*
* Prototypes for local routines.
*/

local void outputdata(void);                    /* write N-body data        */
local void diagnostics(void);                   /* eval N-body diagnostics  */
local void in_int(stream, int *);               /* input integer value      */
local void in_real(stream, real *);             /* input real value         */
local void in_vector(stream, vector);           /* input vector of reals    */
local void out_int(stream, int);                /* output integer value     */
local void out_real(stream, real);              /* output real value        */
local void out_vector(stream, vector);          /* output vector of reals   */

/*
* Diagnositc output variables.
*/

local real mtot;                                /* total mass of system     */
local real etot[3];                             /* Etot, KE, PE of system   */
local matrix keten;                             /* kinetic energy tensor    */
local matrix peten;                             /* potential energy tensor  */
local vector cmpos;                             /* center of mass position  */
local vector cmvel;                             /* center of mass velocity  */
local vector amvec;                             /* angular momentum vector  */

/*
* INPUTDATA: read initial conditions from input file.
*/

void inputdata(void)
{
stream instr;
int ndim;
bodyptr p;

instr = stropen(infile, "r");               /* open input stream        */
in_int(instr, &nbody);                      /* read number of bodies    */
if (nbody < 1)
error("inputdata: nbody = %d is absurd\n", nbody);
in_int(instr, &ndim);                       /* read number of dims      */
if (ndim != NDIM)
error("inputdata: ndim = %d; expected %d\n", ndim, NDIM);
in_real(instr, &tnow);                      /* read starting time       */
bodytab = (bodyptr) allocate(nbody * sizeof(body));
/* allocate body array      */
for (p = bodytab; p < bodytab+nbody; p++)   /* loop over all bodies     */
in_real(instr, &Mass(p));               /* read mass of each        */
for (p = bodytab; p < bodytab+nbody; p++)
in_vector(instr, Pos(p));               /* read position of each    */
for (p = bodytab; p < bodytab+nbody; p++)
in_vector(instr, Vel(p));               /* read velocity of each    */
fclose(instr);                              /* close input stream       */
if (scanopt(options, "reset-time"))         /* reset starting time?     */
tnow = 0.0;                             /* then set it to zero      */
for (p = bodytab; p < bodytab+nbody; p++)   /* loop over new bodies     */
Type(p) = BODY;                         /* initialize type field    */
}

/*
* STARTOUTPUT: begin output to log file.
*/

void startoutput(void)
{
#if defined(USEFREQ)
printf("\n%8s%10s%10s", "nbody", "freq", "eps");
#else
printf("\n%8s%10s%10s", "nbody", "dtime", "eps");
#endif
#if !defined(QUICKSCAN)
printf("%10s", "theta");
#endif
#if defined(USEFREQ)
printf("%8d%10.2f%10.4f", nbody, freq, eps);
#else
printf("%8d%10.5f%10.4f", nbody, dtime, eps);
#endif
#if !defined(QUICKSCAN)
printf("%10.2f", theta);
#endif
#if defined(USEFREQ)
printf("%10s%10.2f%10.4f\n", usequad ? "true" : "false", freqout, tstop);
#else
printf("%10s%10.5f%10.4f\n", usequad ? "true" : "false", dtout, tstop);
#endif
if (! strnull(options))                     /* print options, if any    */
printf("\n\toptions: %s\n", options);
if (! strnull(savefile))                    /* was state file given?    */
savestate(savefile);                    /* save initial data        */
}

/*
* FORCEREPORT: print staristics on tree construction and force calculation.
*/

void forcereport(void)
{
printf("\n\t%8s%8s%8s%8s%10s%10s%8s\n",
"rsize", "tdepth", "ftree",
"actmax", "nbbtot", "nbctot", "CPUfc");
printf("\t%8.1f%8d%8.3f%8d%10d%10d%8.3f\n",
rsize, tdepth, (nbody + ncell - 1) / ((real) ncell),
actmax, nbbcalc, nbccalc, cpuforce);
}

/*
* OUTPUT: compute diagnostics and output body data.
*/

void output(void)
{
real cmabs, amabs, teff;

diagnostics();                              /* compute std diagnostics  */
ABSV(cmabs, cmvel);                         /* find magnitude of cm vel */
ABSV(amabs, amvec);                         /* find magnitude of J vect */
printf("\n    %8s%8s%8s%8s%8s%8s%8s%8s\n",
"time", "|T+U|", "T", "-U", "-T/U", "|Vcom|", "|Jtot|", "CPUtot");
printf("    %8.3f%8.5f%8.5f%8.5f%8.5f%8.5f%8.5f%8.3f\n",
tnow, ABS(etot[0]), etot[1], -etot[2], -etot[1]/etot[2],
cmabs, amabs, cputime());
#if defined(USEFREQ)
teff = tnow + (freq > 0 ? 0.125/freq : 0);  /* anticipate slightly...   */
#else
teff = tnow + dtime/8;                      /* anticipate slightly...   */
#endif
if (! strnull(outfile) && teff >= tout)     /* time for data output?    */
outputdata();
if (! strnull(savefile))                    /* was state file given?    */
savestate(savefile);                    /* save data for restart    */
}

/*
* OUTPUTDATA: output body data.
*/

void outputdata(void)
{
char namebuf[256];
struct stat buf;
stream outstr;
bodyptr p;

sprintf(namebuf, outfile, nstep);           /* construct output name    */
if (stat(namebuf, &buf) != 0)               /* no output file exists?   */
outstr = stropen(namebuf, "w");         /* create & open for output */
else                                        /* else file already exists */
outstr = stropen(namebuf, "a");         /* reopen and append output */
out_int(outstr, nbody);                     /* write number of bodies   */
out_int(outstr, NDIM);                      /* number of dimensions     */
out_real(outstr, tnow);                     /* and current time value   */
for (p = bodytab; p < bodytab+nbody; p++)   /* loop over all bodies     */
out_real(outstr, Mass(p));              /* output mass of each      */
for (p = bodytab; p < bodytab+nbody; p++)
out_vector(outstr, Pos(p));             /* output positions         */
for (p = bodytab; p < bodytab+nbody; p++)
out_vector(outstr, Vel(p));             /* output velocities        */
if (scanopt(options, "out-phi"))            /* potentials requested?    */
for (p = bodytab; p < bodytab+nbody; p++)
out_real(outstr, Phi(p));           /* output potentials        */
if (scanopt(options, "out-acc"))            /* accelerations requested? */
for (p = bodytab; p < bodytab+nbody; p++)
out_vector(outstr, Acc(p));         /* output accelerations     */
fclose(outstr);                             /* close up output file     */
printf("\n\tdata output to file %s at time %f\n", namebuf, tnow);
#if defined(USEFREQ)
tout += 1.0 / freqout;                      /* schedule next output     */
#else
tout += dtout;                              /* schedule next output     */
#endif
}

/*
* DIAGNOSTICS: compute set of dynamical diagnostics.
*/

local void diagnostics(void)
{
register bodyptr p;
real velsq;
vector tmpv;
matrix tmpt;

mtot = 0.0;                                 /* zero total mass          */
etot[1] = etot[2] = 0.0;                    /* zero total KE and PE     */
CLRM(keten);                                /* zero ke tensor           */
CLRM(peten);                                /* zero pe tensor           */
CLRV(amvec);                                /* zero am vector           */
CLRV(cmpos);                                /* zero c. of m. position   */
CLRV(cmvel);                                /* zero c. of m. velocity   */
for (p = bodytab; p < bodytab+nbody; p++) { /* loop over all particles  */
mtot += Mass(p);                        /* sum particle masses      */
DOTVP(velsq, Vel(p), Vel(p));           /* square vel vector        */
etot[1] += 0.5 * Mass(p) * velsq;       /* sum current KE           */
etot[2] += 0.5 * Mass(p) * Phi(p);      /* and current PE           */
MULVS(tmpv, Vel(p), 0.5 * Mass(p));     /* sum 0.5 m v_i v_j        */
OUTVP(tmpt, tmpv, Vel(p));
MULVS(tmpv, Pos(p), Mass(p));           /* sum m r_i a_j            */
OUTVP(tmpt, tmpv, Acc(p));
CROSSVP(tmpv, Vel(p), Pos(p));          /* sum angular momentum     */
MULVS(tmpv, tmpv, Mass(p));
MULVS(tmpv, Pos(p), Mass(p));           /* sum cm position          */
MULVS(tmpv, Vel(p), Mass(p));           /* sum cm momentum          */
}
etot[0] = etot[1] + etot[2];                /* sum KE and PE            */
DIVVS(cmpos, cmpos, mtot);                  /* normalize cm coords      */
DIVVS(cmvel, cmvel, mtot);
}

/*
* IN_INT, IN_REAL, IN_VECTOR: low level input routines.
*/

local void in_int(stream str, int *iptr)
{
#if !defined(BINARYIO)
if (fscanf(str, "%d", iptr) != 1)
error("in_int: input conversion error\n");
#else
if (fread((void *) iptr, sizeof(int), 1, str) != 1)
#endif
}

local void in_real(stream str, real *rptr)
{
double tmp;

#if !defined(BINARYIO)
if (fscanf(str, "%lf", &tmp) != 1)
error("in_real: input conversion error\n");
*rptr = tmp;
#else
if (fread((void *) rptr, sizeof(real), 1, str) != 1)
#endif
}

local void in_vector(stream str, vector vec)
{
double tmpx, tmpy, tmpz;

#if !defined(BINARYIO)
if (fscanf(str, "%lf%lf%lf", &tmpx, &tmpy, &tmpz) != 3)
error("in_vector: input conversion error\n");
vec[0] = tmpx;
vec[1] = tmpy;
vec[2] = tmpz;
#else
if (fread((void *) vec, sizeof(real), NDIM, str) != NDIM)
#endif
}

/*
* OUT_INT, OUT_REAL, OUT_VECTOR: low level output routines.
*/

#define IFMT  " %d"                             /* output format for ints   */
#define RFMT  " %14.7E"                         /* output format for reals  */

local void out_int(stream str, int ival)
{
#if !defined(BINARYIO)
if (fprintf(str, IFMT "\n", ival) < 0)
error("out_int: fprintf failed\n");
#else
if (fwrite((void *) &ival, sizeof(int), 1, str) != 1)
error("out_int: fwrite failed\n");
#endif
}

local void out_real(stream str, real rval)
{
#if !defined(BINARYIO)
if (fprintf(str, RFMT "\n", rval) < 0)
error("out_real: fprintf failed\n");
#else
if (fwrite((void *) &rval, sizeof(real), 1, str) != 1)
error("out_real: fwrite failed\n");
#endif
}

local void out_vector(stream str, vector vec)
{
#if !defined(BINARYIO)
if (fprintf(str, RFMT RFMT RFMT "\n", vec[0], vec[1], vec[2]) < 0)
error("out_vector: fprintf failed\n");
#else
if (fwrite((void *) vec, sizeof(real), NDIM, str) != NDIM)
error("out_vector: fwrite failed\n");
#endif
}

/*
* SAVESTATE: write current state to disk file.
*/

#define safewrite(ptr,len,str)                  \
if (fwrite((void *) ptr, len, 1, str) != 1) \
error("savestate: fwrite failed\n")

void savestate(string pattern)
{
char namebuf[256];
stream str;
int nchars;

sprintf(namebuf, pattern, nstep & 1);       /* construct alternate name */
str = stropen(namebuf, "w!");
nchars = strlen(getargv0()) + 1;
safewrite(&nchars, sizeof(int), str);
safewrite(getargv0(), nchars * sizeof(char), str);
nchars = strlen(getversion()) + 1;
safewrite(&nchars, sizeof(int), str);
safewrite(getversion(), nchars * sizeof(char), str);
#if defined(USEFREQ)
safewrite(&freq, sizeof(real), str);
#else
safewrite(&dtime, sizeof(real), str);
#endif
#if !defined(QUICKSCAN)
safewrite(&theta, sizeof(real), str);
#endif
safewrite(&eps, sizeof(real), str);
nchars = strlen(options) + 1;
safewrite(&nchars, sizeof(int), str);
safewrite(options, nchars * sizeof(char), str);
safewrite(&tstop, sizeof(real), str);
#if defined(USEFREQ)
safewrite(&freqout, sizeof(real), str);
#else
safewrite(&dtout, sizeof(real), str);
#endif
safewrite(&tnow, sizeof(real), str);
safewrite(&tout, sizeof(real), str);
safewrite(&nstep, sizeof(int), str);
safewrite(&rsize, sizeof(real), str);
safewrite(&nbody, sizeof(int), str);
safewrite(bodytab, nbody * sizeof(body), str);
fclose(str);
}

/*
* RESTORESTATE: restore state from disk file.
*/

if (fread((void *) ptr, len, 1, str) != 1) \

void restorestate(string file)
{
stream str;
int nchars;
string program, version;

str = stropen(file, "r");
program = (string) allocate(nchars * sizeof(char));
version = (string) allocate(nchars * sizeof(char));
if (! streq(program, getargv0()) ||         /* check program, version   */
! streq(version, getversion()))
printf("warning: state file may be outdated\n\n");
#if defined(USEFREQ)
#else
#endif
#if !defined(QUICKSCAN)
#endif
options = (string) allocate(nchars * sizeof(char));
#if defined(USEFREQ)