/*
 * Program to create a plummer model dataset
 * for benchmarking gravitational FMM and Tree Codes
 *
 * based on Ruby code by Piet Hut and Jun Makino
 * as found in The Art of Computational Science vol 9
 *
 * C version with binary file support by Tsamouris Panayiotis
 *
 *  gcc plummer.c -std=c99 -lm
 * 
*/ 

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>

#define  IM1 2147483563
#define  IM2 2147483399
#define  AM (1.0/IM1)
#define  IMM1 (IM1-1)
#define  IA1 40014
#define  IA2 40692
#define  IQ1 53668
#define  IQ2 52774
#define  IR1 12211
#define  IR2 3791
#define  NTAB 32
#define  NDIV (1+IMM1/NTAB)
#define  EPS 1.2e-7
#define  RNMX (1.0-EPS)
#define PI 3.14159265

float ran2(long *idum)
/* Long period (> 2 × 1018) random number generator of L’Ecuyer with Bays-Durham shuffle
*  and added safeguards. Returns a uniform random deviate between 0.0 and 1.0 (exclusive of
*  the endpoint values). Call with idum a negative integer to initialize; thereafter, do not alter
*  idum between successive deviates in a sequence. RNMX should approximate the largest floating
*  value that is less than 1.
*/
{
    int j;
    long k;
    static long idum2=123456789;
    static long iy=0;
    static long iv[NTAB];
    float temp;
                                              
    if (*idum <= 0) {				 // Initialize.
       if (-(*idum) < 1) *idum=1;                //Be sure to prevent idum = 0.
       else *idum = -(*idum);     
       idum2=(*idum);     
       for (j=NTAB+7;j>=0;j--) {                 //Load the shuffle table (after 8 warm-ups).
            k=(*idum)/IQ1;
            *idum=IA1*(*idum-k*IQ1)-k*IR1;
            if (*idum < 0) *idum += IM1;
            if (j < NTAB) iv[j] = *idum;
            }                                            
 	iy=iv[0];
    }                                           
    k=(*idum)/IQ1;				//Start here when not initializing.
    *idum=IA1*(*idum-k*IQ1)-k*IR1;              //Compute idum=(IA1*idum) % IM1 without
    if (*idum < 0) *idum += IM1;		//overflows by Schrage’s method.
    k=idum2/IQ2;
    idum2=IA2*(idum2-k*IQ2)-k*IR2;
    if (idum2 < 0) idum2 += IM2;                // Compute idum2=(IA2*idum) % IM2 likewise.
    j=iy/NDIV;					// Will be in the range 0..NTAB-1.
    iy=iv[j]-idum2;				//Here idum is shuffled, idum and idum2 are
    iv[j] = *idum;				//combined to generate output.
    if (iy < 1) iy += IMM1;
    if ((temp=AM*iy) > RNMX) return RNMX;       //Because users don’t expect endpoint values.
    else return temp;
}

#define xx(x, y) x + ran2(&seed)*(y-x)

#ifndef M_PI
#define M_PI           3.14159265358979323846
#endif

// a point in the 3d space
typedef struct {
  double x, y, z;
  double vx , vy , vz;
  double charge;
} d3_point;

size_t fwrite ( const void * ptr, size_t size, size_t count, FILE * stream );
FILE *fopen(const char *filename, const char *mode);
int fclose(FILE *fp);
double sin ( double x );
double cos ( double x );
double acos ( double x );
void plummer_fill( d3_point *matr , long particles);

int main(int argc, char *argv[]){


	printf("Plummer Model Generator\n");
	
	long int particles = atoi(argv[1]);
	double total_charge=1.0;
	double side = 1.0, radius = 1.0;
	double centerx = 0.0, centery = 0.0, centerz = 0.0;

	d3_point *matrix3d;
	matrix3d = (d3_point *) malloc(particles*( sizeof( d3_point )));
	
   	plummer_fill ( matrix3d , particles);
}


void plummer_fill( d3_point *pinakas , long particles){

	int res = 0;
	long int seed=0;
	double radius = 0.0;
	double theta = 0.0;
	double velocity = 0.0;
	double phi = 0.0;
	double x =-2.0/3.0;
	double sum = 0.0;

	for (int i = 0 ; i < particles ; i++ ){
		pinakas[i].charge =ran2(&seed) ;
		sum = sum + pinakas[i].charge;
		radius = pow (ran2(&seed),x) -1;
		radius = 1.0 / sqrt (radius);
		if (radius >1.0) radius = 0.5;
		printf("%f \n", radius);
		theta = acos( xx (-1.0,1.0)) ;
		phi = xx ( 0 , 2*M_PI);
 		pinakas[i].x = radius * sin ( theta ) * cos (phi);
		pinakas[i].y = radius * sin ( theta ) * sin (phi);
		pinakas[i].z = radius * cos ( theta );
 		double x1 = 0.0;
		double y1 = 0.1;
		while ( y1 > (x1 * x1 * pow ( (1.0 - (x1 * x1)), 3.5)    )   )
		{
			x1 = xx(0,1);
			y1 = xx(0,0.1);
		}
		velocity = x1 * sqrt(2.0) * pow(( 1.0 + radius*radius), -0.25);
		theta = acos ( xx (-1, 1));
		phi = xx( 0 , 2 * M_PI);
		pinakas[i].vx = velocity * sin( theta ) * cos( phi );
		pinakas[i].vy = velocity * sin( theta ) * sin( phi );
		pinakas[i].vz = velocity * cos( theta );
	}
	
	for (int i = 0 ; i < particles ; i++ ){
		pinakas[i].charge =pinakas[i].charge/sum ;
	}
	FILE *fp;
	if ((fp = fopen("dump","w"))==NULL) {
        	printf("Cannot open file.\n");
        	exit(1);
	      }

	for (int i = 0 ; i < particles ; i++ ){				// write to file
	fprintf (fp ,"%lf %lf %lf %lf\n" , pinakas[i].x , pinakas[i].y , pinakas[i].z , pinakas[i].charge);// ,pinakas[i].vx,pinakas[i].vy,pinakas[i].vz);
	}
	
	fclose(fp);

	double *output_x;
	double *output_y;
	double *output_z;
	double *output_charge;

	output_x = malloc ( particles * (sizeof(double)));
	output_y = malloc ( particles * (sizeof(double)));
	output_z = malloc ( particles * (sizeof(double)));
	output_charge = malloc ( particles * (sizeof(double)));

	for ( int i = 0 ; i < particles ; i++){
	output_x[i] = pinakas[i].x;
	output_y[i] = pinakas[i].y;
	output_z[i] = pinakas[i].z;
	output_charge[i] = pinakas[i].charge;
	}

	 if ((fp = fopen("dump_x","w"))==NULL) {
              printf("Cannot open file.\n");
              exit(1);
	      }
	
	for (int i = 0 ; i < particles ; i++ ){				// write to file
	if (((i % 4)==0) && (i!=0)) fprintf(fp,"\n"); 	
	fprintf (fp ,"%lf " , pinakas[i].x);
	}
	
	fclose(fp);
								// open file
	   if ((fp = fopen("dump_y","w"))==NULL) {
              printf("Cannot open file.\n");
              exit(1);
	      }
	
	for (int i = 0 ; i < particles ; i++ ){				// write to file
	if (((i % 4)==0) && (i!=0)) fprintf(fp,"\n");	
	fprintf (fp ,"%lf " , pinakas[i].y);
	}
	
	fclose(fp);
								// open file
	   if ((fp = fopen("dump_z","w"))==NULL) {
              printf("Cannot open file.\n");
              exit(1);
	      }
	
	for (int i = 0 ; i < particles ; i++ ){				// write to file
	if (((i % 4)==0) && (i!=0)) fprintf(fp,"\n");
	fprintf (fp ,"%lf " , pinakas[i].z);
	}
	
	fclose(fp);
								// open file
	   if ((fp = fopen("dump_charge","w"))==NULL) {
              printf("Cannot open file.\n");
              exit(1);
	      }
	
	for (int i = 0 ; i < particles ; i++ ){				// write to file
	if (((i % 4)==0) && (i!=0)) fprintf(fp,"\n");
	fprintf (fp ,"%lf " , pinakas[i].charge);
	}
	
	fclose(fp);
	
	if ((fp = fopen("binary_x","wb"))==NULL) {
              printf("Cannot open file.\n");
             exit(1);
	    }
	fwrite( output_x, 1 , particles * sizeof(double),fp);
	fclose( fp );
	
	if ((fp = fopen("binary_y","wb"))==NULL) {
              printf("Cannot open file.\n");
             exit(1);
	    }
	fwrite( output_y, 1 , particles * sizeof(double),fp);
	fclose( fp );
	
	if ((fp = fopen("binary_z","wb"))==NULL) {
              printf("Cannot open file.\n");
             exit(1);
	    }
	fwrite( output_z , 1 , particles * sizeof(double),fp);
	fclose( fp );

	if ((fp = fopen("binary_charge","wb"))==NULL) {
              printf("Cannot open file.\n");
             exit(1);
	    }
	fwrite( output_charge , 1 , particles *sizeof(double),fp);
}































