/* This program uses PAMR to setup the subdomain
 * division of certain region, and then uses
 * standard mpi calls to evolve a set of particles
 * by Ben, 555 spring 2006 assigment 4.
 *
 */
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <pamr.h>
#include <mpi.h>
#include <time.h>

#define DIM  2 // dimensionality of the problem
#define XMIN -1.0//0.0
#define XMAX 1.0
#define YMIN -1.0//0.0
#define YMAX 1.0
#define RHO 2     // refinment ratio
#define NGLEV 1   // number of grid levels
#define NX   25  // total number of grid points in x-direction
#define NY   25  // total number of grid points in y-direction
#define NTLEV 1   // number of timet levels
#define NMAXGRIDS 3 // maximum number of grids (1 for unigrid)
#define NMAXGF 1 // maximum number of grid functions
#define OUTMODULO 1 // save every OUTMODULO steps
///#define NSTEPS 120// number of evolution steps
#define REGRIDSTEP 5 // regriding period
#define NGHOSTS 1 // number of ghost zones
#define HFRAC 0.25 // defines the size of the grids in hierarchy
#define NUMP 5000   //number of particles
#define GF_F 0 // symbolic name of index for grid function "f"
#define max(a, b)  (((a) > (b)) ? (a) : (b))
#define min(a, b)  (((a) < (b)) ? (a) : (b))
#define PI 3.14159265358979323846

// declare the fortran functions - note the trailing underscore
void initialize_(int *Nx,int *Ny,real *x,real *y,real *f);
void calc_center(real *timet,real *arm,real *x0,real *y0,real *xc);
void newposition(int *Nx,int *Ny,int *pn,real *timet,real *pposx,real *pposy,real *x,real *y);
void save_particles_(real *timet, double *savex, double *savey,int *Np,int *rank);
void make_step_(real *dt, real *timet, double *px, double *py,int *vtype,int *Nx, int *Ny, real *x, real*y, real *vmx, real *vmy);


// define global variables
real *f[NTLEV];
int f_gfn[NTLEV];

real *x,*y;
real xmin,xmax,ymin,ymax;
int ghost_width[2*DIM],Nx,Ny,phys_bdy[2*DIM],size;
int Lf,rank;
real base_bbox[2*DIM],bbox[2*DIM],dx,dy,lambda,na,ne;
int base_shape[DIM];
real *pposx[NUMP],*pposy[NUMP]; //*mx[NUMP], *my[NUMP];
// set the gfn's for all the variables
void set_gfns(void)
{
  int lev;
  for(lev=1;lev<=NTLEV;lev++) {
    if ((f_gfn[lev-1]=PAMR_get_gfn("f",PAMR_AMRH,lev))<0) perror("set_gnfs error");
  }
}







// load pointers
void ldptr(void)
{
   int dim,ngfs,lev,shape[DIM];
   real t,*x0[DIM],*gfs[NMAXGF*NTLEV],dx0[DIM];
   static int first=1;

   if (first) // call only once
   {
      first=0;
      set_gfns();
   }
   if (!(PAMR_get_g_attribs(&rank,&dim,shape,bbox,ghost_width,&t,&ngfs,x0,gfs)))
      perror("ldptr: PAMR_get_g_attribs failed");

//   for(lev=0;lev<NTLEV;lev++) {
     f[0]=gfs[f_gfn[0]-1];

//   }
//pointer to local perimeter coordinate arrays
   x=x0[0]; dx=x[1]-x[0];
   y=x0[1]; dy=y[1]-y[0];

//     PAMR_get_g_dim(&dim);
     PAMR_get_g_rank(&rank);
//     PAMR_get_g_shape(base_shape);
     PAMR_get_g_bbox(bbox);



   if ((bbox[0]-base_bbox[0])<dx/2) phys_bdy[0]=1; else phys_bdy[0]=0;
   if ((base_bbox[1]-bbox[1])<dx/2) phys_bdy[1]=1; else phys_bdy[1]=0;
   if ((bbox[2]-base_bbox[2])<dy/2) phys_bdy[2]=1; else phys_bdy[2]=0;
   if ((base_bbox[3]-bbox[3])<dy/2) phys_bdy[3]=1; else phys_bdy[3]=0;

   Nx=shape[0];
   Ny=shape[1];

   size=Nx*Ny;
}

// update the hierarchy bbox around the center defined by (xc[0],xc[1])
void update_gbox(int ngl,real *gbox,real *xc) {
  real width[2*DIM],pdx[2*DIM],temp,factor;
  int i,l;
//printf("dentro upoate rank=%i, xc1=%f,xc2=%f\n",rank,xc[1],xc[2]);


  for(i=0;i<DIM;i++) {
    width[i] = base_bbox[2*i+1]-base_bbox[2*i]; // width of the global bounding box
    pdx[i] = width[i]/(base_shape[i]-1); // mesh spacing of the parent grid (to be used
                                         // in the next loop)
    gbox[2*i] = base_bbox[2*i];          // initialize base grid
    gbox[2*i+1] = base_bbox[2*i+1];
  }
  factor = 0.5*HFRAC; // relative half width of the bounding box
  for(l=1;l<ngl;l++) { // loop over levels (starting from second level)
    for(i=0;i<DIM;i++) { // loop over dimensions
      temp = max(base_bbox[2*i],xc[i]-width[i]*factor);
      // align with the parent grid
      gbox[2*DIM*l+2*i] = base_bbox[2*i]+(int)((temp-base_bbox[2*i])/pdx[i]+0.5)*pdx[i];
      temp = min(base_bbox[2*i+1],xc[i]+width[i]*factor);
      // align with the parent grid
      gbox[2*DIM*l+2*i+1] = base_bbox[2*i+1]-(int)((base_bbox[2*i+1]-temp)/pdx[i]+0.5)*pdx[i];
      // adjust pdx for next level
      pdx[i] /=RHO;
    }
    factor *=HFRAC;
  }

}
// save variables (all levels)
void save_vars(real timet) {
  char *name[NMAXGF];
  int num_save,i,l;
  num_save = 1;
  name[0] = (char *) strdup("f");
  for(l=1;l<=Lf;l++) {
    for (i=0;i<num_save;i++) {
      if (!PAMR_save_gfn(name[i],PAMR_AMRH,1,l,timet,"","")) {
        perror("save_vars: error saving grid function");
      }
    }
  }
}
//////////////////////////////////////////////////////functions for array management
typedef struct {
   int n;
   void *arreglo;
} Arreglo;


Arreglo creaArreglo(int n) {
   Arreglo nuevo;
   int i;
   float *fa;

   nuevo.n = n;
   nuevo.arreglo = (void *) malloc(n * sizeof(float));
   fa = (float *) nuevo.arreglo;
   for (i=0; i < n; i++) {
     fa[i] = 0.1;
   }
   return nuevo;
}



void insertaArreglo(Arreglo a, int pos, float ele) {
   float *af;

   if (pos < a.n) {
     af = a.arreglo;
     af[pos] = ele;
   }
}

void MuestraArreglo(Arreglo a) {
   int i;
   float *af;

   af = a.arreglo;
   for (i=0; i < a.n; i++) {
     printf("%f %i \n", af[i],i);
   }
   printf("\n");
}




//MuestraElemento
float gete(Arreglo a,int n) {
           float *af;

           af = a.arreglo;
//	   printf("%f\n", af[n]);
//	   printf("\n");
	   return af[n];
}





void borraArreglo(Arreglo a) {
   free(a.arreglo);
}


/////////////////////////////////////////////////////////////////////////////////////////main

int main( int argc, char *argv[])
{
  int cp_rank,pamr_context,processors,Np,pn,totalN,dim,filtered,tss;
  real q,q_size,timet,dt;
  char *cp_file;
  int lev[NMAXGRIDS];
  char *varname;
  int in_amrh,in_mgh,num_tl,amr_inject,amr_interp;
  int amr_bdy_interp,amr_sync,mg_inject,mg_interp;
  int mg_sync,mg_noinj_to_amr,regrid_transfer,phys_bdy_type[2*DIM];

  int g_dim,g_shape[DIM],g_ghost_width[2*DIM],g_ngfs;

  int valid,Ng[DIM],min_lev,max_lev,num,bcast_flag,mig,newvalue,miga;

  int step,nsteps,outmodulo,tlev1,i,j,l,p,k,u,entero;
  float elemento,elem;

  real gbox[2*DIM*NMAXGRIDS],x0,y0,xc[DIM],arm,delta,out[2],pposx[NUMP],pposy[NUMP],total[NUMP],posxi,posyi,addy;

  Arreglo posx,posy,mx,my,temp,bax,bay,nuevasx,nuevasy,tsx,tsy;
  int *nRecv, *displ, *migsrecv, vtype, NUMP2, NSTEPS;
  float *recvArrayx, *recvArrayy,*mpimx, *mpimy, *newcx, *newcy;
  real  *savex, *savey, px, py, *vmx, *vmy;
  double starttime, endtime;

  char message[20];
  entero=0;
  base_bbox[0] = XMIN;
  base_bbox[1] = XMAX;
  base_bbox[2] = YMIN;
  base_bbox[3] = YMAX;
  base_shape[0] = NX;
  base_shape[1] = NY;

//  printf("a=%g,%g\n",XMIN,base_bbox[1]);
//  printf("b=%g,%g\n",NY,base_shape[1]);


  q = 0;
  q_size = 0;
  cp_file = (char *) NULL;
  cp_rank = -1;

  
  MPI_Init(&argc,&argv);

//timer for benchmarking
  starttime = MPI_Wtime();
  //
  //

///////////////////////////////////////////////////////////////reading parameters from file


if(!get_int_param(argv[1],"vtype",&vtype,1));
if(!get_int_param(argv[1],"NUMP2",&NUMP2,1));
if(!get_int_param(argv[1],"NSTEPS",&NSTEPS,1));

///////////////////////////////////////////////////////////////////
//
MPI_Comm_size(MPI_COMM_WORLD,&processors);
// MPI_Comm_rank(MPI_COMM_WORLD,&rank);
//printf("a=%g\n",base_bbox[0]);
//printf("b=%g\n",base_shape[1]);
Np=NUMP2/processors;
pn=NUMP2;


/////////////////////////////////////////////////////////////////////
pamr_context = PAMR_init_context(0,q_size,DIM,base_shape,base_bbox,cp_file,cp_rank);
  if (!pamr_context) {
    perror("Error initializing pamr_context");
    exit(-1);
  }
  // define variable attributes for unigrid (no transfers necessary)
  varname = (char *) strdup("f");
  in_amrh = 1;
  in_mgh = 0;
  num_tl = NTLEV;
  amr_inject = PAMR_STRAIGHT_INJECT;
//  amr_interp = PAMR_SECOND_ORDER;
//  amr_bdy_interp = PAMR_SECOND_ORDER;
  amr_interp = PAMR_NO_INTERP;
   amr_bdy_interp = PAMR_NO_INTERP;
  amr_sync = PAMR_SYNC;
  mg_inject = PAMR_NO_INJECT;
  mg_interp = PAMR_NO_INTERP;
  mg_sync = PAMR_NO_SYNC;
  mg_noinj_to_amr = 0;
regrid_transfer = PAMR_NO_TRANSFER;
//  regrid_transfer = PAMR_SECOND_ORDER;

  // note that sgfn starts at 1 therefore we subtract -1 for C-style arrays
  PAMR_def_var_full(varname,in_amrh,in_mgh,num_tl,amr_inject,
              amr_interp,amr_bdy_interp,amr_sync,mg_inject,mg_interp,
              mg_sync,mg_noinj_to_amr,regrid_transfer,phys_bdy_type)-1;
  Ng[0] = NGHOSTS;
  Ng[1] = NGHOSTS;
  PAMR_set_ghost_width(Ng);
  // for unigrid we have only level = 1 grid
  min_lev = 1;
  max_lev = NGLEV;
//  num = NGLEV;
  num=1;

//  for(i=0;i<NGLEV;i++) {
//    lev[i] = i+1;
lev[0]=1;

  timet = 0.0;
  x0 = 0.5*(XMIN+XMAX);
  y0 = 0.5*(YMIN+YMAX);
  arm = 0.25*(XMAX-XMIN);
  delta = 0.1*(XMAX-XMIN);

  // calculate the center of the gaussian xc
 //     printf("antes xc=%i, xc1=%f,xc2=%f\n",rank,xc[1],xc[2]);

//  calc_center_(&timet,&arm,&x0,&y0,xc);
//        printf("despues xc=%i, xc1=%f,xc2=%f\n",rank,xc[1],xc[2]);
  // create hierarchy gbox centered around xc
  update_gbox(NGLEV,gbox,xc);

  // compose the hierarchy
  PAMR_compose_hierarchy(min_lev,max_lev,num,lev,gbox,timet);
  // initialize all levels
  tlev1 = 0;
  Lf =  PAMR_get_max_lev(PAMR_AMRH);

  lambda = 0.01;
  dx = (XMAX-XMIN)/(NX-1);
  dy = (YMAX-YMIN)/(NY-1);
  dt = 0.01;//lambda*dx;

//  printf("texas Lf=%i\n",Lf);

  for(l=1;l<=Lf;l++) { // loop over levels;
    valid = PAMR_init_s_iter(l,PAMR_AMRH,0);
    while (valid) { // loop over grids at each level
      ldptr();

      initialize_(&Nx,&Ny,x,y,f[tlev1]);//,new);

// printf("texas Nx=%i,f=%f,rank=%i\n",Nx,f[tlev1],rank);



      valid = PAMR_next_g();
       }
//  }



  /*
  // print the boundaries assigned to each processor
      printf("rank=%i\n",rank);
      printf("box0=%f\n",bbox[0]);
      printf("box1=%f\n",bbox[1]);
      printf("box2=%f\n",bbox[2]);
      printf("box3=%f\n",bbox[3]);
  */
////////////////////////////////////////////////////Initialize particle positions in the subdomain randomly

//We need rank, otherwise each node has the same time and then it will
//use the same seed and we would have the same pseudorandom sequence
//i.e. same positions for the particles relative to their respective subdomians

srand48( time(NULL) +rank);

// Create two arrays that will contain the positions of the particles

 posx = creaArreglo(Np);
 posy = creaArreglo(Np);

// Initialize random positions in this subdomain

 for(i=0;i<=(Np-1);i++){

	posxi=(bbox[1]-bbox[0])*drand48()+bbox[0];
	posyi=(bbox[3]-bbox[2])*drand48()+bbox[2];
	insertaArreglo(posx,i,posxi);
	insertaArreglo(posy,i,posyi);

//  printf("initial positions: i=%i,posxi=%f,posyi=%f\n",i,posxi,posyi);

  //gete(posx,i),gete(posy,i));

//printf("test\n");
////printf("initpos %i %i %g %f %f\n",rank,i,timet,gete(posx,i),gete(posy,i));

//   posxi=0.0;
//   posyi=0.0;

 }
//// printf("end\n\n");

 //MuestraArreglo(posx);
 //MuestraArreglo(posy);


 mig=0;
 step = 0;

////////////////////////////////////////////////////////setup vel in a grid
// each processor calculates the velocity grid from its assigned subgrid
/// rxX[i],Nx, Ny, y[j]..
//	we use the same function as the analitic part, to create a grid.
	
/* old code no ghost cells

   vmx=(real *)malloc((Ny)*sizeof(real));
   vmy=(real *)malloc((Nx)*sizeof(real));


   for(j=0;j<=(Ny-1);j++){

           vmx[j]=2*PI*(-y[j]);

   }//end i for

   for(i=0;i<=(Nx-1);i++){
                  vmy[i]=2*PI*x[i];
   }//end j for

   /printf
   
   */  

 
 
 // two ghost zones on each side
vmx=(real *)malloc((Ny+4)*sizeof(real));
vmy=(real *)malloc((Nx+4)*sizeof(real));

vmx[0]=2*PI*(-(y[0]-2*dy));
vmx[1]=2*PI*(-(y[0]-dy));
vmx[Ny+2]=2*PI*(-(y[Ny-1]+dy));
vmx[Ny+3]=2*PI*(-(y[Ny-1]+2*dy));

for(j=2;j<=(Ny+1);j++){
		vmx[j]=2*PI*(-y[j-2]);

}//end i for

/*for(j=0;j<=(Ny+3);j++){
//printf("oaxaca %i, %f\n", j,vmx[j]);
}
*/
vmy[0]=2*PI*(x[0]-2*dx);
vmy[1]=2*PI*(x[0]-dx);
vmy[Nx+2]=2*PI*(x[Nx-1]+dx);
vmy[Nx+3]=2*PI*(x[Nx-1]+2*dx);

for(i=2;i<=(Nx+1);i++){
                vmy[i]=2*PI*x[i-2];
}//end j for

//for(i=0;i<=(Nx+3);i++){
//  printf("oaxaca %i, %f\n", i,vmy[i]);
//  }
 
/////////////////////////////////////////////each processor has its own mesh Nx,Ny,
//each one was initialized separately, this code shows it:

//printf("step=%i,timet=%g\n",step,timet);

/*
for(i=0;i<=(Nx-1);i++){
  printf("meshx Nx=%i.i=%i,x=%f,rank=%i\n",Nx,i,x[i],rank);
 }//endfor
 for(i=0;i<=(Ny-1);i++){
  printf("meshy Ny=%i.i=%i,y=%f,rank=%i\n",Ny,i,y[i],rank);
 }//endfor
*/

// for(i=0;i<=(Nx-1);i++){
// for(j=0;j<=(Ny-1);j++){
//  printf("meshx Nx=%i.i=%i,f=%f,rank=%i\n",Nx,i,f[tlev1],rank);
// }
// }
////printf("mesh for processor=%i done",rank);


  while (step<=NSTEPS) { //////////////////////////////////////////////////////////main '(time) evolution' loop





//////////////////////////////////////////////////////////////

//          PAMR_get_g_rank(&rank);
//          printf("rank afuera=%i\n",rank);



  if (rank ==0)
    //  	printf("step=%d, rank=%i\n",step,rank);
  if ((step % OUTMODULO) == 0){
////  if (rank ==0)
//         printf("SAVING at timet=%g\n",timet);
////   printf("test\n");
 //   save_vars(timet);
    }

////////////////////////////////////////////////////////////////







    //evolution sequence: time updated at the end of the loop

//////////////////////////filter out newcomers and reconstruct posx,poy, locals + newcomers, change Np too
//////////////////////////doesnt make sense for t=0, cuz we receive no migration array from root


if (timet!=0){// dont filter t=0

//printf("arreglos posx y posy antes del merge\n");

////	MuestraArreglo(posx);
////	MuestraArreglo(posy);

	//check the state of the arrays
	for(j=0;j<totalN;j++){

//	printf("grrr=%i,newcx=%f,newcy=%f,%g,%i,%i\n",j,newcx[j],newcy[j],timet,rank,totalN);

	}

	//any newcomers belong to this domain?
	filtered=0;
	for(i=0;i<totalN;i++){

		if(  newcx[i] > bbox[0] && newcx[i] < bbox[1] && newcy[i] > bbox[2] && newcy[i] < bbox[3])
		  	{
		   	filtered=filtered+1;
		  	}

	}//end for


// if newcomers belong to our domain
if (filtered!=0){

			tsx=creaArreglo(filtered);
			tsy=creaArreglo(filtered);
			tss=0;


	for(i=0;i<totalN;i++){

	        if(  newcx[i] > bbox[0] && newcx[i] < bbox[1] && newcy[i] > bbox[2] && newcy[i] < bbox[3])
		  	{
			insertaArreglo(tsx,tss,newcx[i]);
			insertaArreglo(tsy,tss,newcy[i]);
			tss=tss+1;
	 	}

	}//end for

// put newcomers belonging to us in array nuevasx
////	printf("here filtered=%i\n",filtered);


	nuevasx=creaArreglo(filtered+Np);
        nuevasy=creaArreglo(filtered+Np);


	for(i=0;i<=(filtered-1);i++){

                   	insertaArreglo(nuevasx,i,gete(tsx,i));
                        insertaArreglo(nuevasy,i,gete(tsy,i));


	}//end for


	for(i=(filtered);i<=(Np+filtered-1);i++){


	insertaArreglo(nuevasx,i,gete(posx,(i-filtered)));
	insertaArreglo(nuevasy,i,gete(posy,(i-filtered)));


	}//end for

        ////printf("here posx,nuevasx,posy,nuevasy,filtered=%i,Np=%i\n",filtered,Np);


	////MuestraArreglo(posx);
	////MuestraArreglo(nuevasx);

	////MuestraArreglo(posy);
	////MuestraArreglo(nuevasy);



	borraArreglo(posx);
	borraArreglo(posy);
	borraArreglo(tsx);
	borraArreglo(tsy);

	posx=creaArreglo(filtered+Np);
	posy=creaArreglo(filtered+Np);


		for(i=0;i<=(Np+filtered-1);i++){

			insertaArreglo(posx,i,gete(nuevasx,i));
			insertaArreglo(posy,i,gete(nuevasy,i));

		}//end for

	Np=filtered+Np;

	borraArreglo(nuevasx);
        borraArreglo(nuevasy);

}//end if filtered !=0


	free(newcx);
	free(newcy);

}///////////////////////////////////////////////////////////////////end filter out  if timet!=0


//subdomain partition by PAMR
  ////  printf("----------------------\n");
  ////  printf("rank=%i,timet=%g\n\n",rank,timet);

/*
    printf("box_0xmin=%f\n",bbox[0]);
    printf("box_1xmax=%f\n",bbox[1]);
    printf("box_2ymin=%f\n",bbox[2]);
    printf("box_3ymax=%f\n",bbox[3]);
*/
////printf("checkrank=%i,Np=%i,time=%g\n",rank,Np,timet);



/////////////////////////////////////////////////////////////////////save particles


	savex=(real *)calloc(Np, sizeof(real));
	savey=(real *)calloc(Np, sizeof(real));
	for(i=0;i<=(Np-1);i++){
		savex[i]=gete(posx,i);
		savey[i]=gete(posy,i);
		////printf("hola %f %f %g holarank=%i\n",savex[i],savey[i],timet,rank);
	}//end for
	save_particles_(&timet, savex, savey,&Np, &rank);
	free(savex);
	free(savey);
////////////////////////////////////////////////////



    for(i=0;i<=(Np-1);i++){ /////////////////////////////////////////begin first cycle thru particles


////////////////////////////////update

px=gete(posx,i);
py=gete(posy,i);

//printf("make %i %i %g %f %f\n",rank,i,timet,gete(posx,i),gete(posy,i));
	make_step_(&dt,&timet, &px, &py,&vtype, &Nx, &Ny, x, y, vmx, vmy);
			insertaArreglo(posx,i,px);
			insertaArreglo(posy,i,py);
//printf("make %i %i %g %f %f\n",rank,i,timet,gete(posx,i),gete(posy,i));

////////////////////////////////////////////////////


//      insertaArreglo(posx,i,gete(posx,i)+0.01);
//      printf("%i %i %g %f %f\n",rank,i,timet,gete(posx,i),gete(posy,i));
      //printf("after update rank=%i,particle=%i,timet=%g,pposx=%f,pposy=%f,\n\n\n",rank,i,timet,gete(posx,i),gete(posy,i));

      //check if the particles migrate and count their number

		if(gete(posx,i) < bbox[0] || gete(posx,i) > bbox[1] || gete(posy,i) < bbox[2] || gete(posy,i) > bbox[3])
		  {
		    mig=mig+1;
		  }



    }////////////////////////////////////////////////////////////////////////////end first cycle thru particles
    ////printf("mig= %i\n",mig);



	mx=creaArreglo(1);
	my=creaArreglo(1);
	    insertaArreglo(mx,0,XMAX+1);
	    insertaArreglo(my,0,YMAX+1);
dim=1;

    // if there are migrants we create an array to store the entries
    // of posx and posy that will migrate and the resize those
    if(mig!=0){
        borraArreglo(mx);
	borraArreglo(my);
	temp=creaArreglo(mig);
	mx=creaArreglo(mig);
	my=creaArreglo(mig);


	l=0.0;

	for(i=0;i<=(Np-1);i++){//2nd circle thru updated particles to fetch migrants

	 	if(gete(posx,i) < bbox[0] || gete(posx,i) > bbox[1] || gete(posy,i) < bbox[2] || gete(posy,i) > bbox[3])
		  {
		    insertaArreglo(temp,l,i);

		    l=l+1;
		  }

        }//end 2nd cycle particles
	//printf("Muestra arreglo temp,mig=%i,l=%i\n\n",mig,l);
	//MuestraArreglo(temp);

	for(j=0;j<=(mig-1);j++){
	  addy=gete(temp,j);
	  //	    printf("gete=%f\n",addy);

	    insertaArreglo(mx,j,gete(posx,addy));
	    insertaArreglo(my,j,gete(posy,addy));

			 }//end j loop

	////MuestraArreglo(mx);
	////MuestraArreglo(my);
	////printf("Np=%i,mig=%i\n",Np,mig);

	//now i adjust the posx and posy arrays to their new size excluding migrants
	////////////////////////////////////////////////////////////////////////////////resize
	bax=creaArreglo(Np-mig);
	p=0;
	k=0;
	for(i=0;i<=(Np-1);i++) {
       //printf("i=%i\n",i);
             for(u=0;u<=(mig-1);u++) {
	       //printf("i=%i,u=%i, muestra=%f\n",i,u,gete(temp,u));
                     if (i==gete(temp,u))
	             {
                     p=1;
		     }

	     }//end u loop
	     if (p==1)//means if there isa  match and this entry will migrate and be excluded from posx
			{
				////printf("do nothing\n");
			}
	     else//we keep this entry for the new array be which will become the new posx
			{
			elem=gete(posx,i);

				insertaArreglo(bax,k,elem);	//ba[k]==posx[u];
				k=k+1;
				//Debug// printf("elem %f, k=%i\n",elem,k);
			}

	    p=0;

     }//end i loop
	//////////////////////////////////////////////////////////////////////////////end resize
	////////////////////////////////////////////////////////////////////////////////resize
	bay=creaArreglo(Np-mig);
	p=0;
	k=0;
	for(i=0;i<=(Np-1);i++) {
       //printf("i=%i\n",i);
             for(u=0;u<=(mig-1);u++) {
	       //printf("i=%i,u=%i, muestra=%f\n",i,u,gete(temp,u));
                     if (i==gete(temp,u))
	             {
                     p=1;
		     }

	     }//end u loop
	     if (p==1)//means if there is a  match and this entry will migrate and be excluded from posx
			{
				////printf("do nothing\n");
			}
	     else//we keep this entry for the new array be which will become the new posx
			{
			elem=gete(posy,i);

				insertaArreglo(bay,k,elem);	//ba[k]==posy[u];
				k=k+1;
				//Debug//printf("elem %f, k=%i\n",elem,k);
			}

	    p=0;

     }//end i loop
	//////////////////////////////////////////////////////////////////////////////end resize
	////MuestraArreglo(bax);
	////MuestraArreglo(bay);


	Np=Np-mig;//We redefine the number of particles in the domain, Np-mig

	borraArreglo(posx);//erasing posx
	posx=creaArreglo(Np);//we create it again (reallocate memory)
	     for(i=0;i<=(Np-1);i++) {//refill it with bax
	        insertaArreglo(posx,i,gete(bax,i));
	     }
	     //same for posy
	borraArreglo(posy);//erasing posy
	posy=creaArreglo(Np);//we create it again (reallocate memory)
	     for(i=0;i<=(Np-1);i++) {//refill it with bay
	        insertaArreglo(posy,i,gete(bay,i));
	     }


   ////printf("new posx and posy\n");
   ////MuestraArreglo(posx);
   ////MuestraArreglo(posy);

//   printf("arreglo mx");
//   MuestraArreglo(mx);

dim=mig;
}//end if there are migrants... (mig)/////////////////////////////////////////////////

  nRecv=(int *)calloc(processors, sizeof(int));
  migsrecv=(int *)calloc(processors, sizeof(int));
  displ=(int *)calloc(processors, sizeof(int));

  ////printf("newrank=%i,mig=%i,dim=%i,time=%g\n",rank,mig,dim,timet);

  MPI_Gather(&dim, 1, MPI_INT, nRecv, 1, MPI_INT, 0,MPI_COMM_WORLD);

  MPI_Gather(&mig, 1, MPI_INT, migsrecv, 1, MPI_INT, 0,MPI_COMM_WORLD);

// calculate the displacement
if (rank==0){
  displ[0] = 0;
  for (i=1;i<processors;i++){
    displ[i] = nRecv[i-1] + displ[i-1];
  }
  totalN = displ[processors-1] + nRecv[processors-1];
  ////printf("totalN=%i\n",totalN);
 recvArrayx=(float *)calloc(totalN,sizeof(float));
 recvArrayy=(float *)calloc(totalN,sizeof(float));

}
MPI_Bcast(&totalN, 1, MPI_INT, 0, MPI_COMM_WORLD);

mpimx = (float *)calloc(dim,sizeof(float));
mpimy = (float *)calloc(dim,sizeof(float));

for(i=0;i<=(dim-1);i++){

mpimx[i]=gete(mx,i);
mpimy[i]=gete(my,i);
////printf("mpimx=%f,mpimy=%f\n",mpimx[i],mpimy[i]);

  }


MPI_Gatherv(mpimx, dim, MPI_FLOAT, recvArrayx, nRecv, displ, MPI_FLOAT, 0, MPI_COMM_WORLD);
MPI_Gatherv(mpimy, dim, MPI_FLOAT, recvArrayy, nRecv, displ, MPI_FLOAT, 0, MPI_COMM_WORLD);
newcx = (float *)calloc(totalN,sizeof(float));
newcy = (float *)calloc(totalN,sizeof(float));

if(rank==0){

for(i=0;i<totalN;i++){
  newcx[i]=recvArrayx[i];
  newcy[i]=recvArrayy[i];
  ////printf("recvArrayxy: time=%g,%i,%f,newcx=%f,%f,newcy=%f,%f,i=%i,\n",timet,totalN,recvArrayx[i],newcx[i],recvArrayy[i],newcy[i],i);


  }

for(i=0;i<=(processors-1);i++){
  ////printf("migrants from proc %i  are dim=%i,mig=%i,time=%g\n",i,nRecv[i],migsrecv[i],timet);
  }


free(recvArrayx);
free(recvArrayy);
free(nRecv);
free(migsrecv);
}


free(mpimx);
free(mpimy);
borraArreglo(mx);
borraArreglo(my);


MPI_Bcast(newcx, totalN, MPI_FLOAT, 0, MPI_COMM_WORLD);
MPI_Bcast(newcy, totalN, MPI_FLOAT, 0, MPI_COMM_WORLD);


////for(i=0;i<totalN;i++){

  ////printf("newnewcx=%f,%i,%g,%i,%i\n",newcx[i],i,timet,totalN,rank);
  ///printf("newnewcy=%f,%i,%g,%i,%i\n",newcy[i],i,timet,totalN,rank);
////  }
//free(newcx);


    mig=0;//reset mig for next time step
    step++;
    timet += dt;
  }  ///////////////////////////////////////////////////////////////////////////////////close evolution time-loop

     endtime   = MPI_Wtime();
     printf("That took %f seconds for processor %i\n",endtime-starttime,rank);



////////////////////////////////////////////////////////////////////////////////////
//Lf =  PAMR_get_max_lev(PAMR_AMRH);
    //for(l=1;l<=Lf;l++) {
      //valid = PAMR_init_s_iter(l,PAMR_AMRH,0);
      //while (valid) {
     //   ldptr();
		//initialize_(&Nx,&Ny,x,y,f[tlev1],xc,&delta);
//suma=x+y;
//	printf("res=%i,%f,%f,%f\n",rank,x,y,suma);

	//valid = PAMR_next_g();
      //}//end while
    //}//end for
    // not really necessary here, but for further reference
////    PAMR_sync(1,tlev1+1,PAMR_AMRH,0);
  }
  PAMR_free_context(pamr_context);

  
  }