#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "tools.h"
#include "Parameters.h"

void find_planet_particles(int *flag_planet, int ID0, double R0){
  int surface_index[np];
  int Npla = 0;
  int Nsur = 0;
  int pnum = -1;
  for(int i=0;i<np;i++) if(particles[i].pid == ID0) pnum=i;
  struct particle p0 = particles[pnum];

  for(int i=0; i<np; i++){
    struct particle p = particles[i];
    double dx = p.x - p0.x;
    double dy = p.y - p0.y;
    double dz = p.z - p0.z;
    double dr= sqrt(dx*dx + dy*dy + dz*dz);
    
    if(dr <= 0.98*R0){
      flag_planet[i] = 1;
      Npla++;
      if(dr>=0.975*R0){
	surface_index[Nsur] = i;
	Nsur++;
      }	
    }else{
      flag_planet[i] = 0;
    }
  }

  for(int i=0; i<Nsur; i++){
    for(int j=0; j<np; j++){
      if(flag_planet[j] == 1) continue;
      if( check_if_particles_are_in_touch( particles[surface_index[i]], particles[j] ) == 1){
	surface_index[Nsur] = j;
	Nsur++;
	Npla++;
	flag_planet[j] = 1;
      } 
    }
  }
}

void clump_finding(int count_c[NSIZEX], int c_index[][NSIZEX], double m_clump[NSIZEX]){
  int    N_for_clump=0;
  double m_clump_tmp=0;
  int    n_clump_tmp=0;
  int *ID_for_clump  = (int *)calloc(np+10,sizeof(int));
  int *ID_clump_tmp  = (int *)calloc(np+10,sizeof(int));
  int *flag_new      = (int *)calloc(np+10,sizeof(int));
  int *flag_done     = (int *)calloc(np+10,sizeof(int));
  int (*clump_ID_flag)[NSIZEX] = malloc(sizeof(int) * (np+10) * NSIZEX);
  for(int ii=0; ii<NSIZEX; ii++) m_clump[ii]   = 0;
  for(int ii=0; ii<NSIZEX; ii++) count_c[ii]   = 0;
  for(int ii=0; ii<np;     ii++) flag_done[ii] = 0;
  
  for(int i=0; i<np; i++){
      ID_for_clump[N_for_clump]  = i;
      N_for_clump               += 1;
  }
  N_for_clump = np; 

  for(int k=0; k<N_for_clump; k++){
    if(flag_done[ID_for_clump[k]] == 1) continue;

    ID_clump_tmp[0] = ID_for_clump[k];
    n_clump_tmp  = 0;
    m_clump_tmp  = 0;
    for(int l=0; l<N_for_clump; l++) flag_new[ID_for_clump[l]] = 0;

    for(int i=0; i<=n_clump_tmp; i++){
      for(int j=0; j<N_for_clump; j++){
	if(flag_new[ID_for_clump[j]] == 1 || flag_done[ID_for_clump[j]] == 1) continue;

	if( check_if_particles_are_in_touch(particles[ID_clump_tmp[i]], particles[ID_for_clump[j]]) == 1 ){
	  n_clump_tmp  ++;                              
	  m_clump_tmp  += particles[ID_for_clump[j]].m; 
	  flag_new[ID_for_clump[j]]  = 1;               
	  flag_done[ID_for_clump[j]] = 1;               
	  ID_clump_tmp[n_clump_tmp]  = ID_for_clump[j]; 
	}
	
      }
    }
    
    for(int ii=0; ii<NSIZEX; ii++){
      if(m_clump_tmp > m_clump[ii]){
      
	for(int jj=NSIZEX-1; jj>ii; jj--){
	  m_clump[jj] = m_clump[jj-1];
	  count_c[jj] = count_c[jj-1];
	  for(int kk=0; kk<N_for_clump; kk++) clump_ID_flag[ID_for_clump[kk]][jj] = clump_ID_flag[ID_for_clump[kk]][jj-1]; 
	}
	
	m_clump[ii] = m_clump_tmp;
	count_c[ii] = n_clump_tmp;
        for(int kk=0; kk<N_for_clump; kk++){ 
          if(flag_new[ID_for_clump[kk]]==1) clump_ID_flag[ID_for_clump[kk]][ii]=1;
          else                              clump_ID_flag[ID_for_clump[kk]][ii]=0;	
	}
	
	ii=NSIZEX; 
      }
    }

  }

  for(int i=0;i<NSIZEX;i++) count_c[i]=0;  
  for(int jj=0; jj<NSIZEX; jj++){
    for(int ii=0; ii<N_for_clump; ii++){
      if(clump_ID_flag[ID_for_clump[ii]][jj] == 1){
	c_index[count_c[jj]][jj] = ID_for_clump[ii];
	count_c[jj]++;
      }
    }
  }

  free(ID_for_clump); free(ID_clump_tmp); free(flag_new); free(flag_done); free(clump_ID_flag);
}

int check_if_particles_are_in_touch(struct particle p1, struct particle p2){
  int check=0;

  double dx = p1.x - p2.x;
  double dy = p1.y - p2.y;
  double dz = p1.z - p2.z;
  double dr = sqrt(dx*dx + dy*dy + dz*dz);

  double m_ave = 0.5*(p1.m+p2.m);         
  double pr3   = 3.0*m_ave/(4.*M_PI*rho_clump_check);
  double pr    = pow(pr3, (1./3.));       

  double dr_actual = 2.0*pr;  
  double dr_alpha  = 0;       

  if(dr <= (dr_actual + dr_alpha)) check=1;
  return check;
}

void output_clumps(int count_c[NSIZEX], int c_index[][NSIZEX], double m_clump[NSIZEX]){  
  FILE *output;
  char filename[30];

  for(int i=0; i<NSIZEX; i++){
    sprintf(filename,"clump_%04d.ascii",i);
    output=fopen(filename,"w");
    for(int j=0; j<count_c[i]; j++){
      struct particle p = particles[c_index[j][i]];
      fprintf(output,"%d\t%e\t%e\t%e\t%e\t%e\t%e\n",p.pid,p.x,p.y,p.z,p.vx,p.vy,p.vz); 
    }
    fclose(output);
  }

}

void output_mass_of_clumps(int count_c[NSIZEX], int c_index[][NSIZEX], double m_clump[NSIZEX]){
  FILE *output;
  char filename[30];
  sprintf(filename,"clump_mass_distribution.ascii");
  output=fopen(filename,"w");

  for(int i=0; i<NSIZEX; i++) fprintf(output,"N_clump[%d]= %d\tm_clump[%d]= %e\n",i, count_c[i], i, m_clump[i]);

  fclose(output);
}

void output_COM_of_clumps(int count_c[NSIZEX], int c_index[][NSIZEX], double m_clump[NSIZEX]){
  FILE *output;
  char filename[30];
  sprintf(filename,"COM_of_clumps.ascii");
  output=fopen(filename,"w");
  
  for(int i=0; i<NSIZEX; i++){
      double xc=0, yc=0, zc=0, tot_m=0;
      for(int j=0; j<count_c[i]; j++){
	struct particle p = particles[c_index[j][i]];
	xc += p.m*p.x;  yc += p.m*p.y;  zc += p.m*p.z;
	tot_m += p.m;
      xc /= tot_m;  yc /= tot_m;  zc /= tot_m;
      fprintf(output,"xc= %e\tyc= %e\tzc= %e\n",xc,yc,zc);
    }
  }

  fclose(output);
}