#include <stdio.h>
#include <math.h>
#include <string.h>
#include <stdlib.h>
#include <time.h>
#include "mpi.h"

#define NNMAX 200

main(int argc, char** argv) {
  //init MPI stuff here..
  int my_pe;
  int npe;
  int source;
  int dest;
  int tag = 50;
  char message[100];
  char hostname[1024];
  MPI_Status status;
  
  MPI_Init(&argc, &argv);
  MPI_Comm_rank(MPI_COMM_WORLD, &my_pe);
  MPI_Comm_size(MPI_COMM_WORLD, &npe);
  FILE *fp1,*fpE, *fpt,*fout;
  hostname[1023] = '\0';
  if(my_pe == 0) {
    fout=fopen("host","w");
  }
  if(my_pe != 0) {
    fout=fopen("node","w");
  }
 double year,sm,day,oneau,twoau,kpc,dt,myr,tmax,Mdisc,Mcore;
 double mv,mvc,alpha,rdiskmin,alpha2,g,A;
 double m0,rc,mcdm0,PI,rv,t2,t1;
 double dr,dro,math,rkpc,delvt,mass,vm;
 double xsq,ysq,restart,galaxy,elapsed;
 double delv,rad,rn,rmin,drmin,En,En0,En2,En02,EnT0,EnT;
 double vm2,U,t,fr,ex,ey,f,rsq,vxo,vyo;
 double Fintx,Finty,fint,vxg2,vxg,vyg,vyg2;
 double exi,eyi,Fx,Fy,dx,dy,cdm;
 double Mdisc2, Mcore2, m02, mcdm02, dx2, dy2; 
 double  th;
 double vs[120],rvs[120];
 double vs2[120],rvs2[120];

 int  n0,idump,ik,iv,nvs,starran,flag,idumpmax,inn;
 int  i,j,k,n,ind,ict,ierr[0];
 int i1,i2,i3,i4;
 int nstar_pe;

 double *x,  *y,  *vx,  *vy,  *Ax,  *Ay;
 double *m,  *r,  *mt,  *mcdm,  *mtb;
 double *x2,  *y2,  *vx2,  *vy2,  *Ax2,  *Ay2;
 double *m2,  *r2,  *mt2,  *mcdm2,  *mt2b;
 double *Axg, *Ayg;
 int  *nnA,  *nn1,  *jnA,  *jn1;
 int *nnB,  *nn2,  *jnB,  *jn2;

char fname[20];
time_t timeraw;
       struct tm * timeinfo;
/* File pointers for data output */
 
 int id = 0, id2 = 1;

 /* ============================ START of Program ====================== */

 /* inn=1 is Nearest Neighbor */
 /* inn=0 is N-Body */
inn=1;
 ict = 1;
// WITH CDM (cdm=1)
// cdm=1.0;
 // NO CDM!!
 cdm=1.0;
 
 /* galaxy = 2 is NGC 2403 */
 /* galaxy = 1 is NGC 3198 */
 /* galaxy = 0 is ANDROMEDA */
 galaxy = 0;
 
 /*version info */
 // version = "galaxy_v5.1.c";
 
 /* one solar mass in kg*/
 sm = 1.98892e30;
 /* Length of day in seconds*/
 day = 24.0*3600.0;
 /* Length of a year in seconds*/
 year = 365.25*day;
 /* 1 million years */
 myr=1.0E6*year;
 /* one AU in meters*/
 oneau = 149597871.0E3;
 twoau = 229388960.0E3;
/* one kilo parsec*/
kpc = 2.06e8*oneau;
/*-----------------------------------------*/
dt = year*200.0;
/* Maximum simulation time  (seconds)*/
tmax = myr*12000.0;
// for dt=year*100, 100000/10 is every 1myr
idumpmax=100000/20;
idump = 0;

/* Gravity cuttoff parameters */
rmin=1.0*kpc;
drmin=rmin*rmin;

/*------------------------------*/
/*DISC AND CORE MASSES in kg!!*/
 if(galaxy == 2) {
   Mdisc=1.18e10*sm*.88;
   Mcore=.72e10*sm*0.75;
 }
 if(galaxy == 1) {
   Mdisc=.98e10*sm;
   Mcore=1.03e10*sm;
 }
 if(galaxy == 0) {
   Mdisc=7.0e10*sm;
   Mcore=2.5e10*sm;
 }

 Mdisc2 = 3.528957327e10*sm;
 Mcore2 = 1.18852474e10*sm;
/*------------------------------*/
 n=8000;
 n0=7800;
 nstar_pe = n/npe;
 if (n0 > n ) {
   printf("ERROR: n0 must be less than n\n");
   exit(1);
 }


 /*===========================================*/
 /* ALLOCATE MEMORY FOR ARRAYS */
 /*===========================================*/

  Axg=(double *) malloc(2 * sizeof(double));
  if (Axg== NULL){
    printf("ERROR: Failed to allocate memory \n");
    exit (1);}

  Ayg=(double *) malloc(2 * sizeof(double));
  if (Ayg == NULL){
    printf("ERROR: Failed to allocate memory \n");
    exit (1);}


  /* Galaxy 1 */
  x=(double *) malloc(n * sizeof(double));
  if (x == NULL){
    printf("ERROR: Failed to allocate memory \n");
    exit (1);}

  y=(double *) malloc(n * sizeof(double));
  if (y == NULL){
    printf("ERROR: Failed to allocate memory \n");
    exit (1);}

  vx=(double *) malloc(n * sizeof(double));
  if (vx == NULL){
    printf("ERROR: Failed to allocate memory \n");
    exit (1);}

  vy=(double *) malloc(n * sizeof(double));
  if (vy == NULL){
    printf("ERROR: Failed to allocate memory \n");
    exit (1);}

  Ax=(double *) malloc(n * sizeof(double));
  if (Ax == NULL){
    printf("ERROR: Failed to allocate memory \n");
    exit (1);}

  Ay=(double *) malloc(n * sizeof(double));
  if (Ay == NULL){
    printf("ERROR: Failed to allocate memory \n");
    exit (1);}

  m=(double *) malloc(n * sizeof(double));
  if (m == NULL){
    printf("ERROR: Failed to allocate memory \n");
    exit (1);}

  r=(double *) malloc(n * sizeof(double));
  if (r == NULL){
    printf("ERROR: Failed to allocate memory \n");
    exit (1);}

  mt=(double *) malloc(n * sizeof(double));
  if (mt == NULL){
    printf("ERROR: Failed to allocate memory \n");
    exit (1);}

  mcdm=(double *) malloc(n * sizeof(double));
  if (mcdm == NULL){
    printf("ERROR: Failed to allocate memory \n");
    exit (1);}

  mtb=(double *) malloc(n * sizeof(double));
  if (mtb == NULL){
    printf("ERROR: Failed to allocate memory \n");
    exit (1);}




  /* Galaxy 2 */
  x2=(double *) malloc(n * sizeof(double));
  if (x2 == NULL){
    printf("ERROR: Failed to allocate memory \n");
    exit (1);}

  y2=(double *) malloc(n * sizeof(double));
  if (y2 == NULL){
    printf("ERROR: Failed to allocate memory \n");
    exit (1);}

  vx2=(double *) malloc(n * sizeof(double));
  if (vx2 == NULL){
    printf("ERROR: Failed to allocate memory \n");
    exit (1);}

  vy2=(double *) malloc(n * sizeof(double));
  if (vy2 == NULL){
    printf("ERROR: Failed to allocate memory \n");
    exit (1);}

  Ax2=(double *) malloc(n * sizeof(double));
  if (Ax2 == NULL){
    printf("ERROR: Failed to allocate memory \n");
    exit (1);}

  Ay2=(double *) malloc(n * sizeof(double));
  if (Ay2 == NULL){
    printf("ERROR: Failed to allocate memory \n");
    exit (1);}

  m2=(double *) malloc(n * sizeof(double));
  if (m2 == NULL){
    printf("ERROR: Failed to allocate memory \n");
    exit (1);}

  r2=(double *) malloc(n * sizeof(double));
  if (r2 == NULL){
    printf("ERROR: Failed to allocate memory \n");
    exit (1);}

  mt2=(double *) malloc(n * sizeof(double));
  if (mt2 == NULL){
    printf("ERROR: Failed to allocate memory \n");
    exit (1);}

  mcdm2=(double *) malloc(n * sizeof(double));
  if (mcdm2 == NULL){
    printf("ERROR: Failed to allocate memory \n");
    exit (1);}

  mt2b=(double *) malloc(n * sizeof(double));
  if (mt2b == NULL){
    printf("ERROR: Failed to allocate memory \n");
    exit (1);}


  /*----------------------------------------*/
  /* Pointer Arrays */
  /*----------------------------------------*/
  /* Galaxy 1 */
  nnA=(int *) malloc(n * sizeof(int));
  if (nnA== NULL){
    printf("ERROR: Failed to allocate memory \n");
    exit (1);}

  nn1=(int *) malloc(n * sizeof(int));
  if (nn1== NULL){
    printf("ERROR: Failed to allocate memory \n");
    exit (1);}

  jnA=(int *) malloc(NNMAX*n*sizeof(int));
  if (jnA== NULL){
    printf("ERROR: Failed to allocate memory \n");
    exit (1);}

  jn1=(int *) malloc(NNMAX*n*sizeof(int));
  if (jn1== NULL){
    printf("ERROR: Failed to allocate memory \n");
    exit (1);}

  /* Galaxy 2 */

  nnB=(int *) malloc(n * sizeof(int));
  if (nnB== NULL){
    printf("ERROR: Failed to allocate memory \n");
    exit (1);}

  nn2=(int *) malloc(n * sizeof(int));
  if (nn2== NULL){
    printf("ERROR: Failed to allocate memory \n");
    exit (1);}

  jnB=(int *) malloc(NNMAX*n*sizeof(int));
  if (jnB== NULL){
    printf("ERROR: Failed to allocate memory \n");
    exit (1);}

  jn2=(int *) malloc(NNMAX*n*sizeof(int));
  if (jn2== NULL){
    printf("ERROR: Failed to allocate memory \n");
    exit (1);}



/*------------------------------*/
/*Mass of each "star" in DISC AND CORE in kg!!*/
/*-------------------*/

 printf("galaxy_v5.2.c\n");
 printf("==================================\n");
 printf("galaxy=%f - (0 = And - 1 = NGC3198)\n",galaxy);
 printf("==================================\n");
 printf("\n");
 printf("==========Init Variables==========\n");
 printf("n=%d\n",n);
 printf("n0=%d\n",n0);
 printf("dt=%f | %f (year)\n",dt,dt/year);
 printf("Mdisc=%e (sm)\n",Mdisc/sm);
 printf("Mcore=%e (sm)\n",Mcore/sm); 
 printf("tmax=%e (year)\n",tmax/year);
 printf("idump=%d\n",idump);
 printf("idumpmax=%d\n",idumpmax);
 printf("rmin= %e (kpc) drmin= %e (kpc^2)\n", rmin/kpc,drmin/(kpc*kpc));
/*Mass of each "star"*/

/*----------------------*/
/*BLACK HOLE mass + 1% of the core*/
m0 = 1.0e8*sm;
//m02 MILKY WAY
m02 = 1.0e8*sm;

/* DM core density*/
 if(galaxy == 2) {
   mcdm0=3.2e11*sm*0.8;
 }
 if(galaxy == 1) {
   mcdm0=3.2e11*sm;
     }
 if(galaxy == 0) {
   mcdm0=7.2e11*sm;
 }

 mcdm02 = 7.339172915e11*sm / 4;

 // 4x more CDM in galaxy 2
 // mcdm02=4.0*mcdm02;

 printf("m0=%e (sm) m02=%e (sm)\n",m0/sm,m02/sm);
 printf("mcdm0=%e (sm) mcdm02= %e (sm)\n",mcdm0/sm,mcdm02/sm);

/*Constant in NFW CDM distribution (units of kpc)*/
A=15.0;

/*gravity constant: N m^2/kg-2*/
/* so masses are in kg, and distances are in meters*/
g = 6.673e-11;

/* constant PI */
PI = 4.0*atan(1.0);


/* Initialize star masses in DISK and CORE */
//printf("n=%d n0=%d\n",n,n0);
// printf("setting masses..");



/* ===========================================*/
/* INitialize star locations and velocities*/
/* ===========================================*/

/* Minimum radius for DISK mass */
rdiskmin = 3.0*kpc;

/*DISK region goes from rdiskmin out to around 50kpc */
/* This is for an exponential: exp(-alpha*r) mass distribution */
// alpha = 0.125;
alpha = 0.15;
alpha2 = 0.75;

 printf("alpha=%f\n",alpha);
 printf("alpha2=%f\n",alpha2);
 printf("==================================\n");
 printf("\n");
 printf("==========Unit Variables==========\n");
 printf("sm=%e\n",sm);
 printf("day=%e\n",day);
 printf("year=%e\n",year);
 printf("au=%e\n",oneau);
 printf("kpc=%e\n",kpc);
 printf("==================================\n");
 printf("\n");
 printf("==========Flag Variables==========\n");
 printf("inn=%d\n",inn);
 printf("cdm=%f\n",cdm);
 printf("==================================\n");
 printf("\n");

 //Galaxy conditions at t= 4.000000e+03 (myr) and separation r= 1.475701e+02 (kpc)
 //x1= -111.215118 (kpc) y1=-45.009486 (kpc) vx1= 8.351124 (kpc/100myr) vy1= 0.709167 (kpc/100myr)
 //x2= 18.814076 (kpc) y2=24.771124 (kpc) vx2= -17.277489 (kpc/100myr) vy2= -3.585101 (kpc/100myr)


 dx = -111.215118*kpc;
 dy = -45.009486*kpc;

 vxg=8.351124*kpc/(myr*100.0);
 vxg2=-17.277489*kpc/(myr*100.0);

vyg=0.709167*kpc/(myr*100.0);
 vyg2=-3.585101*kpc/(myr*100.0);

 //vxg=5.0*kpc/(myr*100.0);
 //vxg2=-vxg;

 //vyg=5.0*kpc/(myr*100.0);
 //vyg2=-vyg;

 //vxg=0.0;
 //vyg=0.0;

vxg2=0.0;
 vyg2=0.0;

 dx2 = 18.814076*kpc;
 dy2 = 24.771124*kpc;
 if(my_pe == 0) { //host only does the initial distributions
 placeGalaxy(id,cdm,n,n0,dx,dy,r,x,y,alpha,alpha2,A,rdiskmin,m,mt,vx,vy,m0,Mdisc,Mcore,mcdm,mcdm0,vxg,vyg,drmin);

 placeGalaxy(id2,cdm,n,n0,dx2,dy2,r2,x2,y2,alpha,alpha2,A,rdiskmin,m2,mt2,vx2,vy2,m02,Mdisc2,Mcore2,mcdm2,mcdm02,vxg2,vyg2,drmin);
 

/*==============================================*/
/* Dump Mass distributons for plotting */
/*==============================================*/
     fp1=fopen("mass.mtv","w");
     fprintf(fp1,"$DATA=CURVE2D\n");
     fprintf(fp1,"%%toplabel= \"Time = %4.2f (myr)\"\n",0.0);
     fprintf(fp1,"%%xmin= -100\n");
     fprintf(fp1,"%%xmax= +100\n");
     fprintf(fp1,"%%ymin= -100\n");
     fprintf(fp1,"%%ymax= +100\n");
     fprintf(fp1,"%%equalscale= T\n");
     fprintf(fp1,"%%xlabel= \"x (kpc)\"\n");
     fprintf(fp1,"%%ylabel= \"y (kpc)\"\n");
     fprintf(fp1,"%%markertype= 13\n");
     fprintf(fp1,"%%markercolor= 3\n");
     fprintf(fp1,"%%linetype= 0\n");
     for (i=0;i<n;i++)

       {
	         fprintf(fp1,"%f %f\n",x[i]/kpc,y[i]/kpc);
		 fprintf(fp1,"%f %f\n",x2[i]/kpc,y2[i]/kpc);
       }
     fprintf(fp1,"\n");
     fclose(fp1);
 }
    
 MPI_Barrier(MPI_COMM_WORLD);
 printf("initial broadcast\n");
 printf("Host broadcasting arrays to nodes...\n");
 MPI_Bcast(x,n,MPI_DOUBLE,0,MPI_COMM_WORLD);
 MPI_Bcast(y,n,MPI_DOUBLE,0,MPI_COMM_WORLD);
 MPI_Bcast(vx,n,MPI_DOUBLE,0,MPI_COMM_WORLD);
 MPI_Bcast(vy,n,MPI_DOUBLE,0,MPI_COMM_WORLD);
 MPI_Bcast(r,n,MPI_DOUBLE,0,MPI_COMM_WORLD);
 MPI_Bcast(mcdm,n,MPI_DOUBLE,0,MPI_COMM_WORLD);
 MPI_Bcast(m,n,MPI_DOUBLE,0,MPI_COMM_WORLD);
 MPI_Bcast(x2,n,MPI_DOUBLE,0,MPI_COMM_WORLD);
 MPI_Bcast(y2,n,MPI_DOUBLE,0,MPI_COMM_WORLD);
 MPI_Bcast(vx2,n,MPI_DOUBLE,0,MPI_COMM_WORLD);
 MPI_Bcast(vy2,n,MPI_DOUBLE,0,MPI_COMM_WORLD);
 MPI_Bcast(r2,n,MPI_DOUBLE,0,MPI_COMM_WORLD);
 MPI_Bcast(mcdm2,n,MPI_DOUBLE,0,MPI_COMM_WORLD);
 MPI_Bcast(m2,n,MPI_DOUBLE,0,MPI_COMM_WORLD);

/*==============================================*/
/* Compute Effective central mass = All Mass within a radius rkpc */
/*==============================================*/
/* If nearest neighbor (inn=1), then */
/* update mass within radius r[i]. */
/* In order to speed up calculations, we */
/* don't do this every time step, only every idump */

    if (inn == 1) {
          for(i=0;i<n;i++)
	  {
	       rv=sqrt((x[i]-dx)*(x[i]-dx) + (y[i]-dy)*(y[i]-dy))/kpc;
	       rn = 0.5*kpc*(1.0 + rv/5.0);     

               mt[i]=0.0;
               for(j=0;j<n;j++)
	       {
		 xsq=(x[j] - x[i])*(x[j] - x[i]);
		 ysq=(y[j] - y[i])*(y[j] - y[i]);
		 dr=sqrt(xsq + ysq);
// if mass j is inside r[i] and OUTSIDE nearest neighbor radius rn
// then include it in the total mass within radius r[i]
                    if(r[j] < r[i] && dr > rn)
		    {
                         mt[i] = m[j] + mt[i];
		    }
	       }
	  }
	  for(i=0;i<n;i++)
	  {
	       rv=sqrt((x2[i]-dx2)*(x2[i]-dx2) + (y2[i]-dy2)*(y2[i]-dy2))/kpc;
	       rn = 0.5*kpc*(1.0 + rv/5.0);     

               mt2[i]=0.0;
               for(j=0;j<n;j++)
	       {
		 xsq=(x2[j] - x2[i])*(x2[j] - x2[i]);
		 ysq=(y2[j] - y2[i])*(y2[j] - y2[i]);
		 dr=sqrt(xsq + ysq);
// if mass j is inside r[i] and OUTSIDE nearest neighbor radius rn
// then include it in the total mass within radius r[i]
                    if(r2[j] < r2[i] && dr > rn)
		    {
                         mt2[i] = m2[j] + mt2[i];
		    }
	       }
	  }
	  for(i=0;i<n;i++)
	  {
               rv=sqrt((x[i]-dx2)*(x[i]-dx2) + (y[i]-dy2)*(y[i]-dy2));
	       rn = 0.5*kpc*(1.0 + rv/5.0);     

               mtb[i]=0.0;
               for(j=0;j<n;j++)
	       {
		 xsq=(x2[j] - x[i])*(x2[j] - x[i]);
		 ysq=(y2[j] - y[i])*(y2[j] - y[i]);
		 dr=sqrt(xsq + ysq);

// if mass j is inside r[i] and OUTSIDE nearest neighbor radius rn
// then include it in the total mass within radius r[i]
                    if(r2[j] < rv && dr > rn)
		    {
                         mtb[i] = m2[j] + mtb[i];
		    }
	       }
	  }
	  for(i=0;i<n;i++)
	  {
               rv=sqrt((x2[i]-dx)*(x2[i]-dx) + (y2[i]-dy)*(y2[i]-dy));
	       rn = 0.5*kpc*(1.0 + rv/5.0);     

               mt2b[i]=0.0;
               for(j=0;j<n;j++)
	       {
		 xsq=(x[j] - x2[i])*(x[j] - x2[i]);
		 ysq=(y[j] - y2[i])*(y[j] - y2[i]);
		 dr=sqrt(xsq + ysq);

// if mass j is inside r[i] and OUTSIDE nearest neighbor radius rn
// then include it in the total mass within radius r[i]
                    if(r[j] < rv && dr > rn)
		    {
                         mt2b[i] = m[j] + mt2b[i];
		    }
	       }
	  }
    }


/*================================================*/
/* Compute Nearest Neighbor interactions if inn = 1*/
/*================================================*/

if (inn == 1 ) {
  printf("calling neighbor\n");
  neighbor(n,nnA,jnA,x,y,dx,dy,&ierr);
  if (*ierr != 0) {
    printf("1-ierr= %d\n",*ierr);
    exit(1);  }

  neighbor(n,nnB,jnB,x2,y2,dx2,dy2,&ierr);
  if (*ierr != 0) {
    printf("2-ierr= %d\n",*ierr);
    exit(1);  }

  neighborCollide(n,nn1,jn1,x,y,dx2,dy2,&ierr, x2, y2);
   if (*ierr != 0) {
    printf("3-ierr= %d\n",*ierr);
    exit(1);  }
 neighborCollide(n,nn2,jn2,x2,y2,dx,dy,&ierr,x,y);
  if (*ierr != 0) {
    printf("4-ierr= %d\n",*ierr);
    exit(1);  }
 }



/*================================================*/
/* Compute velocity rotation curve*/
/*================================================*/


 velrot(n,r,vx,vy,&vs,&rvs,x,y,dx,dy,vxg,vyg);
 velrot(n,r2,vx2,vy2,&vs2,&rvs2,x2,y2,dx2,dy2,vxg2,vyg2);

     fp1=fopen("velrot.mtv","w");
     fprintf(fp1,"%%xmin= 0\n");
     fprintf(fp1,"%%xmax= +30\n");
     fprintf(fp1,"%%ymin= 0\n");
     if(galaxy == 0) {
       fprintf(fp1,"%%ymax= +300\n");
     }
     if(galaxy == 1) {
       fprintf(fp1,"%%ymax= +200\n");
     }
     if(galaxy == 2) {
       fprintf(fp1,"%%ymax= +200\n");
     }
     fprintf(fp1,"%%xlabel= \"r (kpc)\"\n");
     fprintf(fp1,"%%ylabel= \"v (km/s)\"\n");
     fprintf(fp1,"%%markertype= 13\n");
     fprintf(fp1,"%%markercolor= 3\n");
     fprintf(fp1,"%%linecolor= 3\n");
     fprintf(fp1,"%%linetype= 1\n");

     for (i=0;i<41;i++)
       {
	 fprintf(fp1,"%f %f\n",rvs[i],vs[i]);
	 //printf("velrot: %f %f\n",rvs[i], vs[i]);
       }
     fprintf(fp1,"\n");
     fclose(fp1);
     
     fp1=fopen("velrot2.mtv","w");
     fprintf(fp1,"%%xmin= 0\n");
     fprintf(fp1,"%%xmax= +30\n");
     fprintf(fp1,"%%ymin= 0\n");
     if(galaxy == 0) {
       fprintf(fp1,"%%ymax= +300\n");
     }
     if(galaxy == 1) {
       fprintf(fp1,"%%ymax= +200\n");
     }
     if(galaxy == 2) {
       fprintf(fp1,"%%ymax= +200\n");
     }
     fprintf(fp1,"%%xlabel= \"r (kpc)\"\n");
     fprintf(fp1,"%%ylabel= \"v (km/s)\"\n");
     fprintf(fp1,"%%markertype= 13\n");
     fprintf(fp1,"%%markercolor= 3\n");
     fprintf(fp1,"%%linecolor= 3\n");
     fprintf(fp1,"%%linetype= 1\n");

     for (i=0;i<41;i++)
       {
	 fprintf(fp1,"%f %f\n",rvs2[i],vs2[i]);
       }
     fprintf(fp1,"\n");
     fclose(fp1);

/* Compute total energy */
     eng(&En0,vx,vy,n,m0,mt,mcdm,r,x,y,drmin,m,dx,dy);
     eng(&En02,vx2,vy2,n,m02,mt2,mcdm2,r2,x2,y2,drmin,m2,dx2,dy2);
     engT(&EnT0,vx,vy,n,m0,mt,mcdm,r,x,y,drmin,m,dx,dy,vx2,vy2,m02,mt2,mcdm2,r2,x2,y2,m2,dx2,dy2,vxg,vyg,vxg2,vyg2,A,cdm,mcdm0,mcdm02);

printf("Initial Energy En0= %e\n",En0);
printf("Initial Energy En02= %e\n",En02);
printf("Initial Energy En0 + En02= %e\n",En02+En0);
printf("Initial TOTAL Energy EnT= %e\n",EnT0);

/* Open file to store energy vs time */
     fpE=fopen("En.mtv","w");
     fprintf(fpE,"%%xmin= 0\n");
     fprintf(fpE,"%%xmax= +1000\n");
     fprintf(fpE,"%%ymin= 0\n");
     fprintf(fpE,"%%ymax= +2.0\n");
     fprintf(fpE,"%%xlabel= \"t (myr)\"\n");
     fprintf(fpE,"%%ylabel= \"E/E0\"\n");
     fprintf(fpE,"%%markertype= 13\n");
     fprintf(fpE,"%%markercolor= 0\n");
     fprintf(fpE,"%%linecolor= 0\n");
     fprintf(fpE,"%%linetype= 1\n");

     fpt=fopen("traj.mtv","w");
     fprintf(fpt,"$DATA=CURVE2D\n");
     fprintf(fpt,"%%toplabel= \"CDM Trajectories\"\n");
     fprintf(fpt,"%%xmin= -100\n");
     fprintf(fpt,"%%xmax= +100\n");
     fprintf(fpt,"%%ymin= -100\n");
     fprintf(fpt,"%%ymax= +100\n");
     fprintf(fpt,"%%equalscale= T\n");
     fprintf(fpt,"%%xlabel= \"x (kpc)\"\n");
     fprintf(fpt,"%%ylabel= \"y (kpc)\"\n");


/*===========================================*/
/* Start Galaxy evolution loop */
/*===========================================*/
t=0;
ik=0;

/* Initial Energy ratio = 1.0 */
fprintf(fpE,"%f %f\n",t,En0/En0);
fprintf(fpE,"%f %f\n",t,En02/En02);

 MPI_Barrier(MPI_COMM_WORLD);
 printf("starting initial getaccels..\n");
/* Compute initial Acclerations of all mass points */
 getaccel(n, inn, cdm, mcdm0, m0, A, r, x, y, Ax, Ay, mcdm, mt, m, jnA, nnA, drmin,dx,dy,Axg,Ayg,id,my_pe,nstar_pe);
 getaccel(n, inn, cdm, mcdm02, m02, A, r2, x2, y2, Ax2, Ay2, mcdm2, mt2, m2, jnB, nnB, drmin,dx2,dy2,Axg,Ayg,id2,my_pe,nstar_pe);

   getaccelCore(n,ict,cdm, mcdm0,m0,A, r, x, y, Ax, Ay, mcdm, mt, m, jnA, nnA, drmin,dx, dy,mcdm02, m02, r2,x2, y2, Ax2, Ay2, mcdm2, mt2, m2, jnB, nnB, dx2, dy2, Axg, Ayg);

   getaccelCollide(n,ict,cdm, mcdm0,m0,A, r, x, y, Ax, Ay, mcdm, mt, m, jnA, nnA, drmin,dx, dy,mcdm02, m02, r2,x2, y2, Ax2, Ay2, mcdm2, mt2, m2, jnB, nnB, dx2, dy2,Axg,Ayg,id2,my_pe,nstar_pe);
   getaccelCollide(n,ict,cdm, mcdm02,m02,A, r2, x2, y2, Ax2, Ay2, mcdm2, mt2, m2, jnB, nnB, drmin,dx2, dy2,mcdm0, m0, r,x, y, Ax, Ay, mcdm, mt, m, jnA, nnA, dx, dy,Axg,Ayg,id,my_pe,nstar_pe);


   //   printf("calling neighbor collide.. \n");
   getaccelNeighborCollide(n,ict,cdm, mcdm0,m0,A, r, x, y, Ax, Ay, mcdm, mtb, m, jn1, nn1, drmin,dx, dy,mcdm02, m02, r2,x2, y2, Ax2, Ay2, mcdm2, mt2b, m2, jn2, nn2, dx2, dy2,my_pe, nstar_pe);
   getaccelNeighborCollide(n,ict,cdm, mcdm02,m02,A, r2, x2, y2, Ax2, Ay2, mcdm2, mt2b, m2, jn2, nn2, drmin,dx2, dy2,mcdm0, m0, r,x, y, Ax, Ay, mcdm, mtb, m, jn1, nn1, dx, dy, my_pe, nstar_pe);

   if(my_pe != 0) {
     MPI_Send(&Ax[my_pe*nstar_pe], nstar_pe, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD);
     MPI_Send(&Ay[my_pe*nstar_pe], nstar_pe, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD);
     MPI_Send(&Ax2[my_pe*nstar_pe], nstar_pe, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD);
     MPI_Send(&Ay2[my_pe*nstar_pe], nstar_pe, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD);
     MPI_Recv(&Ax[0], nstar_pe, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD, &status);
     MPI_Recv(&Ay[0], nstar_pe, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD, &status);
     MPI_Recv(&Ax2[0], nstar_pe, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD, &status);
     MPI_Recv(&Ay2[0], nstar_pe, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD, &status);
   }
   if(my_pe == 0) {
     MPI_Recv(&Ax[(my_pe+1)*nstar_pe], nstar_pe,  MPI_DOUBLE, 1, 0, MPI_COMM_WORLD, &status);
     MPI_Recv(&Ay[(my_pe+1)*nstar_pe], nstar_pe,  MPI_DOUBLE, 1, 0, MPI_COMM_WORLD, &status);
     MPI_Recv(&Ax2[(my_pe+1)*nstar_pe], nstar_pe,  MPI_DOUBLE, 1, 0, MPI_COMM_WORLD, &status);
     MPI_Recv(&Ay2[(my_pe+1)*nstar_pe], nstar_pe,  MPI_DOUBLE, 1, 0, MPI_COMM_WORLD, &status);
     MPI_Send(&Ax[0], nstar_pe, MPI_DOUBLE, 1, 0, MPI_COMM_WORLD);
     MPI_Send(&Ay[0], nstar_pe, MPI_DOUBLE, 1, 0, MPI_COMM_WORLD);
     MPI_Send(&Ax2[0], nstar_pe, MPI_DOUBLE, 1, 0, MPI_COMM_WORLD);
     MPI_Send(&Ay2[0], nstar_pe, MPI_DOUBLE, 1, 0, MPI_COMM_WORLD);
   }
   MPI_Barrier(MPI_COMM_WORLD);   
// printf("done with init getaccel");
/*===================================*/
/* Start time evolution */
/*===================================*/


printf(">running program..\n");
//for(i=0;i<n;i++) {
// printf("Ax2[i]= %f   Ax[i]= %f \n",Ax2[i],Ax[i]);
// printf("Ay2[i]= %f   Ay[i]= %f \n",Ay2[i],Ay[i]);
//}

while (t<tmax)
{
     t=t+dt;
     //printf("vxN: %f\n",vx[n-1]);
/* idump determines how often to write results to disk */
     idump = idump + 1;

     for(j=0;j<n;j++)
     {
/* update the position of mass j */
          x[j] = x[j] + vx[j]*dt + 0.5*Ax[j]*dt*dt;
	  y[j] = y[j] + vy[j]*dt + 0.5*Ay[j]*dt*dt;
/* form 1st half of velocity update */
	  vx[j]=vx[j] + 0.5*Ax[j]*dt;
	  vy[j]=vy[j] + 0.5*Ay[j]*dt;
/* update the position of mass j */
          x2[j] = x2[j] + vx2[j]*dt + 0.5*Ax2[j]*dt*dt;
	  y2[j] = y2[j] + vy2[j]*dt + 0.5*Ay2[j]*dt*dt;
/* form 1st half of velocity update */
	  vx2[j]=vx2[j] + 0.5*Ax2[j]*dt;
	  vy2[j]=vy2[j] + 0.5*Ay2[j]*dt;
     }
     // printf("galaxy velocity update 1/2\n");
      dx = dx + vxg*dt + 0.5*Axg[0]*dt*dt;
      dy = dy + vyg*dt + 0.5*Ayg[0]*dt*dt;
      vxg=vxg + 0.5*Axg[0]*dt;
      vyg=vyg + 0.5*Ayg[0]*dt;
      dx2 = dx2 + vxg2*dt + 0.5*Axg[1]*dt*dt;
      dy2 = dy2 + vyg2*dt + 0.5*Ayg[1]*dt*dt;
      vxg2=vxg2 + 0.5*Axg[1]*dt;
      vyg2=vyg2 + 0.5*Ayg[1]*dt;

      //      printf("dx= %e dy = %e dx2 = %e dy2 = %e \n",dx,dy,dx2,dy2);
      //      printf("vxg= %e vyg = %e vxg2 = %e vyg2 = %e \n",vxg,vyg,vxg2,vyg2);
      //      printf("1st Half\n");
      
      if(my_pe != 0) {
	MPI_Send(&x[my_pe*nstar_pe], nstar_pe, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD);
	MPI_Send(&y[my_pe*nstar_pe], nstar_pe, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD);
	MPI_Send(&x2[my_pe*nstar_pe], nstar_pe, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD);
	MPI_Send(&y2[my_pe*nstar_pe], nstar_pe, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD);
	MPI_Recv(&x[0], nstar_pe, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD, &status);
	MPI_Recv(&y[0], nstar_pe, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD, &status);
	MPI_Recv(&x2[0], nstar_pe, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD, &status);
	MPI_Recv(&y2[0], nstar_pe, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD, &status);
      }
      if(my_pe == 0) {
	MPI_Recv(&x[(my_pe+1)*nstar_pe], nstar_pe,  MPI_DOUBLE, 1, 0, MPI_COMM_WORLD, &status);
	MPI_Recv(&y[(my_pe+1)*nstar_pe], nstar_pe,  MPI_DOUBLE, 1, 0, MPI_COMM_WORLD, &status);
	MPI_Recv(&x2[(my_pe+1)*nstar_pe], nstar_pe,  MPI_DOUBLE, 1, 0, MPI_COMM_WORLD, &status);
	MPI_Recv(&y2[(my_pe+1)*nstar_pe], nstar_pe,  MPI_DOUBLE, 1, 0, MPI_COMM_WORLD, &status);
	MPI_Send(&x[0], nstar_pe, MPI_DOUBLE, 1, 0, MPI_COMM_WORLD);
	MPI_Send(&y[0], nstar_pe, MPI_DOUBLE, 1, 0, MPI_COMM_WORLD);
	MPI_Send(&x2[0], nstar_pe, MPI_DOUBLE, 1, 0, MPI_COMM_WORLD);
	MPI_Send(&y2[0], nstar_pe, MPI_DOUBLE, 1, 0, MPI_COMM_WORLD);
      }
      MPI_Barrier(MPI_COMM_WORLD);
      
      
/* -------------------------------------------- */
/* Now compute the Ax and Ay at the new x and y 
/* So that the Velocity can be updated (velocity Verlet algorithm)
/* The effective accleration is the average of A(t) and A(t+dt)
/* This method is second order accurate and has good stability for galaxy type
/* simulations (better stability than Runge-Kutta and predictor corrector)
/* -------------------------------------------- */

      getaccel(n, inn, cdm, mcdm0, m0, A, r, x, y, Ax, Ay, mcdm, mt, m, jnA, nnA, drmin,dx,dy,Axg,Ayg,id,my_pe,nstar_pe);
      getaccel(n, inn, cdm, mcdm02, m02, A, r2, x2, y2, Ax2, Ay2, mcdm2, mt2, m2, jnB, nnB, drmin,dx2,dy2,Axg,Ayg,id2,my_pe,nstar_pe);

      getaccelCore(n,ict,cdm, mcdm0,m0,A, r, x, y, Ax, Ay, mcdm, mt, m, jnA, nnA, drmin,dx, dy,mcdm02, m02, r2,x2, y2, Ax2, Ay2, mcdm2, mt2, m2, jnB, nnB, dx2, dy2, Axg, Ayg);

      getaccelCollide(n,ict,cdm, mcdm0,m0,A, r, x, y, Ax, Ay, mcdm, mt, m, jnA, nnA, drmin,dx, dy,mcdm02, m02, r2,x2, y2, Ax2, Ay2, mcdm2, mt2, m2, jnB, nnB, dx2, dy2, Axg,Ayg,id2,my_pe,nstar_pe);
      getaccelCollide(n,ict,cdm, mcdm02,m02,A, r2, x2, y2, Ax2, Ay2, mcdm2, mt2, m2, jnB, nnB, drmin,dx2, dy2,mcdm0, m0, r,x, y, Ax, Ay, mcdm, mt, m, jnA, nnA, dx, dy,Axg,Ayg,id,my_pe,nstar_pe);
    
     //     printf("calling neighbor collide.. \n");
      getaccelNeighborCollide(n,ict,cdm, mcdm0,m0,A, r, x, y, Ax, Ay, mcdm, mtb, m, jn1, nn1, drmin,dx, dy,mcdm02, m02, r2,x2, y2, Ax2, Ay2, mcdm2, mt2b, m2, jn2, nn2, dx2, dy2,my_pe,nstar_pe);
      getaccelNeighborCollide(n,ict,cdm, mcdm02,m02,A, r2, x2, y2, Ax2, Ay2, mcdm2, mt2b, m2, jn2, nn2, drmin,dx2, dy2,mcdm0, m0, r,x, y, Ax, Ay, mcdm, mtb, m, jn1, nn1, dx, dy,my_pe,nstar_pe);

     if(my_pe != 0) {
       MPI_Send(&Ax[my_pe*nstar_pe], nstar_pe, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD);
       MPI_Send(&Ay[my_pe*nstar_pe], nstar_pe, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD);
       MPI_Send(&Ax2[my_pe*nstar_pe], nstar_pe, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD);
       MPI_Send(&Ay2[my_pe*nstar_pe], nstar_pe, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD);
       MPI_Recv(&Ax[0], nstar_pe, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD, &status);
       MPI_Recv(&Ay[0], nstar_pe, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD, &status);
       MPI_Recv(&Ax2[0], nstar_pe, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD, &status);
       MPI_Recv(&Ay2[0], nstar_pe, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD, &status);
     }
     if(my_pe == 0) {
       MPI_Recv(&Ax[(my_pe+1)*nstar_pe], nstar_pe,  MPI_DOUBLE, 1, 0, MPI_COMM_WORLD, &status);
       MPI_Recv(&Ay[(my_pe+1)*nstar_pe], nstar_pe,  MPI_DOUBLE, 1, 0, MPI_COMM_WORLD, &status);
       MPI_Recv(&Ax2[(my_pe+1)*nstar_pe], nstar_pe,  MPI_DOUBLE, 1, 0, MPI_COMM_WORLD, &status);
       MPI_Recv(&Ay2[(my_pe+1)*nstar_pe], nstar_pe,  MPI_DOUBLE, 1, 0, MPI_COMM_WORLD, &status);
       MPI_Send(&Ax[0], nstar_pe, MPI_DOUBLE, 1, 0, MPI_COMM_WORLD);
       MPI_Send(&Ay[0], nstar_pe, MPI_DOUBLE, 1, 0, MPI_COMM_WORLD);
       MPI_Send(&Ax2[0], nstar_pe, MPI_DOUBLE, 1, 0, MPI_COMM_WORLD);
       MPI_Send(&Ay2[0], nstar_pe, MPI_DOUBLE, 1, 0, MPI_COMM_WORLD);
     }
     MPI_Barrier(MPI_COMM_WORLD);   
   
     for(j=0;j<n;j++)
     {
	  vx[j]=vx[j] + 0.5*Ax[j]*dt;
	  vy[j]=vy[j] + 0.5*Ay[j]*dt;
	  vx2[j]=vx2[j] + 0.5*Ax2[j]*dt;
	  vy2[j]=vy2[j] + 0.5*Ay2[j]*dt;

     }

     /* Update galaxy center position  */
      vxg=vxg + 0.5*Axg[0]*dt;
      vyg=vyg + 0.5*Ayg[0]*dt;
      vxg2=vxg2 + 0.5*Axg[1]*dt;
      vyg2=vyg2 + 0.5*Ayg[1]*dt;


      //      printf("dx= %e dy = %e dx2 = %e dy2 = %e \n",dx,dy,dx2,dy2);
      //      printf("vxg= %e vyg = %e vxg2 = %e vyg2 = %e \n",vxg,vyg,vxg2,vyg2);
      //      printf("Axg0 = %e Ayg0= %e Axg1= %e Ayg1= %e \n",Axg[0],Ayg[0],Axg[1],Ayg[1]);
      //      printf("2nd Half\n"); 


/* Dump data it idumpmax count is reached */
     if (idump == idumpmax)
     {
       idump=0;
       //make sure host has the most current x y vx and vy for saving to file
       if(my_pe != 0) {
	MPI_Send(&x[my_pe*nstar_pe], nstar_pe, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD);
	MPI_Send(&y[my_pe*nstar_pe], nstar_pe, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD);
	MPI_Send(&x2[my_pe*nstar_pe], nstar_pe, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD);
	MPI_Send(&y2[my_pe*nstar_pe], nstar_pe, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD);
	MPI_Recv(&x[0], nstar_pe, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD, &status);
	MPI_Recv(&y[0], nstar_pe, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD, &status);
	MPI_Recv(&x2[0], nstar_pe, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD, &status);
	MPI_Recv(&y2[0], nstar_pe, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD, &status);
	MPI_Send(&vx[my_pe*nstar_pe], nstar_pe, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD);
	MPI_Send(&vy[my_pe*nstar_pe], nstar_pe, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD);
	MPI_Send(&vx2[my_pe*nstar_pe], nstar_pe, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD);
	MPI_Send(&vy2[my_pe*nstar_pe], nstar_pe, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD);
	MPI_Recv(&vx[0], nstar_pe, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD, &status);
	MPI_Recv(&vy[0], nstar_pe, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD, &status);
	MPI_Recv(&vx2[0], nstar_pe, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD, &status);
	MPI_Recv(&vy2[0], nstar_pe, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD, &status);
      }
      if(my_pe == 0) {
	MPI_Recv(&x[(my_pe+1)*nstar_pe], nstar_pe,  MPI_DOUBLE, 1, 0, MPI_COMM_WORLD, &status);
	MPI_Recv(&y[(my_pe+1)*nstar_pe], nstar_pe,  MPI_DOUBLE, 1, 0, MPI_COMM_WORLD, &status);
	MPI_Recv(&x2[(my_pe+1)*nstar_pe], nstar_pe,  MPI_DOUBLE, 1, 0, MPI_COMM_WORLD, &status);
	MPI_Recv(&y2[(my_pe+1)*nstar_pe], nstar_pe,  MPI_DOUBLE, 1, 0, MPI_COMM_WORLD, &status);
	MPI_Send(&x[0], nstar_pe, MPI_DOUBLE, 1, 0, MPI_COMM_WORLD);
	MPI_Send(&y[0], nstar_pe, MPI_DOUBLE, 1, 0, MPI_COMM_WORLD);
	MPI_Send(&x2[0], nstar_pe, MPI_DOUBLE, 1, 0, MPI_COMM_WORLD);
	MPI_Send(&y2[0], nstar_pe, MPI_DOUBLE, 1, 0, MPI_COMM_WORLD);
	MPI_Recv(&vx[(my_pe+1)*nstar_pe], nstar_pe,  MPI_DOUBLE, 1, 0, MPI_COMM_WORLD, &status);
	MPI_Recv(&vy[(my_pe+1)*nstar_pe], nstar_pe,  MPI_DOUBLE, 1, 0, MPI_COMM_WORLD, &status);
	MPI_Recv(&vx2[(my_pe+1)*nstar_pe], nstar_pe,  MPI_DOUBLE, 1, 0, MPI_COMM_WORLD, &status);
	MPI_Recv(&vy2[(my_pe+1)*nstar_pe], nstar_pe,  MPI_DOUBLE, 1, 0, MPI_COMM_WORLD, &status);
	MPI_Send(&vx[0], nstar_pe, MPI_DOUBLE, 1, 0, MPI_COMM_WORLD);
	MPI_Send(&vy[0], nstar_pe, MPI_DOUBLE, 1, 0, MPI_COMM_WORLD);
	MPI_Send(&vx2[0], nstar_pe, MPI_DOUBLE, 1, 0, MPI_COMM_WORLD);
	MPI_Send(&vy2[0], nstar_pe, MPI_DOUBLE, 1, 0, MPI_COMM_WORLD);
      }
      MPI_Barrier(MPI_COMM_WORLD);
      // if (my_pe == 0) {
       
       if(my_pe==0) {
       time (&timeraw);
       timeinfo = localtime (&timeraw);
       // printf("Time is: %s \n", asctime(timeinfo));

       t2=clock();
       elapsed = ((double) (t2-t1))/CLOCKS_PER_SEC;
       
	  printf("t (myr)= %f %%= %f\n",t/myr,100.0*t/tmax);
	  printf("Elapsed= %f (sec) | %f (min) || Localtime= %s\n",elapsed,elapsed/60,asctime(timeinfo));
       }


/*   -------------------------------------------*/
 //Write CDM trajectories to file

	  fprintf(fpt," \n");
	  fprintf(fpt,"%%markercolor=3\n");
	  fprintf(fpt,"%%markertype=13\n");
          fprintf(fpt,"%f %f\n",dx/kpc,dy/kpc);
	  fprintf(fpt," \n");
	  fprintf(fpt,"%%markercolor=4\n");
	  fprintf(fpt,"%%markertype=13\n");
          fprintf(fpt,"%f %f\n",dx2/kpc,dy2/kpc);
	  fprintf(fpt," \n");
	  fflush(fpt);

	  // TEST
	  //     getaccelCore(n,ict,cdm, mcdm0,m0,A, r, x, y, Ax, Ay, mcdm, mt, m, jnA, nnA, drmin,dx, dy,mcdm02, m02, r2,x2, y2, Ax2, Ay2, mcdm2, mt2, m2, jnB, nnB, dx2, dy2, Axg, Ayg);
	  //     printf("Axg= %e Ayg= %e Axg2= %e Ayg2= %e\n",Axg[0],Ayg[0],Axg[1],Ayg[1]);

/*   -------------------------------------------*/

	  //printf("t (myr)= %f %%= %f\n",t/myr,100.0*t/tmax);
          fprintf(fout,"Updating Mass...\n");
/* file name counter update */	  
      ik=ik+1;

/* If nearest neighbor (inn=1), then */
/* update mass within radius r[i]. */
/* In order to speed up calculations, we */
/* don't do this every time step, only every idump */
    if (inn == 1) {
          for(i=0;i<n;i++)
	  {
	    rv=sqrt((x[i]-dx)*(x[i]-dx) + (y[i]-dy)*(y[i]-dy))/kpc;
	       rn = 0.5*kpc*(1.0 + rv/5.0);     

               mt[i]=0.0;
               for(j=0;j<n;j++)
	       {
		 xsq=(x[j] - x[i])*(x[j] - x[i]);
		 ysq=(y[j] - y[i])*(y[j] - y[i]);
		 dr=sqrt(xsq + ysq);
// if mass j is inside r[i] and OUTSIDE nearest neighbor radius rn
// then include it/ in the total mass within radius r[i]
                    if(r[j] < r[i] && dr > rn)
		    {
                         mt[i] = m[j] + mt[i];
		    }
	       }
	  }
	  for(i=0;i<n;i++)
	  {
	    rv=sqrt((x2[i]-dx2)*(x2[i]-dx2) + (y2[i]-dy2)*(y2[i]-dy2))/kpc;
	       rn = 0.5*kpc*(1.0 + rv/5.0);     

               mt2[i]=0.0;
               for(j=0;j<n;j++)
	       {
		 xsq=(x2[j] - x2[i])*(x2[j] - x2[i]);
		 ysq=(y2[j] - y2[i])*(y2[j] - y2[i]);
		 dr=sqrt(xsq + ysq);
// if mass j is inside r[i] and OUTSIDE nearest neighbor radius rn
// then include it/ in the total mass within radius r[i]
                    if(r2[j] < r2[i] && dr > rn)
		    {
                         mt2[i] = m2[j] + mt2[i];
		    }
	       }
	  }
	  for(i=0;i<n;i++)
	  {
	       rv=sqrt((x[i]-dx2)*(x[i]-dx2) + (y[i]-dy2)*(y[i]-dy2));
	       rn = 0.5*kpc*(1.0 + rv/5.0);     

               mtb[i]=0.0;
               for(j=0;j<n;j++)
	       {
		 xsq=(x2[j] - x[i])*(x2[j] - x[i]);
		 ysq=(y2[j] - y[i])*(y2[j] - y[i]);
		 dr=sqrt(xsq + ysq);

// if mass j is inside r[i] and OUTSIDE nearest neighbor radius rn
// then include it in the total mass within radius r[i]
                    if(r2[j] < rv && dr > rn)
		    {
                         mtb[i] = m2[j] + mtb[i];
		    }
	       }
	  }
	  for(i=0;i<n;i++)
	  {
               rv=sqrt((x2[i]-dx)*(x2[i]-dx) + (y2[i]-dy)*(y2[i]-dy));
	       rn = 0.5*kpc*(1.0 + rv/5.0);     

               mt2b[i]=0.0;
               for(j=0;j<n;j++)
	       {
		 xsq=(x[j] - x2[i])*(x[j] - x2[i]);
		 ysq=(y[j] - y2[i])*(y[j] - y2[i]);
		 dr=sqrt(xsq + ysq);

// if mass j is inside r[i] and OUTSIDE nearest neighbor radius rn
// then include it in the total mass within radius r[i]
                    if(r[j] < rv && dr > rn)
		    {
                         mt2b[i] = m[j] + mt2b[i];
		    }
	       }
	  }
    }

/*================================================*/
/* UPdate Nearest Neighbors every idump */
/*================================================*/
if (inn == 1 ) {
  neighbor(n,nnA,jnA,x,y,dx,dy,&ierr);
  if (*ierr != 0) {
    printf("1-ierr= %d\n",*ierr);
    exit(1);}

  neighbor(n,nnB,jnB,x2,y2,dx2,dy2,&ierr);
  if (*ierr != 0) {
    printf("2-ierr= %d\n",*ierr);
    exit(1);}

   neighborCollide(n,nn1,jn1,x,y,dx2,dy2,&ierr,x2,y2);
  if (*ierr != 0) {
    printf("3-ierr= %d\n",*ierr);
    exit(1);}

  neighborCollide(n,nn2,jn2,x2,y2,dx,dy,&ierr,x,y);
  if (*ierr != 0) {
    printf("4-ierr= %d\n",*ierr);
    exit(1);}
  

 }
/*=========================================*/
/* compute and save total energy*/
/* Total energy conservation      */    

 eng(&En,vx,vy,n,m0,mt,mcdm,r,x,y,drmin,m,dx,dy);
 eng(&En2,vx2,vy2,n,m02,mt2,mcdm2,r2,x2,y2,drmin,m2,dx2,dy2);
 engT(&EnT,vx,vy,n,m0,mt,mcdm,r,x,y,drmin,m,dx,dy,vx2,vy2,m02,mt2,mcdm2,r2,x2,y2,m2,dx2,dy2,vxg,vyg,vxg2,vyg2,A,cdm,mcdm0,mcdm02);

 fprintf(fout,"Energy En/En0= %f\n",En/En0);
//fprintf(fpE,"%f %f\n",t/myr,En/En0);
 fprintf(fout,"Energy En2/En02= %f\n",En2/En02);
//fprintf(fpE,"%f %f\n",t/myr,En2/En02);
 fprintf(fout,"TOTAL EnT/EnT0= %f\n",EnT/EnT0);
fprintf(fpE,"%f %f\n",t/myr,EnT/EnT0);
// FLUSH the file to disk (to keep it up to date)
 fflush(fpE);


/*===============================================*/
/*Compute rotational velocity at this time step*/
 velrot(n,r,vx,vy,&vs,&rvs,x,y,dx,dy,vxg,vyg);
 velrot(n,r2,vx2,vy2,&vs2,&rvs2,x2,y2,dx2,dy2,vxg2,vyg2);
 if (my_pe == 0) {
     sprintf(fname,"data/mass_%04d",ik);
     printf("%s\n",fname);
     fp1=fopen(fname,"w");
     fprintf(fp1,"$DATA=CURVE2D\n");
     fprintf(fp1,"%%toplabel= \"Time = %4.2f (myr)\"\n",t/myr);
     fprintf(fp1,"%%xmin= -100\n");
     fprintf(fp1,"%%xmax= +100\n");
     fprintf(fp1,"%%ymin= -100\n");
     fprintf(fp1,"%%ymax= +100\n");
     fprintf(fp1,"%%equalscale= T\n");
     fprintf(fp1,"%%xlabel= \"x (kpc)\"\n");
     fprintf(fp1,"%%ylabel= \"y (kpc)\"\n");
     fprintf(fp1,"%%markertype= 13\n");
     fprintf(fp1,"%%markercolor= 3\n");
     fprintf(fp1,"%%linetype= 0\n");
     for (i=0;i<n;i++)
       {
	 fprintf(fp1,"%f %f\n",x[i]/kpc,y[i]/kpc);
       }
     fprintf(fp1,"\n");
     fprintf(fp1,"%%markertype= 13\n");
     fprintf(fp1,"%%markercolor= 6\n");
     fprintf(fp1,"%%linetype= 0\n");
     for (i=0;i<n;i++)
       {
	 fprintf(fp1,"%f %f\n",x2[i]/kpc,y2[i]/kpc);
       }
     fprintf(fp1,"\n");
     fflush(fp1);
     fclose(fp1);


     sprintf(fname,"data/vdat_%04d",ik);
     printf("%s\n",fname);

     fp1=fopen(fname,"w");
     fprintf(fp1,"%d %d\n",n,n);
     fprintf(fp1,"%e %e %e %e\n",dx/kpc,dy/kpc,dx2/kpc,dy2/kpc);

     for (i=0;i<n;i++)
       {
	 fprintf(fp1,"%e %e %e %e\n",x[i]/kpc,y[i]/kpc,vx[i],vy[i]);
       }
     for (i=0;i<n;i++)
       {
	 fprintf(fp1,"%e %e %e %e\n",x2[i]/kpc,y2[i]/kpc,vx2[i],vy2[i]);
       }

     fprintf(fp1,"\n");
     fclose(fp1);

    
     sprintf(fname,"data/vrot_%04d",ik);
     printf("%s\n",fname);
     fp1=fopen(fname,"w");
     fprintf(fp1,"$DATA=CURVE2D\n");
     fprintf(fp1,"%%toplabel= \"Time = %4.2f (myr)\"\n",t/myr);
     fprintf(fp1,"%%xmin= 0\n");
     fprintf(fp1,"%%xmax= +30\n");
     fprintf(fp1,"%%ymin= 0\n");
     fprintf(fp1,"%%ymax= +300\n");
     fprintf(fp1,"%%xlabel= \"r (kpc)\"\n");
     fprintf(fp1,"%%ylabel= \"v (km/s)\"\n");
     fprintf(fp1,"%%markertype= 13\n");
     fprintf(fp1,"%%linetype= 1\n");
     for (i=0;i<41;i++)
       {
	 fprintf(fp1,"%f %f\n",rvs[i],vs[i]);
       }
     fprintf(fp1,"\n");
     fclose(fp1);

     sprintf(fname,"data/vrt2_%04d",ik);
     printf("%s\n",fname);

     fp1=fopen(fname,"w");
     fprintf(fp1,"$DATA=CURVE2D\n");
     fprintf(fp1,"%%toplabel= \"Time = %4.2f (myr)\"\n",t/myr);
     fprintf(fp1,"%%xmin= 0\n");
     fprintf(fp1,"%%xmax= +30\n");
     fprintf(fp1,"%%ymin= 0\n");
     fprintf(fp1,"%%ymax= +300\n");
     fprintf(fp1,"%%xlabel= \"r (kpc)\"\n");
     fprintf(fp1,"%%ylabel= \"v (km/s)\"\n");
     fprintf(fp1,"%%markertype= 13\n");
     fprintf(fp1,"%%linetype= 1\n");
     for (i=0;i<41;i++)
       {
	 fprintf(fp1,"%f %f\n",rvs2[i],vs2[i]);
       }
     fprintf(fp1,"\n");
     fclose(fp1);
 }

/*========================================*/
/* Dump data to disk here */
/*========================================*/

     

/* Dump velocities  to disk */
     


/*========================================*/
 MPI_Barrier(MPI_COMM_WORLD);
     }
}


/* Compute velocity rotation curve AT END OF SIMULATION*/
/* Average all the star velocities with +0.5 and -0.5 kpc*/ 
/* of radius = rad*/
 velrot(n,r,vx,vy,&vs,&rvs,x,y,dx,dy,vxg,vyg);

     fp1=fopen("velrotF.mtv","w");
     fprintf(fp1,"%%xmin= 0\n");
     fprintf(fp1,"%%xmax= +30\n");
     fprintf(fp1,"%%ymin= 0\n");
     fprintf(fp1,"%%ymax= +300\n");
     fprintf(fp1,"%%xlabel= \"r (kpc)\"\n");
     fprintf(fp1,"%%ylabel= \"v (km/s)\"\n");
     fprintf(fp1,"%%markertype= 13\n");
     fprintf(fp1,"%%markercolor= 3\n");
     fprintf(fp1,"%%linetype= 1\n");
     fprintf(fp1,"%%linecolor= 3\n");
     for (i=0;i<41;i++)
       {
	 fprintf(fp1,"%f %f\n",rvs[i],vs[i]);
       }
     fprintf(fp1,"\n");
     fclose(fp1);

/*==========================================*/
/*END OF MAIN PROGRAM*/
/*==========================================*/

}


/* ======================================================== *
/* SUBROUTINES BELOW */
/* ======================================================== */

/* Computes Total Energy returned in En*/
eng(double *En,double *vx,double *vy,int n,double m0, double *mt, double *mcdm,double *r, double *x,double *y,double drmin,double *m,double dx,double dy)
{
  int i,j;
  double vm2,vx2,vy2,xsq,ysq;
  double mass,g,U,dr,sm;

sm = 1.98892e30;
/*gravity constant: N m^2/kg-2*/
g = 6.673e-11;

// printf("m[0]= %f\n",*(m+0));

*En=0.0;
for(j=0;j<n;j++)
{
  vm2=vx[j]*vx[j] + vy[j]*vy[j];

/*CENTRAL mass = BLACK hole + CDM*/
  mass=m0 + mcdm[j];

  //  printf("eng: j= %d mass= %e (sm)\n",j,mass/sm);

/* Potential energy due to large central mass */
/* use drmin to limit force/potential when r -> 0 */
/* central force drmin is 1/4 drmin to better model core */
  U=-g*m[j]*mass/sqrt(r[j]*r[j] + drmin*0.1);

  //  printf("vm2= %f, mass= %f, U= %f\n",vm2,mass,U);

/* Add up potential energies due to all other mass points */
     for(i=0;i<n;i++)
     {
          if (i != j)
	  {
	 xsq=(x[j] - x[i])*(x[j] - x[i]);
  	 ysq=(y[j] - y[i])*(y[j] - y[i]);
	    dr=xsq + ysq + drmin;
            U=U-g*m[j]*m[i]/sqrt(dr);
	    //            printf("i=%d U= %f\n",i,U);
	  }
     }
     *En=*En + 0.5*m[j]*vm2 + U;
 }
}

// Computes the total energy of BOTH galaxies
engT(double *EnT,double *vx,double *vy,int n,double m0, double *mt, double *mcdm, double *r, double *x,double *y,double drmin,double *m,double dx,double dy,double *vx2,double *vy2, double m02, double *mt2,double *mcdm2, double *r2, double *x2, double *y2, double *m2, double dx2, double dy2, double vxg, double vyg, double vxg2, double vyg2,double A, double cdm, double mcdm0, double mcdm02)
{
  int i,j;
  double vm,vm2,xsq,ysq;
  double mass,mass2,g,U,dr,kpc;
  double oneau,rkpc,rkpc2,dr2;

/*gravity constant: N m^2/kg-2*/
g = 6.673e-11;
/* one kilo parsec*/
oneau = 149597871.0E3;
kpc = 2.06e8*oneau;

// printf("m[0]= %f\n",*(m+0));

*EnT=0.0;
for(j=0;j<n;j++)
{
  vm=vx[j]*vx[j] + vy[j]*vy[j];
  vm2=vx2[j]*vx2[j] + vy2[j]*vy2[j];

/*CENTRAL mass = BLACK hole + CDM*/
//  mass=m0 + mcdm[j];
//  mass2=m02 + mcdm2[j];

/* Potential energy due to large central mass */
/* use drmin to limit force/potential when r -> 0 */
/* central force drmin is 1/4 drmin to better model core */
  dr=sqrt(r[j]*r[j] + drmin*0.1);
  dr2=sqrt(r2[j]*r2[j] + drmin*0.1);

  /* Central Black hole energy */
  U=-g*m[j]*m0/dr;
  U=U -g*m2[j]*m02/dr2;

  /* Halo  potential energy */
    rkpc=dr/kpc;
    U=U-g*m[j]*mcdm0*log(1.0 + rkpc/A)/rkpc;

    rkpc2=dr2/kpc;
    U=U-g*m2[j]*mcdm02*log(1.0 + rkpc2/A)/rkpc2;

/* Add up potential energies due to all other mass points */
     for(i=0;i<n;i++)
     {
          if (i != j)
	  {
	    // GALAXY 1
	 xsq=(x[j] - x[i])*(x[j] - x[i]);
  	 ysq=(y[j] - y[i])*(y[j] - y[i]);
	    dr=xsq + ysq + drmin;
            U=U-g*m[i]*m[j]/sqrt(dr);
	    // GALAXY 2
	 xsq=(x2[j] - x2[i])*(x2[j] - x2[i]);
  	 ysq=(y2[j] - y2[i])*(y2[j] - y2[i]);
	    dr=xsq + ysq + drmin;
            U=U-g*m2[i]*m2[j]/sqrt(dr);
	  }
     }

// Don't add in potential yet
     *EnT=*EnT + 0.5*m[j]*vm + 0.5*m2[j]*vm2;
 }


/* Now add in potential energy between galaxies */

for(j=0;j<n;j++) {
      // GALAXY 1 star with GALAXY 2 core CDM
	 xsq=(x[j] - dx2)*(x[j] - dx2);
  	 ysq=(y[j] - dy2)*(y[j] - dy2);
	    dr=xsq + ysq;
            rkpc=sqrt(dr)/kpc;
//Compute mass of other galaxy's CDM + central black hole m0 at this distance
             U=U-g*m[j]*m02/sqrt(dr+drmin);
	     U=U-g*m[j]*mcdm02*log(1.0 + rkpc/A)/rkpc;

      // GALAXY 2 star with GALAXY 1 core CDM
	 xsq=(x2[j] - dx)*(x2[j] - dx);
  	 ysq=(y2[j] - dy)*(y2[j] - dy);
	    dr=xsq + ysq;
            rkpc=sqrt(dr)/kpc;
//Compute mass of other galaxy's CDM + central black hole m0 at this distance
             U=U-g*m2[j]*m0/sqrt(dr+drmin);
	     U=U-g*m2[j]*mcdm0*log(1.0 + rkpc/A)/rkpc;
}


for(j=0;j<n;j++) {
  for(i=0;i<n;i++) {
      // GALAXY 1 star with GALAXY 2 star
	 xsq=(x[j] - x2[i])*(x[j] - x2[i]);
  	 ysq=(y[j] - y2[i])*(y[j] - y2[i]);
         dr=xsq + ysq + drmin;
         U=U-g*m[j]*m2[i]/sqrt(dr);
  }
}

 printf("dx= %f dy=%f dx2=%f dy2=%f \n",dx/kpc,dy/kpc,dx2/kpc,dy2/kpc);
 
// finally, add in potential energy btween the two CDM cores
/* Use fixed size of Halo here to get effective total mass of halo for kinetic energy*/
rkpc=50.0;
mass=m0 + mcdm0*(-rkpc/(A + rkpc) - log(A) + log(A+rkpc));
mass2=m02 + mcdm02*(-rkpc/(A + rkpc) - log(A) + log(A+rkpc));

/* Compute distance between cdm halos */
 dr=(dx-dx2)*(dx-dx2) + (dy-dy2)*(dy-dy2);
 rkpc=sqrt(dr)/kpc;
 printf("engT: core to core rkpc= %e\n",rkpc);

 /* account for varying mass of halo at shorter distances when they overlap */
 U=U - g*mcdm0*mcdm02*(1.0/(A+rkpc) + log(1+rkpc/A)*log(1+rkpc/A)/rkpc);
 /* Black hole in cores */
 U=U - g*m0*m02/sqrt(dr);

 // Now add in total potential energy and core kinetic energy to get total energy= KE + PE
 vm=vxg*vxg + vyg*vyg;
 vm2=vxg2*vxg2 + vyg2*vyg2;

 *EnT=*EnT +  0.5*mass*vm + 0.5*mass2*vm2 + U;
}

/*================================================*/
/* Compute Nearest Neighbor lists */
/* nn(i) is the number of nearest neighbors about ith mass point
/* jn stores the index for each neighbor about ith mass point
/*================================================*/

neighbor(int n,int *nn,int *jn,double *x, double *y, double dx, double dy, int *ierr)
{

#define NNMAX 200

 int i,j,ind,nnx;
 double rv,kpc,rn,oneau,xsq,ysq,dr;

/* one kilo parsec*/
oneau = 149597871.0E3;
kpc = 2.06e8*oneau;
*ierr=0;

nnx=0;

for(i=0;i<n;i++)
{
/*make NN range increase with radius*/     
  rv=sqrt((x[i]-dx)*(x[i]-dx) + (y[i]-dy)*(y[i]-dy))/kpc;
  rn = 0.2*kpc*(1.0 + rv/5.0);     

  nn[i]=0;

     for(j=0;j<n;j++)
     {
          if(j != i)
	  {
	    xsq=(x[j] - x[i])*(x[j] - x[i]);
    	    ysq=(y[j] - y[i])*(y[j] - y[i]);
	    dr=sqrt(xsq + ysq);
               if(dr < rn)
	       {
		 ind=i*NNMAX + nn[i];
		 jn[ind]=j;
                 nn[i]++;
		 if (nn[i] > NNMAX) {
		   printf("ERROR: nn > NNMAX! nn=%d\n", nn[i]);
                   printf("Terminating..\n");
		   *ierr=1;
                   exit(1);
                 }
                   
	       }
	  }
     }

     //     printf("i= %d nn= %d\n",i,nn[i]);

     if (nn[i] > nnx) {
       nnx=nn[i];
          }
}
 printf("Max nn = %d\n",nnx);
}

velrot(int n,double *r,double *vx,double *vy,double *vs,double *rvs,double *x,double *y,double dx, double dy,double vxg,double vyg)
{
int j,i,nvs,jp;
double rad,oneau,kpc,dnvs,v2,v2x,v2y;
double vxv,vyv;

/* one kilo parsec*/
oneau = 149597871.0E3;
kpc = 2.06e8*oneau;

/* use a finer grid from 1 - 5 kpc */
for(j=0;j<16;j++)
{
/* Add one to put on 1 - 5 kpc range */
  rad = ((double) j + 1.0)*0.25;
  rvs[j]=rad;
/* nvs counts the number of stars */
     nvs=0;
/* initialize velocity to zero */
     vs[j]=0;
/* loop over all stars and test if within +/- 0.5kpc of rad */
     for(i=0;i<n;i++)
     {
/* is star j within +/- 0.125 kpc of rad?*/
       if ((r[i] / kpc > rad - 0.125) && (r[i] / kpc < rad + 0.125)){
/* Project velocity of star j onto tangent direction */
// first subtract out translational motion of galaxy
	 vxv=vx[i]-vxg;
	 vyv=vy[i]-vyg;
	 v2x=(vxv)*(y[i]-dy)/(r[i]);
	 v2y=-(vyv)*(x[i]-dx)/(r[i]);
	 v2=v2x + v2y;
/* Add up the tangential velocities */
	 vs[j]=vs[j]+v2;
/* count the number of stars in this range*/
	 nvs=nvs+1;
        }
     }

     if(nvs > 0)
     {
/* COMPUTE AVERAGE tangential VELOCITY at radius rad by dividing by the total number stars there*/
/* divide by 1000 to get km/sec*/
          dnvs = (double) nvs;
          vs[j]=vs[j]/(1000.0*dnvs);
     }
     else {
/* set to previous value if no particles were inside this radius range */
          vs[j]=vs[j-1];
      }
}

/* use a coarser grid from 5-30kpc */
for(j=0;j<26;j++)
{
  jp=j+16;
  rad = ((double) j)*1.0 + 5.0;
  rvs[jp]=rad;
/* nvs counts the number of stars */
     nvs=0;
/* initialize velocity to zero */
     vs[jp]=0;
/* loop over all stars and test if within +/- 0.5kpc of rad */
     for(i=0;i<n;i++)
     {
/* is star j within +/- 0.125 kpc of rad?*/
       if ((r[i] / kpc > rad - 0.5) && (r[i] / kpc < rad + 0.5)){
/* Project velocity of star j onto tangent direction */
// first subtract out translational motion of galaxy
	 vxv=vx[i]-vxg;
	 vyv=vy[i]-vyg;
	 v2x=(vxv)*(y[i]-dy)/(r[i]);
	 v2y=-(vyv)*(x[i]-dx)/(r[i]);
	 v2=v2x + v2y;
/* Add up the tangential velocities */
	 vs[jp]=vs[jp]+v2;
/* count the number of stars in this range*/
         nvs=nvs+1;
        }
     }

     if(nvs > 0)
     {
/* COMPUTE AVERAGE tangential VELOCITY at radius rad by dividing by the total number stars there*/
/* divide by 1000 to get km/sec*/
          dnvs = (double) nvs;
          vs[jp]=vs[jp]/(1000.0*dnvs);
     }
     else {
/* Set to previous value if no stars are in this range of rad */
       vs[jp]=vs[jp-1];
     }
 }

}


getaccel(int n, int inn, double cdm, double mcdm0, double m0,double A, double *r, double *x, double *y, double *Ax, double *Ay, double *mcdm, double *mt, double *m, int *jn, int *nn, double drmin,double dx, double dy,double *Axg, double *Ayg, int id,int my_pe, int nstar_pe)
{

#define NMAX 500

      int i,j,iv,ind;
      double kpc,rkpc,rsq,oneau,mass,g,f;
      double ex,ey,exi,eyi,dr,Fintx,Finty,fint;
      double Fx,Fy,sm,fc,masc;

      double mass_50;

sm = 1.98892e30;

/*gravity constant: N m^2/kg-2*/
g = 6.673e-11;
/* one kilo parsec*/
oneau = 149597871.0E3;
kpc = 2.06e8*oneau;


// Evaluate total CDM mass at 50kpc for computing A=F/mass //
 rkpc=50.0;
 mass_50=m0+cdm*mcdm0*(-rkpc/(A + rkpc) - log(A) + log(A+rkpc));


// printf ("getaccel: n= %d cdm= %f mcdm0= %e (sm)\n",n,cdm,mcdm0/sm);

 Axg[id]=0.0;
 Ayg[id]=0.0;

 //for(j=0;j<n;j++)
 for(j=my_pe*nstar_pe;j<(my_pe+1)*nstar_pe;j++)
     {
       rsq = ((x[j]-dx)*(x[j]-dx)) + ((y[j]-dy)*(y[j]-dy));
          r[j]=sqrt(rsq);
          rkpc=r[j]/kpc;


/*THIS DARK MATTER DIST is NFW*/
/*static, spherically symmetric mass*/
/* mcdm is the integrated NFW CDM density x volume = total CDM within radius rkpc */
/* this gives an approximately linear mass increase in DISK region */
/* it becomes flat at very large radius to give a finite mass*/
          mcdm[j]=cdm*mcdm0*(-rkpc/(A + rkpc) - log(A) + log(A+rkpc));

/*Total mass inside radius rkpc= Black hole + core + disK + CDM*/
/* if nearest neighbor (inn=1), then central mass includes mt */
/* if N-body (inn=0), then don't include it since it is in the full interaction term */
          mass=m0 + mt[j] + mcdm[j]; 
// Core mass only
	  masc=m0 + mcdm[j];

	  //	  printf("rkpc= %f mass=%e (sm) mcdm= %e (sm)\n",rkpc,mass/sm,mcdm[j]/sm);

/* force due to gravity between mass point m[j] and central total mass, use drmin to cut off force if rsq ->0 */
/* central force drmin is 1/4 drmin to better model core */
	  f = (g*m[j]*mass)/(rsq+drmin*0.1);
// Compute force that core sees from jth star (back reaction)
	  fc= (g*m[j]*masc)/(rsq+drmin*0.1);
/* cosine and sine */
	  ex = (x[j]-dx) / r[j];
	  ey = (y[j]-dy) / r[j];

/* Compute NEAREST NEIGHBOR force INTERACTIONS if inn == 1 else do N-Body (all interactions)*/
          Fintx=0.0;
          Finty=0.0;
          if (nn[j] != 0)
	  {
               for(i=0;i<nn[j];i++)
	       {
                    ind=j*NNMAX + i;
                    iv=jn[ind];
                    dr=(x[j] - x[iv])*(x[j] - x[iv]) + (y[j]-y[iv])*(y[j]-y[iv]);
                    fint=g*m[j]*m[iv]/(dr+drmin);
                    dr=sqrt(dr);
                    exi = (x[j]-x[iv]) / dr;
                    eyi = (y[j]-y[iv]) / dr;
                    Fintx=-fint*exi + Fintx;
                    Finty=-fint*eyi + Finty;
		    //		    printf("i= %d ind= %d iv= %d fint/f= %f dr= %f\n",i,ind,iv,fint/f,dr/kpc);
	       }
	  }
	  
/* Total force due to Gravity of central mass + NN interactions*/
	  Fx = -f * ex   + Fintx;
	  Fy = -f * ey   + Finty;

//TEST centeral only....BKK
//           Fx = -f * ex;
//           Fy = -f * ey;
/* Total acceleration of mass m[j] */
	  Ax[j] = Fx / m[j];
	  Ay[j] = Fy / m[j];

// Back reaction of stars on CDM halo if they are OUTSIDE (r> 5kpc) galaxy 
// assume central core region is dense/symmetric enough to cancel out star forces on cdm core 
	  if (rkpc > 5.0) {
	  Axg[id] = +fc*ex/mass_50 + Axg[id];
	  Ayg[id] = +fc*ey/mass_50 + Ayg[id];
	     	  }
 
     }
}      

placeGalaxy(int id,double cdm, int n, int n0, double dx, double dy, double *r,double *x, double *y, double alpha, double alpha2, double A, double rdiskmin, double *m, double *mt, double *vx, double *vy,double m0, double Mdisc, double Mcore,double *mcdm,double mcdm0,double vxg, double vyg, double drmin)
{
 
  double oneau, kpc, rv, th, PI, dro, dro2, dr,rkpc,mv,mvc,delv,delvt,mass,vm;
  double g;
int flag,i,j;

oneau = 149597871.0E3;
kpc = 2.06e8*oneau;
PI = 4.0*atan(1.0);


g = 6.673e-11;
 
 printf("placeGalaxycalled with x= %f y= %f id= %d \n",dx/kpc, dy/kpc, id);
 printf("     >n= %d n0= %d \n",n,n0);
 printf("\n");
mv = (Mdisc / n0);
mvc = (Mcore / (n-n0));

for (i=0;i<n0;i++){
/* set all star masses the same in kg in DISC*/
    m[i] = mv;
}

for (i=n0;i<n;i++){
  /* set all star masses the same in kg in CORE*/
     m[i] = mvc;
}

 printf("     >Initializing mass distribution (GalaxyID= %d): Disk\n",id);
     for (i=0;i<n0;i++)
     {
       rv = ((double) (rand() %10000))/10000.1;
       //	  printf("i= %d rv= %f\n",i,rv);
       r[i]=rdiskmin + kpc*(-1.0/alpha)*log(1.0-rv);
	  //	  printf("i= %d r= %f\n",i,r[i]);
          rv = ((double) (rand() %10000))/10000.1;
          th=rv*2.0*PI;
	  x[i]=r[i]*cos(th) + dx;
          y[i]=r[i]*sin(th) + dy;
	  // printf("x[i] = %f \n",x[i]/kpc);
	  //	  printf("%f  %f \n",x[i]/kpc,y[i]/kpc);
     }
     
/* Minimum distance in meters between two mass points in disk region*/
     dro=0.05*kpc;
     flag = 1;
     
     while(flag==1)
     {
          flag=0;
          printf("     >Scanning stars if too close\n");

          for(i=0;i<n0;i++)
	  {
               for(j=0;j<n0;j++)
	       {
                    if(j!=i)
		    {
                         dr=sqrt((x[i]-x[j])*(x[i]-x[j]) + (y[i]-y[j])*(y[i]-y[j]));

                         if(dr<dro)
			 {
			      rv = ((double) (rand() %10000))/10000.1;
                              r[i]=rdiskmin + kpc*(-1.0/alpha)*log(1.0-rv);
			      rv = ((double) (rand() %10000))/10000.1;
                              th=rv*2.0*PI;
                              x[i]=r[i]*cos(th) + dx;
                              y[i]=r[i]*sin(th) + dy;
                              flag=1;
			 }
		    }
	       }
	  }
     }

     
/* CORE region starts at 0.25kpc and goes out to 2.0kpc */
/* This is for an exponential: exp(-alpha*r) mass distribution */
     printf("     >Initializing mass distribution (GalaxyID= %d): Core\n",id);
     
     for(i=n0;i<n;i++)
     {
          rv = ((double) (rand() %10000))/10000.1;
	  //	  printf("i= %d rv= %f\n",i,rv);
	  //          r[i]=0.25*kpc + kpc*(-1.0/alpha2)*log(1.0-rv);
	  r[i]=0.25*kpc + kpc*(-1.0/alpha2)*log(1.0-rv);
          rv = ((double) (rand() %10000))/10000.1;
          th=rv*2.0*PI;
          x[i]=r[i]*cos(th) + dx;
          y[i]=r[i]*sin(th) + dy;
     }

  //Minimum distance between two mass points in core region     
     //     dro=0.07*kpc;
     dro=0.01*kpc;
     flag = 1;
     while(flag==1)
     {
          flag=0;
          printf("     >Scanning stars if too close\n");

          for(i=n0;i<n;i++)
	  {
               for(j=n0;j<n;j++)
	       {
                    if(j!=i)
		    {
                         dr=sqrt((x[i]-x[j])*(x[i]-x[j]) + (y[i]-y[j])*(y[i]-y[j]));

                         if(dr<dro)
			 {
			      rv = ((double) (rand() %10000))/10000.1;
			      r[i]=0.25*kpc + kpc*(-1.0/alpha2)*log(1.0-rv);
			      rv = ((double) (rand() %10000))/10000.1;
                              th=rv*2.0*PI;
                              x[i]=r[i]*cos(th) + dx;
                              y[i]=r[i]*sin(th) + dy;
                              flag=1;
			 }
		    }
	       }
	  }
	  
     }
     
for(i=0;i<n;i++) {
  mt[i]=0.0;
/*THIS DARK MATTER DIST is NFW*/
/*static, spherically symmetric mass*/
/* mcdm is the integrated NFW CDM density x volume = total CDM within radius rkpc */
     rkpc=r[i]/kpc;
     /* this gives an approximately linear mass increase in DISK region */
     /* it becomes flat at very large radius*/
     mcdm[i]=cdm*mcdm0*(-rkpc/(A + rkpc) - log(A) + log(A+rkpc));

/* Compute total galaxy mass within radius r[i] */     
/*n0 to n is only CORE!!
/*0 to n0 is disk only
/*0 to n is BOTH core + disc*/
     for(j=0;j<n;j++)
     {
          if(r[j]<r[i])
	  {
	       mt[i]=m[j]+mt[i];
	  }
     }
}

for(i=0;i<n;i++)
{
  ///* total mass is black hole + core + disk + CDM*/
  mass=m0 + mt[i] + mcdm[i];

/* DARK MATTER ONLY*/
  //  mass=mcdm[i];
/* BLACK HOLE + DISK AND/OR CORE*/
//  mass=m0+mt[i];
/* DISK only*/
//  mass=mt[i];

/* Equilibrium velocity for uniform circular motion
   about the total mass */
//     vm=sqrt(g*mass/r[i]);
/* include drmin in order to set better initial velocities near center */
     dr=r[i]*r[i] + drmin*0.1;
     vm=sqrt(r[i]*g*mass/dr);

/* Now Add in random radial velocity dispersion*/
/* to DISK only */
     rv = ((double) (rand() %10000))/10000.1;
     rv = rv * 2.0 - 1.0;
     if(r[i]>3.0*kpc && r[i]<60.0*kpc)
     {
       delv=rv*40000.0*(60.0*kpc - r[i] )/(57.0*kpc);
     }
     else
     {
          delv=0.0;
     }
/* Tangential dispersion */
/* to DISK only */
     rv = ((double) (rand() %10000))/10000.1;
     rv = rv * 2.0 -1.0;
     if(r[i]>3.0*kpc && r[i]<60.0*kpc)
     {
       delvt=rv*40000.0*(60.0*kpc-r[i])/(57.0*kpc);
     }
     else
     {
          delvt=0.0;
     }

/* Total velocity + random radial + random tangential dispersion */

     vx[i]=(vm*(y[i]-dy)/r[i] + delv*(x[i]-dx)/r[i] + delvt*(y[i]-dy)/r[i]) + vxg;  
     vy[i]=(-vm*(x[i]-dx)/r[i] + delv*(y[i]-dy)/r[i] - delvt*(x[i]-dx)/r[i]) + vyg;
     //     printf("vx[i] = %f vm= %f r[i] = %f delv= %f\n",vx[i],vm,r[i],delv);
     //     printf("vy[i] = %f \n",vy[i]);
}
//     printf("\n");
}


getaccelCollide(int n, int ict, double cdm, double mcdm0, double m0,double A, double *r, double *x, double *y, double *Ax, double *Ay, double *mcdm, double *mt, double *m, int *jn, int *nn, double drmin,double dx, double dy,double mcdm02, double m02, double *r2, double *x2, double *y2, double *Ax2, double *Ay2, double *mcdm2, double *mt2, double *m2, int *jn2, int *nn2, double dx2, double dy2, double *Axg, double *Ayg,int id,int my_pe,int nstar_pe)
{

#define NNMAX 200

      int i,j,iv,ind;
      double kpc,rkpc,rsq,oneau,mass,g,f;
      double ex,ey,exi,eyi,dr,Fintx,Finty,fint;
      double Fx,Fy,Fx2,Fy2,ex2,ey2,Fintx2,Finty2;
      double mass2,mass_50;
      int inn = 1;

/*gravity constant: N m^2/kg-2*/
g = 6.673e-11;
/* one kilo parsec*/
oneau = 149597871.0E3;
kpc = 2.06e8*oneau;

// Evaluate total CDM mass at 50kpc for computing A=F/mass //
// the Other galaxy here
 rkpc=50.0;
 mass_50=m02+cdm*mcdm02*(-rkpc/(A + rkpc) - log(A) + log(A+rkpc));

	  if (id != 0 && id != 1) {
          printf("ERROR in getaccelCollide id out of range: id= %d\n",id);
	  exit(1);}

// Loop over j stars in galaxy 1 and compute interactions with other galaxy's CDM
	  // for(j=0;j<n;j++)
	  for(j=my_pe*nstar_pe;j<(my_pe+1)*nstar_pe;j++)
	       {
// Compute distance to center of other galaxy
		 dr=(x[j]-dx2)*(x[j]-dx2) + (y[j]-dy2)*(y[j]-dy2);
                 rkpc=sqrt(dr)/kpc;
//Compute mass of other galaxy's CDM + central black hole m0 at this distance
		 mass=m02+cdm*mcdm02*(-rkpc/(A + rkpc) - log(A) + log(A+rkpc));

  		 fint=g*mass*m[j]/(dr + drmin);
                    dr=sqrt(dr);
                    exi = (x[j]-dx2) / dr;
                    eyi = (y[j]-dy2) / dr;
// Total force due to Gravity between mass m[j] and other galaxy's CDM
                    Fx=-fint*exi;
                    Fy=-fint*eyi;

/* Total acceleration of mass m[j] add galaxy interaction term */
	  Ax[j] = Fx / m[j] + Ax[j];
	  Ay[j] = Fy / m[j] + Ay[j];
/* Add equal but opposite (negative) force to the core acceleration of galaxy 2 */
/* but we need to divide by galaxy core mass here */
	  Axg[id]=Axg[id]-Fx/mass_50;
          Ayg[id]=Ayg[id]-Fy/mass_50;

	       }
}

neighborCollide(int n,int *nn1,int *jn1,double *x, double *y, double dx2, double dy2, int *ierr, double *x2, double *y2)
{

#define NNMAX 200

 int i,j,ind,nnx;
 double rv,kpc,rn,oneau,xsq,ysq,dr;

/* one kilo parsec*/
oneau = 149597871.0E3;
kpc = 2.06e8*oneau;
*ierr=0;

nnx=0;

for(i=0;i<n;i++)
{
/*make NN range increase with radius*/     
/* here rv is defined relative to center of OTHER galaxy */
  rv=sqrt((x[i]-dx2)*(x[i]-dx2) + (y[i]-dy2)*(y[i]-dy2))/kpc;
  rn = 0.5*kpc*(1.0 + rv/5.0);     

  nn1[i]=0;

     for(j=0;j<n;j++)
     {

	    xsq=(x2[j] - x[i])*(x2[j] - x[i]);
    	    ysq=(y2[j] - y[i])*(y2[j] - y[i]);
	    dr=sqrt(xsq + ysq);
               if(dr < rn)
	       {
		 ind=i*NNMAX + nn1[i];
		 jn1[ind]=j;
                 nn1[i]++;
		 if (nn1[i] > NNMAX) {
		   printf("ERROR IN NEIGHBORCOLLIDE: nn > NNMAX\n");
                   printf("Terminating..\n");
		   *ierr=1;
                   exit(1);
                 }
	  }
     }
     if (nn1[i] > nnx) {
       nnx=nn1[i];
          }
}
 printf("Max nn = %d\n",nnx);
}   

getaccelNeighborCollide(int n, int ict, double cdm, double mcdm0, double m0,double A, double *r, double *x, double *y, double *Ax, double *Ay, double *mcdm, double *mtb, double *m, int *jn1, int *nn1, double drmin,double dx, double dy,double mcdm02, double m02, double *r2, double *x2, double *y2, double *Ax2, double *Ay2, double *mcdm2, double *mt2b, double *m2, int *jn2, int *nn2, double dx2, double dy2,int my_pe, int nstar_pe)
{

#define NNMAX 200

      int i,j,iv,ind;
      double kpc,rkpc,rsq,oneau,mass,g,f;
      double ex,ey,exi,eyi,dr,Fintx,Finty,fint;
      double Fx,Fy,Fx2,Fy2,ex2,ey2,Fintx2,Finty2;
      double mass2;
      int inn = 1;

/*gravity constant: N m^2/kg-2*/
g = 6.673e-11;
/* one kilo parsec*/
oneau = 149597871.0E3;
kpc = 2.06e8*oneau;
//for(j=0;j<n;j++) {
 for(j=my_pe*nstar_pe;j<(my_pe+1)*nstar_pe;j++) {
   rsq = ((x[j]-dx2)*(x[j]-dx2)) + ((y[j]-dy2)*(y[j]-dy2));
   f = (g*m[j]*mtb[j])/(rsq+drmin*0.1);
    /* cosine and sine */
    ex = (x[j]-dx2) / sqrt(rsq);
    ey = (y[j]-dy2) / sqrt(rsq);
    
    Fintx2=0.0;
    Finty2=0.0;
    if (nn1[j] != 0)
      {
	for(i=0;i<nn1[j];i++)
	  {
		  ind=j*NNMAX + i;
		  iv=jn1[ind];
		  dr=(x[j] - x2[iv])*(x[j] - x2[iv]) + (y[j]-y2[iv])*(y[j]-y2[iv]);
		  fint=g*m[j]*m2[iv]/(dr+drmin);
		  dr=sqrt(dr);
		  exi = (x[j]-x2[iv]) / dr;
		  eyi = (y[j]-y2[iv]) / dr;
		  Fintx2=-fint*exi + Fintx2;
		  Finty2=-fint*eyi + Finty2;
		}
	    }
	  Fx = -f * ex + Fintx2;
	  Fy = -f * ey + Finty2;
	  Ax[j] = Fx / m[j] + Ax[j];
	  Ay[j] = Fy / m[j] + Ay[j]; 
}
}

getaccelCore(int n, int ict, double cdm, double mcdm0, double m0,double A, double *r, double *x, double *y, double *Ax, double *Ay, double *mcdm, double *mt, double *m, int *jn, int *nn, double drmin,double dx, double dy,double mcdm02, double m02, double *r2, double *x2, double *y2, double *Ax2, double *Ay2, double *mcdm2, double *mt2, double *m2, int *jn2, int *nn2, double dx2, double dy2,double *Axg, double *Ayg)
{

#define NNMAX 200

      int i,j,iv,ind;
      double kpc,oneau,mass,mass2,g,f;
      double exi,eyi,dr,fint,rkpc;
      double Fx,Fy;
      double mass_50,mass2_50;
      int inn = 1;

/*gravity constant: N m^2/kg-2*/
g = 6.673e-11;
/* one kilo parsec*/
oneau = 149597871.0E3;
kpc = 2.06e8*oneau;
// Evaluate total CDM mass at 50kpc for computing A=F/mass //
 rkpc=50.0;
 mass_50=m0+cdm*mcdm0*(-rkpc/(A + rkpc) - log(A) + log(A+rkpc));
 mass2_50=m02+cdm*mcdm02*(-rkpc/(A + rkpc) - log(A) + log(A+rkpc));

// Compute distance to center of other galaxy
 dr=(dx-dx2)*(dx-dx2) + (dy-dy2)*(dy-dy2);
 rkpc=sqrt(dr)/kpc;
// Compute variable masses for force calculation
 mass=m0+cdm*mcdm0*(-rkpc/(A + rkpc) - log(A) + log(A+rkpc));
 mass2=m02+cdm*mcdm02*(-rkpc/(A + rkpc) - log(A) + log(A+rkpc));
 fint=g*mass*mass2/(dr+drmin);

 dr=sqrt(dr);
 exi = (dx-dx2) / dr;
 eyi = (dy-dy2) / dr;
// Total force due to Gravity between mass m[j] and other galaxy's CDM
// printf("fint= %e exi= %e eyi= %e\n",fint,exi,eyi);
 Fx=-fint*exi;
 Fy=-fint*eyi;
 /* Total acceleration of mass m[j] add galaxy interaction term */
 /* Already divided by mass above */
 Axg[0] = Fx/mass_50 + Axg[0];
 Ayg[0] = Fy/mass_50 + Ayg[0];
 // printf("Axg0= %e Ayg0= %e \n",Axg[0],Ayg[0]);
 Axg[1] = -Fx/mass2_50 + Axg[1];
 Ayg[1] = -Fy/mass2_50 + Ayg[1];
 //  printf("Axg1= %e Ayg1= %e \n",Axg[1],Ayg[1]);
  //  exit(1);
 //TEST
 // Axg[0]=0.0;
 // Ayg[0]=0.0;
 // Axg[1]=0.0;
 // Ayg[1]=0.0;
}