#include "/u/data05/idpbbsp/maisie/bin/mayc.h"
#include "Cell.h"
#include "Main.h"
#include "MaisieUtil.h"
#include "Para.h"
#include "Particle.h"


/* **********************************************************************
 * entity driver{argc, argv}
 *   n-body simulation driver
 * **********************************************************************/
entity driver{argc, argv}
  int argc;
  char **argv;
{ 
  /* private states */
  int fNStep;		/* number of steps for sim */
  int fNCompute;	/* number of force computation during the sim */
  int fNParticle;	/* number of particles in the sim */
  float fBound[eNBound];	/* spatial specification */
  ename fQRoot;			/* root of the quadtree */
  Particle *fPList = NULL;	/* particle list for all the particles */

  /* messages */
  message TreeDone{
    int quadId;		/* quadrant id */
    float mass;		/* total mass of this tree branch upto this node */
    Vector cmass;	/* center of mass of this tree node */
  } mTreeDone;		/* return message from tree construction */
  message UpdateP { 
    Particle p; 	/* particle */
  } mUpdateP;		/* message to update the particle */
  message ComputeDone {
    Vector accOnP;	/* acceleration on the particle */
    int nCompute;	/* number of computations */
  } mComputeDone;	/* return message from the force computation */

  /* private functions */
  void InitSim();	/* initialize the simulation */
  void PrintSimConfig();/* print sim configuration */
  void Update();	/* update a particle */

  /* temporary variables */
  int stepI;		/* step index */
  int partId;		/* particle id */
  Vector *accOnP;	/* acceleration on a particle */

  /* input checking */
  if (argc != 7) {
    Error("Usage:  NBody <MinX> <MaxX> <MinY> <MaxY> <N Bodies> <N Steps>");
  }

  /* initialization */
  InitSim();
  PrintSimConfig();

  /* start simulation */
  for (stepI = 0; stepI < fNStep; stepI++) {
    printf("Step:  %d\n", stepI);

    /* build tree */
    fQRoot = new Cell{self, -1};
    for (partId = 0; partId < fNParticle; partId++) {
      PrintParticle(&fPList[partId]);
      invoke fQRoot with AddP{fPList[partId]::sizeof(Particle)};
    }
    invoke fQRoot with Boot{&fBound[0]::eNBound*sizeof(float), fNParticle};
    wait until mtype(TreeDone);

    /* compute acceleration */
    for (partId = 0; partId < fNParticle; partId++) {
      invoke fQRoot with ComputeAcc{fPList[partId]::sizeof(Particle)};
      wait until mComputeDone = mtype(ComputeDone);
      accOnP[partId].x = mComputeDone.accOnP.x;
      accOnP[partId].y = mComputeDone.accOnP.y;
      fNCompute += mComputeDone.nCompute;
    }
    printf("Theta %f, NParticle %d, nCompute:  %d\n", 
	   cTheta, fNParticle, fNCompute);

    /* update */
    for (partId = 0; partId < fNParticle; partId++) {
      Update(partId);
    }

    /* shutdown */
    invoke fQRoot with ShutDown;
  }
  terminate();
}


/* **********************************************************************
 * void driver::InitSim
 *   initializes the privates states from the program input arguments
 * **********************************************************************/
void
driver::InitSim() 
{
  /* initialize the spatial dimension */
  fBound[eMinX] = atof(argv[1]);
  fBound[eMaxX] = atof(argv[2]);
  fBound[eMinY] = atof(argv[3]);
  fBound[eMaxY] = atof(argv[4]);
  fNParticle = (int) atoi(argv[5]);
  fNStep = (int) atoi(argv[6]);
  fNCompute = 0;

  if ((fBound[eMaxX] <= fBound[eMinX]) || 
      (fBound[eMaxY] <= fBound[eMinY]) || 
      (fNParticle <= 0) || 
      (fNStep <= 0)) {
    Error("Input Out of Range");
  }

  /* randomly initializes the particle configurations */
  srand();
  fPList = NewParticle(fNParticle);
  accOnP = (Vector *) Malloc(sizeof(Vector)*fNParticle);
  for (partId = 0; partId < fNParticle; partId++) {
    RandInitParticle(&fPList[partId], partId, 
		     fBound[eMinX], fBound[eMaxX], 
		     fBound[eMinY], fBound[eMaxY]);
  }
}


/* **********************************************************************
 * void driver::PrintSimConfig()
 *   print the simulation configuration 
 * **********************************************************************/
void
driver::PrintSimConfig() 
{
  printf("X:  (%f, %f)\n", fBound[eMinX], fBound[eMaxX]);
  printf("Y:  (%f, %f)\n", fBound[eMinY], fBound[eMaxY]);
  printf("Number of Particiles:  %d\n", fNParticle);
  printf("Simulation Steps:  %d\n", fNStep);
}


/* **********************************************************************
 * void driver::Update
 *   update a single particle with new acceleration, velocity, and 
 *   position
 * **********************************************************************/
void
driver::Update(partId)
  int partId;
{
  fPList[partId].info[eAcc].x = accOnP[partId].x;
  fPList[partId].info[eAcc].y = accOnP[partId].y;

  fPList[partId].info[eVel].x += fPList[partId].info[eAcc].x;
  fPList[partId].info[eVel].y += fPList[partId].info[eAcc].y;

  fPList[partId].info[ePos].x
    += (fPList[partId].info[eVel].x + (1/2)*fPList[partId].info[eAcc].x);
  fPList[partId].info[ePos].y
    += (fPList[partId].info[eVel].y + (1/2)*fPList[partId].info[eAcc].y);
}