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


/* **********************************************************************
 * entity Cell{fParent, fBound}
 *   quad-tree tree node
 * **********************************************************************/
entity Cell{fParent, fQuadId}
  ename fParent;       	/* parent ename */
  int fQuadId;         	/* quadrant id */
{
  /* private states */
  int fNParticle;       /* number of particles in the cell */
  float fMass;		/* total mass in the cell */
  float fBound[eNBound];/* spatial boundary of the cell */
  SubCell fQ[eNQ];	/* sub cells, or sub-branches of the tree */
  Vector fCMass;	/* center of the mass in the cell */

  /* message types */
  message Boot {
    float bound[eNBound];	/* boundary of the cell */
    int nParticle;	/* number of particles */
  } mBoot;		/* boot message from the parent */
  message AddP {
    Particle p;		/* particle */
  } mAddP;		/* add a particle to this cell */
  message TreeDone{
    int quadId;		/* quadrant of the sub-tree */
    float mass;		/* total mass of the sub-tree */
    Vector cmass;	/* center of mass from the sub-tree */
  } mTreeDone;		/* return message from sub-tree construction */
  message ComputeAcc {
    Particle p;		/* particle */
  } mComputeAcc;	/* for computation message to the sub-tree */
  message ComputeDone {
    Vector accOnP;	/* acceleration on the particle */
    int nCompute;	/* number of computation */
  } mComputeDone;	/* return message from force computation */
  message ShutDown {	
  } mShutDown;		/* shutdown message */

  /* temporary variables */
  int boundId;          /* boundary id */
  int quadId;		/* quadrant id */
  int partId;		/* particle id */
  int nPartQ[eNQ];      /* num of particles in a quadrant */
  int nCompute;		/* number of for computations */
  float D;		/* threshold dimension of the bounding box */   
  float r;		/* distance between the particle and the center
				of mass */
  float midX, midY;	/* mid points of current cell bounding box */
  float rangeX, rangeY; /* ranges of current cell bounding box */
  float boundQ[eNQ][eNBound]; 	/* bounding boxes for each sub-quadrant */
  Bool treeDoneFlag = eFalse;	/* to see if the sub-trees are constructed */
  Bool treeDone[eNQ];		/* to see if each sub-tree is constructed */
  Vector accOnP;		/* acceleration vector */

  wait until mBoot = mtype(Boot);

  /* initialization */
  fNParticle = mBoot.nParticle;
  fMass = fCMass.x = fCMass.y = accOnP.x = accOnP.y = 0;

  /* check the number of particles in the simulation */
  switch (fNParticle) {
    case 0: 	/* done */
      invoke fParent with TreeDone{fQuadId, fMass, fCMass::sizeof(Vector)};
      break;
    case 1:	/* initialize a current cell, done */
      wait until mAddP = mtype(AddP);
      fMass = mAddP.p.mass;
      fCMass.x = mAddP.p.info[ePos].x;
      fCMass.y = mAddP.p.info[ePos].y;
      invoke fParent with TreeDone {fQuadId, fMass, fCMass::sizeof(Vector)};
      break;
    default:	/* grow the quad-tree */
      /* initailize the bounding conditions for the sub-quadrants */
      for (quadId = 0; quadId < eNQ; quadId++) {
	nPartQ[quadId] = 0;
	treeDone[quadId] = eFalse;
      }

      for (boundId = 0; boundId < eNBound; boundId++) {
	fBound[boundId] = mBoot.bound[boundId];
      }

      midX = (fBound[eMinX] + fBound[eMaxX])/2;
      midY = (fBound[eMinY] + fBound[eMaxY])/2;

      boundQ[eQ2][eMinX] = boundQ[eQ3][eMinX] = fBound[eMinX];
      boundQ[eQ1][eMaxX] = boundQ[eQ4][eMaxX] = fBound[eMaxX];
      boundQ[eQ3][eMinY] = boundQ[eQ4][eMinY] = fBound[eMinY];
      boundQ[eQ1][eMaxY] = boundQ[eQ2][eMaxY] = fBound[eMaxY];
      
      boundQ[eQ1][eMinX] = boundQ[eQ2][eMaxX] 
	= boundQ[eQ3][eMaxX] = boundQ[eQ4][eMinX] = midX;
      boundQ[eQ1][eMinY] = boundQ[eQ2][eMinY] 
	= boundQ[eQ3][eMaxY] = boundQ[eQ4][eMaxY] = midY;

      for (partId = 0; partId < fNParticle; partId++) {
	wait until mAddP = mtype(AddP);

	/* check to see if the particle has flied out of the boundaries */
	if ((mAddP.p.info[ePos].x < fBound[eMinX]) ||
	    (mAddP.p.info[ePos].x > fBound[eMaxX]) ||
	    (mAddP.p.info[ePos].x < fBound[eMinY]) ||
	    (mAddP.p.info[ePos].x > fBound[eMaxY])) {
	  /* if so, randomly re-initialize the particle */
	  RandInitParticle(&mAddP.p, mAddP.p.id,
			   fBound[eMinX], fBound[eMaxX], 
			   fBound[eMinY], fBound[eMaxY]);
	}

	/* assign the particle to sub-quadrants */
	if (mAddP.p.info[ePos].x > midX) {
	  if (mAddP.p.info[ePos].y > midY) { /* quardrant 1 (+,+) */
	    quadId = eQ1;
	  } else {	        	     /* quardrant 4 (+,-) */
	    quadId = eQ4;
	  }
	} else {
	  if (mAddP.p.info[ePos].y > midY) { /* quardrant 2 (-,+) */
	    quadId = eQ2;
	  } else {		     	     /* quardrant 3 (-,-) */
	    quadId = eQ3;
	  }
	}    

	if (!nPartQ[quadId]) {
	  fQ[quadId].cell = new Cell{self, quadId};
	} 
	nPartQ[quadId]++;
	invoke fQ[quadId].cell with AddP{mAddP.p::sizeof(Particle)};
      }
      
      for (quadId = 0; quadId < eNQ; quadId++) {
	if (nPartQ[quadId]) {
	  invoke fQ[quadId].cell with Boot {
	    &boundQ[quadId][0]::eNBound*sizeof(float), 
	    nPartQ[quadId]
	  };
	} else {
	  treeDone[quadId] = eTrue;
	}
      }
      
      /* check if the sub-tree is completed */
      for (treeDoneFlag = eFalse; !treeDoneFlag;) {
	wait until mTreeDone = mtype(TreeDone) {
	  treeDone[mTreeDone.quadId] = eTrue;
	  fQ[mTreeDone.quadId].mass = mTreeDone.mass;
	  fQ[mTreeDone.quadId].cmass.x = mTreeDone.cmass.x;
	  fQ[mTreeDone.quadId].cmass.y = mTreeDone.cmass.y;
	  fMass += mTreeDone.mass;
	  fCMass.x += mTreeDone.mass*mTreeDone.cmass.x;
	  fCMass.y += mTreeDone.mass*mTreeDone.cmass.y;
	}
	for (treeDoneFlag = eTrue, quadId = 0; quadId < eNQ; quadId++) {
	  treeDoneFlag &= treeDone[quadId];
	}
      }
      /* compute the center of mass for the current cell */
      fCMass.x /= fMass;
      fCMass.y /= fMass;
      invoke fParent with TreeDone{fQuadId, fMass, fCMass::sizeof(Vector)};
      break;
  }

  /* wait for the force computation message, or the shutdown message */
  while (eTrue) {
    wait until {
      mComputeAcc = mtype(ComputeAcc) {
	nCompute = 0;
	if (fNParticle == 1) {
	  /* if the cell contains only 1 particle, compute the force
		between the target particle with the local particle */
	  r = ComputeDist(&mComputeAcc.p.info[ePos], &fCMass);
	  ComputeGAcc(&fCMass, &mComputeAcc.p.info[ePos], 
		      &accOnP, fMass, r);
	  nCompute++;
	} else {
	  /* if the cell contains more than 1 particle, compute the
		force based on the center of mass if the D/r ratio is 
		less than the threshhold */
	  for (quadId = 0; quadId < eNQ; quadId++) {
	    if (nPartQ[quadId]) {
	      rangeX = boundQ[quadId][eMaxX] - boundQ[quadId][eMinX];
	      rangeY = boundQ[quadId][eMaxY] - boundQ[quadId][eMinY];
	      if (rangeX > rangeY) {
		D = rangeX;
	      } else {
		D = rangeY;
	      }
	      r = ComputeDist(&mComputeAcc.p.info[ePos],
			      &fQ[quadId].cmass);	      
	      if (r > 0) {
		if (D/r < cTheta) {
		  ComputeGAcc(&fQ[quadId].cmass, &mComputeAcc.p.info[ePos], 
			      &accOnP, fMass, r);
		  nCompute++;
		} else {
		  invoke fQ[quadId].cell with ComputeAcc {
		    mComputeAcc.p::sizeof(Particle)
		  };
		  wait until mComputeDone = mtype(ComputeDone);
		  accOnP.x += mComputeDone.accOnP.x;
		  accOnP.y += mComputeDone.accOnP.y;
		  nCompute += mComputeDone.nCompute;
		}
	      }
	    }
	  }
	}
	/* computation done */
	invoke fParent with ComputeDone {
	  accOnP::sizeof(Vector), nCompute
	};
      } or mShutDown = mtype(ShutDown) {
	for (quadId = 0; quadId < eNQ; quadId++) {
	  if (nPartQ[quadId]) {
	    invoke fQ[quadId].cell with ShutDown;
	  }
	}
	break;
      }
    }
  }
}