/*
*  A program to create data dumps to be used by various codes
*  running treecodes/FMM simulations. 2D-3D particles having either
*  a charge or mass, or both, by Panayiotis Tsamouris
*  
*
*  The code utilizes the Minimal Standard Random Number Generator 
*  by Parker and Miller, as adapted by Schrage and code by 
*  the writers of "Numeric Recipes in C"
*
*   clear ; gcc 3dCubeElectric.c -lm -std=c99 -Wformat -o cube.o ; ./cube.o X 
*
*/

#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;
}


// a point in the 3d space
typedef struct {
  double x, y, z;
  double charge;
} d3_electric;

//prototypes
size_t fwrite ( const void * ptr, size_t size, size_t count, FILE * stream );
void cube_fill( d3_electric *pinakas , long int part , double mhkos , double x0, double y0, double z0 );

FILE *fopen(const char *filename, const char *mode);
int fclose(FILE *fp);
double sin ( double x );
double cos ( double x );

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

	//int res = 0;
	//res = system( "clear" );

	printf("Running 3d cube particle generator\n");
	printf("Version 1.2 by Panayiotis Tsamouris\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_electric *matrix3d;
	matrix3d = (d3_electric *) malloc(particles*( sizeof( d3_electric )));
	printf("3d matrix allocation for %ld particles\n", particles);
	cube_fill ( matrix3d , particles , side , centerx , centery, centerz );
}


void cube_fill( d3_electric *pinakas ,long int part , double mhkos , double x0, double y0, double z0){

        int res = 0;
	long int seed=0;
	double charge_sum=0.0;
	for (int i = 0 ; i < part ; i++ ){
		pinakas[i].x = x0 + mhkos *  (ran2( &seed ) - 0.5) ;
		pinakas[i].y = y0 + mhkos *  (ran2( &seed ) - 0.5) ;
		pinakas[i].z = z0 + mhkos *  (ran2( &seed ) - 0.5) ;
		//pinakas[i].charge = 1.0/part;		
		pinakas[i].charge = ran2( &seed );
		charge_sum = charge_sum + pinakas[i].charge;
	}
	for (int i = 0 ; i < part ; i++ ){
		pinakas[i].charge = pinakas[i].charge/charge_sum;
	}
	
	FILE *fp ;							// open file
	   if ((fp = fopen("dump","w"))==NULL) {
              printf("Cannot open file.\n");
              exit(1);
	      }
	
	for (int i = 0 ; i < part ; i++ ){				// write to file
	fprintf (fp ,"%E, %E, %E, %E\n" , pinakas[i].x , pinakas[i].y , pinakas[i].z , pinakas[i].charge);
	}
	
	fclose(fp);							//close file
	
        res = system("tr 'E' 'D' < dump > fortran.dump");

	/*
	if ((fp = fopen("dump","w"))==NULL) {
              printf("Cannot open file.\n");
              exit(1);
	      }
	for (int i = 0 ; i < part ; i++ ){				// write to file
	fprintf (fp ,"%lf %lf %lf %lf\n" , pinakas[i].x , pinakas[i].y , pinakas[i].z , pinakas[i].charge);
	}
	//fprintf (fp ,"%f \n" , charge_sum);
	fclose(fp);
	*/
							// open file
	 if ((fp = fopen("dump_x","w"))==NULL) {
              printf("Cannot open file.\n");
              exit(1);
	      }
	
	for (int i = 0 ; i < part ; 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 < part ; 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 < part ; 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 < part ; i++ ){				// write to file
	if (((i % 4)==0) && (i!=0)) fprintf(fp,"\n");
	fprintf (fp ,"%lf " , pinakas[i].charge);
	}
	
	fclose(fp);
	/*	
	double *p;
	p = malloc( 4 * part * (sizeof( double )));

	for ( int i = 0 ; i < part ; i++ ){
	p[i*4] = pinakas[i].x;
	p[i*4+1] = pinakas[i].y;
	p[i*4+2] = pinakas[i].z;
	p[i*4+3] = pinakas[i].charge;
	}
	
	if ((fp = fopen("binary_dump","wb"))==NULL) {
              printf("Cannot open file.\n");
             exit(1);
	    }
	fwrite( p, 1 ,4 * part *sizeof(double),fp);
	fclose( fp );
	*/

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

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

	for ( int i = 0 ; i < part ; 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("binary_x","wb"))==NULL) {
              printf("Cannot open file.\n");
             exit(1);
	    }
	fwrite( output_x, 1 , part * sizeof(double),fp);
	fclose( fp );
	
	if ((fp = fopen("binary_y","wb"))==NULL) {
              printf("Cannot open file.\n");
             exit(1);
	    }
	fwrite( output_y, 1 , part * sizeof(double),fp);
	fclose( fp );
	
	if ((fp = fopen("binary_z","wb"))==NULL) {
              printf("Cannot open file.\n");
             exit(1);
	    }
	fwrite( output_z , 1 , part * sizeof(double),fp);
	fclose( fp );

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

}

