#include <stdlib.h>
#include <string.h>
#include <stdio.h>

#include "paneltree.h"

#include "page.h"
#include "tree.h"
#include "traverse.h"
#include "keycvt.h"
#include "key.h"

panel_t *parray;

/* thetree - holds all the internal state of the tree representing the
   sources and multipoles. */
static tree_t thetree;

/* debug - print a lot of diagnostic info. */
static int debug = 0;
static FILE *dbgfp;

/* ninter - how many interactions computed? */
static int ninter;

/* pgsys - holds all internal state for the paging system. */
PgSys_t pgsys;

/* UnverseBbox - a bounding box large enough to enclose the centers of
   all the sources.  This is certainly "known" to upper levels of code,
   and doesn't change during the course of a simulation, but for the purposes
   of getting something going, it's easier to just recompute it... */
static tbbox UniverseBbox;

/* A 'Barnes-Hut' style theta parameter */
/* Value below is John Salmon initial suggestion. For no multipole, use: #define THETA 0.0 */
#define THETA 0.5
static double theta2 = (THETA*THETA);

/* cmpPanelKeys - A function that compares two panels based on their
   keys.  Passed to qsort to sort the array of panels */
static int
cmpPanelKeys(const void *v1, const void *v2){
  const panel_t *p1 = (const panel_t *)v1;
  const panel_t *p2 = (const panel_t *)v2;
  return KeyCmp(p1->key, p2->key);
}

/* PanelBB - Get the bounding box associated with a panel.  This just
   lets us use some convenient libsw routines to manipulate bounding-boxes.
   They do nothing profound. */
static void
PanelBB(tbbox *bbp, const panel_t *pp){
  /* Violate the bbox abstraction and just set the fields directly... */
  bbp->ndim = NDIM;
  bbp->rmin[0] = pp->xc - pp->dx*0.5;
  bbp->rmax[0] = pp->xc + pp->dx*0.5;
  bbp->rmin[1] = pp->zc - pp->dz*0.5;
  bbp->rmax[1] = pp->zc + pp->dz*0.5;
}

/* mkAccumulator - Assemble an aggregate (aka internal, aka multipole) 
   cell from several terminal and non-terminal daughters.  In this
   particular code all the data structures are identical, so the
   distinction between internal and terminal daughters is a distinction
   without a difference.  In a more general code, the 'terminal' daughters
   are of type 'tbody', while the 'internal' daughters are of type
   '_cellDataIntermediate'.  Yet another function converts
   _cellDataIntermediate into _cellDataFinal (the finalcontents of the
   tree), but memcpy will serve that purpose here. */
static size_t
mkAccumulator(Key_t k, int nterminal, const panel_t **daughter_panels,
		int ninternal, const panel_t **daughter_cells,
		panel_t *to){
  int i;
  double area, tot_area;
  double vel_area;
  const panel_t *dp;
  tbbox dbb, tobb;
  double ctr[NDIM], sz[NDIM];
  double coverage;

  tot_area = 0.;
  vel_area = 0.;
  MakeBbox(&tobb, NDIM);
  for(i=0; i<nterminal; i++){
    dp = daughter_panels[i];
    area = dp->dx * dp->dz;
    tot_area += area;
    vel_area += dp->vel * area;
    PanelBB(&dbb, dp);
    UnionBbox(&tobb, &dbb, &tobb);
  }

  for(i=0; i<ninternal; i++){
    dp = daughter_cells[i];
    area = dp->dx * dp->dz;
    tot_area += area;
    vel_area += dp->vel * area;
    PanelBB(&dbb, dp);
    UnionBbox(&tobb, &dbb, &tobb);
  }

  to->key = k;
  CenterBbox(&tobb, ctr, sz);
  to->xc = ctr[0];
  to->zc = ctr[1];
  to->dx = sz[0];
  to->dz = sz[1];
  /* Sanity check.  The area-weighted average velocity only makes sense
     if the daughters really "cover" the whole area associated with the
     cell.  There should be no gaps.  If there are, the whole approximation
     is questionable and some investigation is DEFINITELY required. */
  to->vel = vel_area/tot_area;
  coverage = tot_area / (to->dx*to->dz);
  if( coverage > 1.001 || coverage < 0.999 ){
    /* The right thing to do here is probably to mark this cell as
       'suspicious', and to propagate that up the tree as we assemble
       its parents.  During tree traversal we could refuse to do any
       interactions with 'suspicious' cells.  But this would force us
       to distinguish between panels (panel_t) and internal tree cells
       (which would carry along this 'suspicious' field), and that's
       an exercise that I'd rather save for another day... */
    if( debug ){
      fprintf(stderr, "Warning!  coverage=%g  Investigate!\n", coverage);
      fprintf(stderr, "center=(%g,%g), sz=(%g,%g), Key=%s\n",
	      to->xc, to->zc, to->dx, to->dz, PrintKeyHex(k));
    }
  }
  return sizeof(panel_t);
}

/* buildtree -- called from FORTRAN.  Actually builds a tree.  The FORTRAN
   data is copied and left untouched.   On return, a static tree (thetree) is
   available for use by other routines in this file.
*/
void
FORTIFY(buildtree)(int *np, double *xc, double *zc, 
		   double *dx, double *dz, double *vel){
  int n = *np;
  int i;
  static int initized = 0;

  if( !initized ){
    int zero;
    initized = 1;
    zero = 0;
    MPMY_Init(&zero, NULL);
    parray = malloc(n*sizeof(panel_t));
    PgInit(&pgsys, 16, 1000, NULL, 4321);
  }
  dbgfp = fopen("dbg", "a");

  /* First put everything into a data structure that's easy to work
     with in C.  This isn't really necessary, but it makes the code
     easier to write (and maybe read). */
  MakeBbox(&UniverseBbox, NDIM);
  for(i=0; i<n; i++){
    tbbox pbb;
    parray[i].xc = xc[i];
    parray[i].zc = zc[i];
    parray[i].dx = dx[i];
    parray[i].dz = dz[i];
    parray[i].vel = vel[i];
    /* Make sure the Universe contains the new panel */
    PanelBB(&pbb, &parray[i]);
    UnionBbox(&UniverseBbox, &pbb, &UniverseBbox);
  }

  for(i=0; i<n; i++){
    parray[i].key = KeyFromDoubles(&parray[i].xc, &UniverseBbox, 
				   BodyBits(NDIM), PH_ORDER);
  }
				   
  /* Now sort on the key. */
  qsort(parray, n, sizeof(panel_t), cmpPanelKeys);

  TreeBuildStart(&thetree, 2, 1, &pgsys, sizeof(panel_t), 
		 mkAccumulator, (makeFinalCell_t)memcpy);
  for(i=0; i<n; i++){
    TreeBuildAdd(&thetree, parray[i].key, &parray[i], sizeof(panel_t));
  }
  TreeBuildFinish(&thetree);

  ninter = 0;
  if( debug ){
    fprintf(dbgfp, "Tree built:\n");
    /* prkPrintTree(&thetree, 1, printf); */
  }
}

/* Need an error-bound?  And a way of combining cells:
   combine velocities with an area-weighted average.
   max/min the bounds and compute a new center. */

static double obsx, obsz;
static double stress;

static void
gemDoLeaves(const sib_group *sg){
  const panel_t *terminals = Tbodies(sg);
  int i;

  for(i=0; i<sg->ntbody; i++){
    const panel_t *panelp = &terminals[i];
    double srcw2, srch2;
    double srcz, diffx;
    double e12;
    srcw2 = panelp->dx*0.5;
    srch2 = panelp->dz*0.5;
    srcz = panelp->zc;
    diffx = obsx - panelp->xc;
    FORTIFY(strn)(&srcw2, &srch2, &srcz, &diffx, &obsz, &e12);
    stress += panelp->vel * e12;
  }
  ninter += sg->ntbody;
}
    
static int
gemDoCell(Key_t key, const tcell *tc){
  const panel_t *panelp = Contents(tc);
  double srcw2, srch2;
  double srcz, diffx, e12;
  double sz2, dist2;

  srcw2 = panelp->dx*0.5;
  srch2 = panelp->dz*0.5;
  srcz = panelp->zc;
  diffx = obsx - panelp->xc;
  /* Is this panel far enough, weak enough, or smooth enough that we
     can do the interaction? */
  sz2 = srcw2*srcw2 + srch2*srch2;
  dist2 = diffx*diffx + (obsz-srcz)*(obsz-srcz);
  if( debug )
    fprintf(dbgfp, "%s %g %g %g %g %g\n", PrintKey(key), sz2, obsx, obsz, panelp->xc, srcz);
  if( sz2 < theta2*dist2 ){
    FORTIFY(strn)(&srcw2, &srch2, &srcz, &diffx, &obsz, &e12);
    stress += panelp->vel * e12;
    ninter++;
    return MAC_ACCEPT;
  }else{
    return MAC_REJECT;
  }
}

/* Greens_function:  call back to STRN */
void
FORTIFY(sumtree)(double *xcp, double *zcp, double *strsp){
  obsx = *xcp;
  obsz = *zcp;
  stress = 0.;

  TraverseSg(&thetree, TreeRoot(&thetree), gemDoCell, gemDoLeaves);
  *strsp += stress;
  if( debug )
    fprintf(dbgfp, "*strsp += %g = %g\n", stress, *strsp);
}

void
FORTIFY(forgettree)(void){
  PagePollTillDone(thetree.pgsys);
  TreeFree(&thetree);
  fprintf(dbgfp, "freetree: ninter=%d\n", ninter);
  fclose(dbgfp);
}