/****************************************************************************/
/* LOAD.C: routines to create tree.  Public routines: maketree().           */
/*                                                                          */
/* Copyright (c) 1993 by Joshua E. Barnes, Honolulu, HI.                    */
/* It's free because it's yours.                                            */
/****************************************************************************/
/* adapted by Bernd Vollmer, CDS, Observatoire de Strasbourg 2004           */
/****************************************************************************/
#include "nrutil.h"
#include "defs.h"
 
local void newtree(void);			/* flush existing tree      */
local cellptr makecell(void);			/* create an empty cell     */
local void expandbox(bodyptr, int);		/* set size of root cell    */
local void loadbody(bodyptr);			/* load body into tree      */
local int subindex(bodyptr, cellptr);		/* compute subcell index    */
local void hackcofm(cellptr, real);		/* find centers of mass     */
local void setrcrit(cellptr, vector, real);	/* set cell's crit. radius  */
local void threadtree(nodeptr, nodeptr);	/* set next and more links  */
			/* compute quad moments     */

/*
 * MAKETREE: initialize tree structure for hierarchical force calculation
 * from body array btab, which contains nbody bodies.
 */
 
local bool bh86, sw93;				/* use alternate criteria   */
local nodeptr pex[100000];
 
void maketree(bodyptr btab, int nbody)
{
    bodyptr p;
 
    newtree();                                  /* flush existing tree, etc */
    root = makecell();				/* allocate the root cell   */
    SETV(Pos(root),offset);		        /* initialize the midpoint  */
    expandbox(btab, nbody);                     /* and expand cell to fit   */
    maxlevel = 0;                               /* init count of levels     */
    for (p = btab; p < btab+nbody; p++)         /* loop over bodies...      */
        if (Mass(p) != 0.0)                     /*   exclude test particles */
            loadbody(p);                        /*     and insert into tree */
    bh86 = scanopt(options, "bh86");		/* set flags for alternate  */
    sw93 = scanopt(options, "sw93");		/* ...cell opening criteria */
    if (bh86 && sw93)				/* can't have both at once  */
	error("maketree: options bh86 and sw93 are incompatible\n");
    hackcofm(root, rsize);                      /*  find c-of-m coordinates  */
    threadtree((nodeptr) root, NULL);           /* add Next and More links  */
}

/*
 * NEWTREE: reclaim cells in tree, prepare to build new one.
 */
 
local nodeptr freecell = NULL;          	/* list of free cells       */
 
local void newtree(void)
{
    permanent bool firstcall = TRUE;
    nodeptr p;
 
    if (! firstcall) {                          /* tree data to reclaim?    */
        p = (nodeptr) root;                     /*   start with the root    */
        while (p != NULL)                       /*   loop scanning tree     */
            if (Type(p) == CELL) {              /*     found cell to free?  */
                Next(p) = freecell;             /*       link to front of   */
                freecell = p;                   /*       ...existing list   */
                p = More(p);                    /*       scan down tree     */
            } else                              /*     skip over bodies     */
                p = Next(p);                    /*       go on to next      */
    } else                                      /* first time through       */
        firstcall = FALSE;                      /*   so just note it        */
    root = NULL;                                /* flush existing tree      */
    cellused = 0;                               /* reset cell count         */
}
 
/*
 * MAKECELL: return pointer to free cell.
 */
 
local cellptr makecell(void)
{
    cellptr c;
    int i;
 
    if (freecell == NULL)                       /* no free cells left?      */
        c = (cellptr) allocate(sizeof(cell));   /*   allocate a new one     */
    else {                                      /* use existing free cell   */
        c = (cellptr) freecell;                 /*   take one on front      */
        freecell = Next(c);                     /*   go on to next one      */
    }
    Type(c) = CELL;                             /* initialize cell type     */
    for (i = 0; i < NSUB; i++)                  /* loop over subcells       */
        Subp(c)[i] = NULL;                      /*   and empty each one     */
    cellused++;                                 /* count one more cell      */
    return (c);
}

/*
 * EXPANDBOX: find range of coordinate values (with respect to root)
 * and expand root cell to fit.  The size is floatd at each step to
 * take advantage of exact representation of powers of two.
 */
 
local void expandbox(bodyptr btab, int nbody)
{
    real xyzmax;
    bodyptr p;
    int k;
 
    xyzmax = 0.0;
    for (p = btab; p < btab+nbody; p++)
	for (k = 0; k < NDIM; k++)
	    xyzmax = MAX(xyzmax, rabs(Pos(p)[k] - Pos(root)[k]));
    while (rsize < 2 * xyzmax)
	rsize = 2 * rsize;
}

/*
 * LOADBODY: descend tree and insert body p in appropriate place.
 */
 
local void loadbody(bodyptr p)
{
    cellptr q, c;
    int qind, lev, k, i;
    real qsize,em1;
 
    q = root;                                   /* start with tree root     */
    qind = subindex(p, q);			/* get index of subcell     */
    qsize = rsize;                              /* keep track of cell size  */
    lev = 0;                                    /* count levels descended   */
    while (Subp(q)[qind] != NULL) {             /* loop descending tree     */
      if (Type(Subp(q)[qind]) == BODY) {      /*   reached a "leaf"?      */
	DISTV(em1,Pos(p),Pos(Subp(q)[qind]));
	if (em1 == 0.0) {
	  if (maketr == 0) {
	    Brothers(Subp(q)[qind])[Nbrothers(Subp(q)[qind])]=(nodeptr) p;
	    Nbrothers(Subp(q)[qind])++;
	    Brothers(p)[Nbrothers(p)]=Subp(q)[qind];
	    Nbrothers(p)++;
	  }
	  bnumber++;	      
	  if (excor > 0) {
	    for (i=0;i<excor-1;i++) {
	      DISTV(em1,Pos(pex[i]),Pos(Subp(q)[qind]));
	      if (em1 == 0.0) {
		if (maketr == 0) {
		  Brothers(pex[i])[Nbrothers(pex[i])]=Subp(q)[qind];
		  Nbrothers(pex[i])++;
		}
		bnumber++;
		cmore++;
	      }
	    }
	  }
	  pex[excor]=Subp(q)[qind];
	  excor++;
	  break;
	}
	c = makecell();                     /*     allocate new cell    */
	for (k = 0; k < NDIM; k++)		/*     initialize midpoint  */
	  Pos(c)[k] = Pos(q)[k] +		/*       offset from parent */
	    (Pos(p)[k] < Pos(q)[k] ? - qsize : qsize) / 4;
	Subp(c)[subindex((bodyptr) Subp(q)[qind], c)] = Subp(q)[qind];
	/*     put body in cell     */
	Subp(q)[qind] = (nodeptr) c;        /*     link cell in tree    */
      }
      q = (cellptr) Subp(q)[qind];		/*   advance to next level  */
      qind = subindex(p, q);			/*   get index to examine   */
      qsize = qsize / 2;                      /*   shrink current cell    */
      lev++;                                  /*   count another level    */
    }
    Subp(q)[qind] = (nodeptr) p;                /* found place, store p     */
    Rcrit2(q)=qsize;
    maxlevel = MAX(maxlevel, lev);		/* remember maximum level   */
}

/*
 * SUBINDEX: compute subcell index for body p in cell q.
 */
 
local int subindex(bodyptr p, cellptr q)
{
    int ind, k;
 
    ind = 0;					/* accumulate subcell index */
    for (k = 0; k < NDIM; k++)			/*   loop over dimensions   */
        if (Pos(q)[k] <= Pos(p)[k])		/*     if beyond midpoint   */
            ind += NSUB >> (k + 1);             /*       skip over subcells */
    return (ind);
}

/*
 * HACKCOFM: descend tree finding center-of-mass coordinates and
 * setting critical cell radii.
 */
 
local void hackcofm(cellptr p, real psize)
{
    vector cmpos, tmpv;
    int i, k;
    nodeptr q;
 
    Mass(p) = 0.0;                              /* init total mass...       */
    CLRV(cmpos);                                /* and center of mass       */
    for (i = 0; i < NSUB; i++)                  /* loop over subnodes       */
	if ((q = Subp(p)[i]) != NULL) {         /*   does subnode exist?    */
	    if (Type(q) == CELL)		/*     and is it a cell?    */
		hackcofm((cellptr) q, psize/2); /*       find subcell cm    */
	    Mass(p) += Mass(q);                 /*     sum total mass       */ 
	    MULVS(tmpv, Pos(q), Mass(q));       /*     weight pos by mass   */
	    ADDV(cmpos, cmpos, tmpv);           /*     sum c-of-m position  */
	}
    DIVVS(cmpos, cmpos, Mass(p));               /* rescale cms position     */
    for (k = 0; k < NDIM; k++)			/* check tree structure...  */
      if (cmpos[k] < Pos(p)[k] - psize/2 ||     /*   if out of bounds       */
	  Pos(p)[k] + psize/2 <= cmpos[k]){	/*   in either direction    */
	if (psize > 1.) printf("hackcofm: tree structure error\n");
      }
    setrcrit(p, cmpos, psize);                  /* set critical radius      */
    /*SETV(Pos(p), cmpos);*/			/*   and center-of-mass pos   */
}

/*
 * SETRCRIT: assign critical radius for cell p, using center-of-mass
 * position cmpos and cell size psize.
 */

local void setrcrit(cellptr p, vector cmpos, real psize)
{
    real rc, bmax2, dmin;
    int k;

    rc=sqrt(2.)*psize+theta/2.;
    Rcrit2(p) = rc*rc;			        /* store square of radius   */
}

/*
 * THREADTREE: do a recursive treewalk starting from node p,
 * with next stop n, installing Next and More links.
 */
 
local void threadtree(nodeptr p, nodeptr n)
{
    int ndesc, i;
    nodeptr desc[NSUB+1];
 
    Next(p) = n;                                /* link to next node        */
    if (Type(p) == CELL) {                      /* any children to thread?  */
        ndesc = 0;                              /*   count extant children  */
        for (i = 0; i < NSUB; i++)              /*   loop over subnodes     */
            if (Subp(p)[i] != NULL)             /*     found a live one?    */
                desc[ndesc++] = Subp(p)[i];     /*       store in table     */
        More(p) = desc[0];                      /*   link to first child    */
        desc[ndesc] = n;                        /*   end table with next    */
        for (i = 0; i < ndesc; i++)             /*   loop over children     */
            threadtree(desc[i], desc[i+1]);     /*     thread each w/ next  */
    }
}