/**
The functions within this file aim to calculate a second-order longitudional structure function from a vector field. We can define a distance vector $\mathbf{l}$ separating two points in space ($\mathbf{x_1,x_1+l}$). Generally speaking, the vectors $\mathbf{u}$ at these points are differt. The difference in the magnitude of their projections on the $\mathbf{l}$ vector may have some interesting properties. We define:  
   
   $$\delta v_{\|}(\mathbf{x_1,l})= \left( \mathbf{v(x_1+l)-v(x_1)}\right) \cdot \frac{\mathbf{l}}{l},$$

where l is $\|\mathbf{l}\|$. By definition, from a statistical perspective, within a homogeneous and isotropic region, one would expect $\|\delta v_{\|}(\mathbf{x_1,l})\|$ only to be dependend on $l$, i.e. the length of the separation vector. Therefore we can define a second order statistic related to a vector field ($\mathbf{u}$) according to:

$$S_2(l)=\langle \delta v_{\|}(\mathbf{x_1,l})^2 \rangle$$ 

where the variable $\langle x \rangle$ represents a region-averaged value of a dummy variable $x$. The code belows aims to evaluate some discretized version of $S_2(l)$ from a vector field on a three-dimensional-tree grid.
*/

int obtain_a_cache(int j,double d[][6],int lev[],vector u){
#if _MPI
  MPI_Barrier(MPI_COMM_WORLD);
#endif
  srand(time(NULL));// In order to obtain a smooth-ish pair distances distribution we add a stogastic component;
  int index[j];
  int n=0;
  foreach()
    n++;
  double meaninc=(double)(n)/(double)(j); //The mean increment to obtain approx j cells.
  index[0]=0;
  int g=1;
  while (g<j){
    double randsteprel =((sign(noise())*pow((rand()/(double)RAND_MAX),0.2))+1.0);
    int inc = (int)((1.0*randsteprel*meaninc)+1.0);
    index[g]=index[g-1]+inc;
    g++;
  }
  n=0;
  int jj=0;
  foreach(){
    if(n==index[jj]){
      d[jj][0]=x; d[jj][1]=y; d[jj][2]=z;
      d[jj][3]=u.x[]; d[jj][4]=u.y[]; d[jj][5]=u.z[];
      lev[jj]=level;
      jj++;
    }
    n++;
  }
  printf("there are %d in cache and %d requested\n",jj,j);
  return jj; //Return jj (<=j) the number of cells in the cache
}

void set_g_and_gg(int ggg[2], int j){
  int total = j*(j-1)/2;
  int start = (pid()*total/npe())+1;
  int m = 0;
  for (int g=0;g<j;g++){
    for (int gg=g+1;gg<j;gg++){
      m++;
      if (m==start){
	ggg[0]=g;
	ggg[1]=gg;
        printf("%d %d %d %d %d %d\n",ggg[0],ggg[1],total,pid(),start,j);
	return;
      }
    }
  }
}

int structure2(double cach[][6],int lev[],int j,double L,int Len,double dvarr[][4]){
  int g=0;
  int gg=1;
#if _MPI
  int ggg[2];
  set_g_and_gg(ggg,j);
  g=ggg[0];
  gg=ggg[1];
#endif
  
  int iter=1;
  int todo=((j*(j-1))/2)/npe();
  double Del = L/(double)Len;
  double l;
  while (iter<=todo){
    l=pow(sq(cach[g][0]-cach[gg][0])+
	  sq(cach[g][1]-cach[gg][1])+
	  sq(cach[g][2]-cach[gg][2]),0.5);
    if (l>0.&&l<L0){ //distance within range
      double weight = (double)(1<<(3*(2*depth()-lev[gg]-lev[g]))); //weigh with cells' volumes
      double dv = (sq((cach[g][3]-cach[gg][3])*(cach[g][0]-cach[gg][0]))+
		   sq((cach[g][4]-cach[gg][4])*(cach[g][1]-cach[gg][1]))+
		   sq((cach[g][5]-cach[gg][5])*(cach[g][2]-cach[gg][2])))/(sq(l));
      int index = (int)((l/Del)+0.5);
      if  (index<Len){
	dvarr[index][1]+=dv*weight;
	dvarr[index][2]+=weight;
	dvarr[index][3]+=1.;
      }
    }
    gg++; 
    if (gg==j){
      g++;
      gg=g+1;
    }
    iter++;
  }
  return 1;
}
void printdvarr(FILE * fp, double dvarr[][4], int Len){
  int it = 0;
  while (it<Len){
    for (int gg=0;gg<4;gg++)
      fprintf(fp,"%g\t",dvarr[it][gg]);
    fprintf(fp,"\n");
    it++;
  }
}

#if _MPI
int combine_MPI_caches(double cach[][6],int lev[],int j){
  int tot=0;
  MPI_Allreduce(&j,&tot,1,MPI_INT,MPI_SUM,MPI_COMM_WORLD);
  double cachsend[j*6];
  double * cachrecv= (double *) malloc(sizeof(double)* tot * 6);
  int  levrecv[tot];
  for (int g=0;g<j;g++){
    for (int gg=0;gg<6;gg++)
      cachsend[6*g+gg]=cach[g][gg];
  }
  MPI_Allgather(&cachsend[0],6*j,MPI_DOUBLE,cachrecv,6*j,MPI_DOUBLE,MPI_COMM_WORLD);
  MPI_Allgather(&lev[0],j,MPI_INT,&levrecv[0],j,MPI_INT,MPI_COMM_WORLD);
  for (int g=0;g<tot;g++){
    lev[g]=levrecv[g];
    for (int gg=0;gg<6;gg++)
      cach[g][gg]=cachrecv[6*g+gg];
  }
  return tot;
  
}
#endif
//The user interface function; 

int long2structure(FILE * fp,vector u,int n,double L,int Len){
  double Del = L/(double)Len;
  double dvarr[Len][4];
  int it=0;
  while (it<Len){
    for (int gg=0;gg<4;gg++)
      dvarr[it][gg]=0.;
    it++;
  }
  
  double cach[n][6];
  int lev[n];
  int locpnts = ((n/npe())+0.5);
  int j=obtain_a_cache(locpnts,cach,lev,u);
#if _MPI
  int d=j;
  j=combine_MPI_caches(cach,lev,d);
#endif
  structure2(cach,lev,j,L,Len,dvarr);
  it=0;
#if _MPI
  double dvarrrecv[Len][4]; // An array to store the reduced version of *dvarr*, according to MPI_SUM
  MPI_Reduce(&dvarr[0][0],&dvarrrecv[0][0],Len*4,MPI_DOUBLE,MPI_SUM,0,MPI_COMM_WORLD);
  if (pid()==0){
    while (it<Len){
      for (int gg=0;gg<4;gg++)
	dvarr[it][gg]=dvarrrecv[it][gg];
      it++;
    }
    it=0;
#endif
    double aa;
    while (it<Len){
      if (dvarr[it][2]>0.){
	aa = dvarr[it][1]/dvarr[it][2];
	dvarr[it][1] = aa;
      }
      dvarr[it][0]=it*Del;
      it++;
    }
    printdvarr(fp,dvarr,Len);
    
#if _MPI
  }
#endif
  return 1; 
}