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


size_t fwrite ( const void * ptr, size_t size, size_t count, FILE * stream );
size_t fread(void *buf, size_t size, size_t count, FILE *stream);
long int ftell(FILE *stream);

FILE *fopen(const char *filename, const char *mode);
int fclose(FILE *fp);
double sqrt(double num); 
double fabs(double num); 
double pow(double base, double exp);
int myComparisonFunction(const void *x, const void *y);

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

	int res;
	size_t result;
	FILE * fp;
  	long int file_size_field;
	long int file_size_potential;

	double * final_field;
	double * direct_field;

	double * final_potential;
	double * direct_potential;



	//res = system("echo Running binary sdiff on points in space");
	//res = system("sdiff -sb points_in_space points_in_space1");
	//res = system("sdiff -sb points_in_space points_in_space2");
//============================================================================//
//									      //
//			FINAL FIELD VALUES				      //
//									      //
//============================================================================//
	if ((fp = fopen("final_field","rb"))==NULL) {
        printf("Cannot open file.\n");
        exit(1);
  	}

	// obtain file size:
  	fseek ( fp , 0 , SEEK_END );
  	file_size_field = ftell ( fp );
  	rewind ( fp );

  	printf("%ld bytes in file final_field\n",file_size_field );
	//printf("%ld sizeof double\n", sizeof(double));
  
  	/* allocate memory to contain the whole file:*/
  	if (( final_field =(double *) malloc(file_size_field*sizeof(double))) == NULL) {
	printf("Cannot allocate final_field.\n");
        exit(1);
 	}

  	/* copy the file into the buffer: */
  	if ((result=fread(final_field,sizeof(double),file_size_field/sizeof(double),fp)) !=(size_t) (file_size_field/sizeof(double)) ) {
  	printf("Read Error fread crash final_field(x,y,z)\n");
        exit(1);
  	}
	long int particles = result;
	//printf("%ld == particles\n\n\n\n\n\n\n\n\n\n",particles);
  	fclose ( fp );
  	//printf("%ld doubles read from final_field\n",(long)result);
	
	
//============================================================================//
//									      //
//			DIRECT FIELD VALUES				      //
//									      //
//============================================================================//

	if ((fp = fopen("direct_field","rb"))==NULL) {
        printf("Cannot open file direct_field.\n");
        exit(1);
  	}

	// obtain file size:
  	fseek ( fp , 0 , SEEK_END );
  	file_size_field = ftell ( fp );
  	rewind ( fp );

  	printf("%ld bytes in file direct_field\n",file_size_field );
//  	printf("%ld sizeof double\n", sizeof(double));
  
  	/* allocate memory to contain the whole file:*/
  	if (( direct_field =(double *) malloc( file_size_field*sizeof(double))) == NULL) {
	printf("Cannot allocate direct field.\n");
        exit(1);
  	}

  	/* copy the file into the buffer: */
  	if ((result=fread(direct_field,sizeof(double),file_size_field/sizeof(double),fp)) !=(size_t) (file_size_field/sizeof(double)) ) {
  	printf("Read Error fread crash in direct_field.\n");
        exit(1);
  	}
  	fclose ( fp );
  	//printf("%ld doubles read from direct_field\n",(long)result);
	
//============================================================================//
//									      //
//			FINAL POTENTIAL					      //
//									      //
//============================================================================//
	if ((fp = fopen("final_potential","rb"))==NULL) {
        printf("Cannot open file final_potential.\n");
        exit(1);
  	}

	// obtain file size:
  	fseek ( fp , 0 , SEEK_END );
  	file_size_potential = ftell ( fp );
  	rewind ( fp );

  	printf("%ld bytes in file final_potential\n",file_size_potential );
//  	printf("%ld sizeof double\n", sizeof(double));
  
  	/* allocate memory to contain the whole file:*/
  	if (( final_potential =(double *) malloc( file_size_potential*sizeof(double))) == NULL) {
	printf("Cannot allocate final_potential file.\n");
        exit(1);
  	}

  	/* copy the file into the buffer: */
  	if ((result=fread(final_potential,sizeof(double),file_size_potential/sizeof(double),fp)) != (size_t) (file_size_potential / sizeof(double)) ) {
  	printf("Read Error fread crash.\n");
        exit(1);
  	}
  	fclose ( fp );
  	//printf("%ld doubles read from final_potential\n",(long)result);

//============================================================================//
//									      //
//			DIRECT POTENTIAL				      //
//									      //
//============================================================================//

	if ((fp = fopen("direct_potential","rb"))==NULL) {
        printf("Cannot open file direct_potential.\n");
        exit(1);
  	}

	// obtain file size:
  	fseek ( fp , 0 , SEEK_END );
  	file_size_potential = ftell ( fp );
  	rewind ( fp );

  	printf("%ld bytes in file direct_potential\n",file_size_potential );
//  	printf("%ld sizeof double\n", sizeof(double));
  
  	/* allocate memory to contain the whole file:*/
  	if (( direct_potential =(double *) malloc( file_size_potential*sizeof(double))) == NULL) {
	printf("Cannot allocate direct_potential.\n");
        exit(1);
  	}

  	/* copy the file into the buffer: */
  	if ((result=fread(direct_potential,sizeof(double),file_size_potential/sizeof(double),fp)) !=(size_t) (file_size_potential/sizeof(double)) ) {
  	printf("Read Error fread crash.\n");
        exit(1);
  	}
  	fclose ( fp );
  	//printf("%ld doubles read from direct_potential\n",(long)result);
//============================================================================//
//									      //
//			POTENTIAL COMPARISON				      //
//									      //
//============================================================================//
	printf("         potential size  = %ld  field size = %ld  \n\n\n",file_size_potential,file_size_field);

	
	double * diff_final_direct;
	double sum = 0.0;
	double mean_abs_diff = 0.0;
	double max_diff = 0.0;
	
  	if (( diff_final_direct =(double *) malloc( file_size_potential*sizeof(double))) == NULL) {
	printf("Cannot allocate diff matrix.\n");
        exit(1);
  	}
	
	for (int i = 0 ; i <= file_size_potential ; i++){
	diff_final_direct[i] = fabs( final_potential[i] - direct_potential[i]); // abs diff
	if ( diff_final_direct[i] > max_diff) max_diff= diff_final_direct[i];
	sum = sum + diff_final_direct[i]; //sum abs diff
	}
	
	mean_abs_diff = (sum / (double) file_size_potential );
	
	printf("max abs difference between FMM and DIRECT POTENTIAL is %E\n" , max_diff);
	printf("mean abs difference between FMM and DIRECT POTENTIAL is %E\n",mean_abs_diff);
	
	


//============================================================================//
//									      //
//			FIELD COMPARISON 				      //
//									      //
//============================================================================//
	
	double * diff_final_direct_f;
	double * diff_final_direct_fx;
	double * diff_final_direct_fy;
	double * diff_final_direct_fz;

  	if (( diff_final_direct_f =(double *) malloc(file_size_field*sizeof(double))) == NULL) {
	printf("Cannot allocate diff matrix.\n");
        exit(1);
  	}	if (( diff_final_direct_fx =(double *) malloc(file_size_potential*sizeof(double))) == NULL) {
	printf("Cannot allocate diff matrix.\n");
        exit(1);
  	}
	if (( diff_final_direct_fy =(double *) malloc(file_size_potential*sizeof(double))) == NULL) {
	printf("Cannot allocate diff matrix.\n");
        exit(1);
  	}
	if (( diff_final_direct_fz =(double *) malloc(file_size_potential*sizeof(double))) == NULL) {
	printf("Cannot allocate diff matrix.\n");
        exit(1);
  	}
	double sum_f = 0.0;
	double sumx = 0.0;
	double sumy = 0.0;
	double sumz = 0.0;

	double mean_abs_diff_f = 0.0;
	double mean_abs_diff_fx = 0.0;
	double mean_abs_diff_fy = 0.0;
	double mean_abs_diff_fz = 0.0;

	double max_diff_f = 0.0;
	double max_diff_fx = 0.0;
	double max_diff_fy = 0.0;
	double max_diff_fz = 0.0;
		
	
	double metro_final = 0.0;
	double metro_direct = 0.0;

	for (int i = 0 ; i < (particles/3) ; i++){
	
	diff_final_direct_fx[i] = fabs( final_field[i*3] - direct_field[i*3]); // abs diff x
	if (diff_final_direct_fx[i] > max_diff_fx) max_diff_fx = diff_final_direct_fx[i];
	sumx = sumx + diff_final_direct_fx[i]; //sum abs diff
	//printf( " particle %d final_x %lf direct_x %lf diff_x %lf \n",i,final_field[i*3],direct_field[i*3],diff_final_direct_fx[i]);
	
	
	diff_final_direct_fy[i] = fabs( final_field [(i*3)+1] - direct_field[(i*3)+1]); // abs diff y
	if (diff_final_direct_fy[i] > max_diff_fy) max_diff_fy = diff_final_direct_fy[i];
	sumy = sumy + diff_final_direct_fy[i];
	//printf( " particle %d final_y %lf direct_y %lf diff_y %lf \n",i,final_field [(i*3)+1],direct_field[(i*3)+1],diff_final_direct_fy[i]);	
	
	
	diff_final_direct_fz[i] = fabs( final_field[(i*3)+2] - direct_field[(i*3)+2]); // abs diff x
	if (diff_final_direct_fz[i] > max_diff_fz) max_diff_fz = diff_final_direct_fz[i];
	sumz = sumz + diff_final_direct_fz[i];
	//printf( " particle %d final_z %lf direct_z %lf diff_z %lf \n",i,final_field [(i*3)+2],direct_field[(i*3)+2],diff_final_direct_fz[i]);	

	
	metro_final = sqrt(pow(final_field[i*3],2)+pow(final_field[(i*3)+1],2)+pow(final_field[(i*3)+2],2));
	metro_direct = sqrt(pow(direct_field[i*3],2)+pow(direct_field[(i*3)+1],2)+pow(direct_field[(i*3)+2],2));

	diff_final_direct_f[i] = fabs( metro_final - metro_direct );
	if (diff_final_direct_f[i] > max_diff_f){
	//printf("			 i is %d 	and %lf		\n\n\n\n",i,diff_final_direct_f[i]);	
	max_diff_f = diff_final_direct_f[i];
	//printf("  @@@@@ %lf \n\n\n\n",max_diff_f);
	}
	sum_f = sum_f + diff_final_direct_f[i];	
	
	//printf(" particle %d diff_x is %E diff_y is %E diff_z is %E \n", i ,diff_final_direct_fx[i] , diff_final_direct_fy[i], diff_final_direct_fz[i]);
	
	}


//	for (int i = 0 ; i < particles ; i++){
//	printf(" %lf %lf %lf %lf\n",metro_final,final_field[i*3],final_field [(i*3)+1],final_field[(i*3)+2]);
//	}
	
//	for (int i = 0 ; i < particles ; i++){
//	printf(" %lf %lf %lf %lf\n",metro_direct,direct_field[i*3],direct_field [(i*3)+1],direct_field[(i*3)+2]);
//	}

	mean_abs_diff_f =  (sum_f / (double) (file_size_field/3) );
	mean_abs_diff_fx = (sumx / (double) (file_size_field/3) );
	mean_abs_diff_fy = (sumy / (double) (file_size_field/3) );
	mean_abs_diff_fz = (sumz / (double) (file_size_field/3) );

	
	printf("max abs difference between FMM and DIRECT field X is %E\n" , max_diff_fx);
	printf("max abs difference between FMM and DIRECT field Y is %E\n" , max_diff_fy);
	printf("max abs difference between FMM and DIRECT field Z is %E\n" , max_diff_fz);
	printf("max abs difference between FMM and DIRECT field is %E \n" , max_diff_f);

	printf("mean abs difference between FMM and DIRECT field X is %E\n",mean_abs_diff_fx);
	printf("mean abs difference between FMM and DIRECT field Y is %E\n",mean_abs_diff_fy);
	printf("mean abs difference between FMM and DIRECT field Z is %E\n",mean_abs_diff_fz);
	printf("mean abs difference between FMM and DIRECT field is %E\n",mean_abs_diff_f);


//============================================================================//
//									      //
//			FREE						      //
//									      //
//============================================================================//

	free(final_field);
	free(direct_field);

	free(final_potential);
	free(direct_potential);

	free(diff_final_direct);
	
	free(diff_final_direct_f);
	free(diff_final_direct_fx);
	free(diff_final_direct_fy);
	free(diff_final_direct_fz);

}



