/*
Create a Sum of Gaussian (SOG) charge distribution from an input file
and print the data in two output files, one in only two columns with radius
and charge distribution, and a second with radius and each Atom's charge
distribution in its own column.
*/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#define PI 3.14159265359
#define e 1.60217646e-19

int main(){

    //Program Introduction
    printf("\n\nThis program takes data from the file SOG_data.dat\n"); 
    printf("and calculates the charge density using the Sum of Gaussian approach.\n"); 
    printf("The charge density can be produced using standard calculation\n");
    printf("or it can be normalized by multiplying the density by A/Z (p*A/Z).\n");  
    printf("Once calculated graphs are created of the density v. radius.\n");
    printf("Finally, if desired, the program will also take the density calculations\n");
    printf("and produce <r^2>=4*pi/Z*int(r^4*p)dr in order to calculate\n"); 
    printf("the nuclear density given by n=(3*A)/(4*Pi*R^3) where R=sqrt(5*<r^2>/3).\n"); 
    printf("A graph will then be created showing n v. A.  At the end of the program\n"); 
    printf("a list of the files created will be shown on the screen.\n\n");
    printf("Hit enter to continue with the program\n\n");
    getchar();

    FILE *count=fopen("SOG_data.dat", "r");
    int atomcount;
    char c;
  
    atomcount=0;
    if (count==NULL) perror("Error opening file");
    else {
        do{                                                     //counts number of atoms (lines) in
            c=fgetc (count);                                      //file in order to seed array that
            if (c == '\n') atomcount++;                           //contains results
        } while (c != EOF);
        fclose(count);
    }

    float p, A_i, r, gamma, rplus, rminus, rms, RP;
    char atom[10], label[10];
    int status, Z, A;
    int a, n, i, j; 
    int counter, option, pick, ND;
    float R[12], Q[12], results[1000][atomcount+1];
    int loopcounter[atomcount], marker;
    float pmax=0.0, rmax=0.0, zero=0.0, az;
    float sum, integral, r2, Ratom, Re, densityR, densityRe;
 
    FILE *inp=fopen("SOG_data.dat", "r");                         //Open input and output
    FILE *outp1=fopen("SOG_chargedistribution.dat", "w");         //files for reading and
    FILE *outp2=fopen("SOG_chargedistribution_columns.dat", "w"); //writing
    FILE *outp3=fopen("titles.tmp", "w");
    FILE *outp4=fopen("SOG_chargedistribution_fence.dat", "w");
    FILE *outp5=fopen("SOG_NuclearDensity.dat", "w");
    FILE *outp6=fopen("Files.txt", "w");

    if (inp==NULL){
        printf("Can't Find File\n");                          //if file doesn't exist
        return 0;                                             //exit program
    }
  
    fprintf(outp6, "SOG_chargedistribution.dat\n");
    fprintf(outp6, "SOG_chargedistribution_columns.dat\n");
    fprintf(outp6, "SOG_chargedistribution_fence.dat\n");
    fprintf(outp6, "SOG_NuclearDensity.dat\n");  
  
    printf("Do you want a graph of the Nuclear Density?\n1-Yes 2-No\n");
    while(1){
        scanf("%d", &ND);
        if (ND==1 || ND==2) break;
        else printf("Not a valid choice, choose Yes (Enter 1) or No (Enter 2)\n");
    }
         
    printf("Would you like to plot charge density or adjusted charge density (p*Z/A)?\n");
    printf("1-Charge Density 2-Adjusted Charge Density\n");
    while(1){  	
  	    scanf("%d", &pick);
  	    if (pick==1 || pick==2) break;
  	    else printf("Not a valid choice choose Charge Density (type 1) or Adjusted Charge Density (type 2)\n");
    }
                                                          //Skips first line of file
    while (c != '\n') c=fgetc(inp);                         //that only contains column titles

    fprintf(outp2, "Radius (fm)");  
    counter=1;                                              //Seed integer for charge dist. array

    while (1){
        status=fscanf(inp, "%d %d %s %f %f", &Z, &A, atom, &rms, &RP);
        if (status==EOF) break;                            //while file can be read, read in atomic number
        r=0.0;                                             //atom name, RMS, and RP.  Break at EOF
        marker=0;
        sum=0.0;
        gamma=sqrt(2./3.)*RP;                              //set r and p to zero for each loop, resetting eq.

        fprintf(outp1, "\n\n%d%s\n", A, atom);             //headers for output files
        fprintf(outp2, "----- %d%s ----", A, atom);
        fprintf(outp3, "%d%s\n", A, atom);                 //print atom names to file for labeling in graphs
        printf("%d - %d%s\n", counter, A, atom);           //Print index of atoms in file, used for selecting
                                                         //what is to be graphed
                                                         
        for (n=0; n<12; n++){                              //Loops over the twelve R's and Q's in data file
    	    fscanf(inp, "%f %f", &R[n], &Q[n]);            //and places them in array for equation
	    }
	   
	    if (pick==2) az=(float)A/Z;
	    else az=1;
	  
        for (i=0; i<1000; i++){                           //Loop for each radius, 0.01 intervals
            A_i=0.0;                                        //Reset both A and p to zero for each new radius
    	    p=0.0;
    	    for (n=0; n<12; n++){                           //loop for each R and Q, making sure that p=sum of
                rplus=pow((r+R[n])/gamma, 2.0);               //p's for each R and Q used.
    		    rminus=pow((r-R[n])/gamma, 2.0);
    		    A_i=(Q[n]*Z)/((2.0*pow(PI,1.5)*pow(gamma,3))*(1+2*pow(R[n]/gamma,2)));
    		    p=(A_i*((exp(-rminus))+exp(-rplus)))+p;
    	    }
    	    results[i][0]=r;                                   //placing r and density into results array
    	    r=r+0.01;                                       //print results to SOG_chargedistribution.dat
    	    if (p>0.0001){
    		    marker++;
    		    results[i][counter]=p*az;                          //setting charge density per unit charge
    		    fprintf(outp1, "%f %g\n", results[i][0], results[i][counter]);
    		}
            if (ND==1){
    	  	    integral=r*r*r*r*p*0.01;
    	  	    sum+=integral;
    	  	}
        }                                                //adds 0.01 to r and begin another loop
        
        loopcounter[counter-1]=marker;
        if (ND==1){
            r2=pow(4*PI/Z*sum, 0.5);
            Re=pow(5*r2*r2/3, 0.5);
            densityRe=(3.*A)/(4.*PI*Re*Re*Re);
            fprintf(outp5, "%d %lf %lf %lf\n", A, densityRe, r2, rms);
        }
        counter++;                                        //places next density results in next columns of array
    }

    fprintf(outp2, "\n\n");
    for (n=0; n<1000; n++){                               //Loop for printing results into columns
        fprintf(outp2, "%-.3f        ", results[n][0]);     //Can be used for printing individual atoms
        for(i=1; i<counter; i++){                           //to GNUplot for graphing
    	    if(loopcounter[i-1]>n) fprintf(outp2, "%-15.6g", results[n][i]);
    	    else fprintf(outp2, "%-15.6g", 0.0/zero);
        }
        fprintf(outp2, "\n");
    }
    
    fclose(outp5);
    printf("%d - all\n", counter);                        //option for printing all atoms on one graph
    printf("%d - done\n", counter+1);
    printf("Which elements do you want an individual graph of?\n");
    printf("(You will be able to group atoms later in the program)\n");
    printf("Type index number (1, 2, 3...) and hit return.\n");
    printf("do not type Name of atom!, type %d when completed\n", counter+1);

    int choice[counter];                                 //array for scanning in indexes for graphs
    a=0;                                                 //counter variable used to place info into choice[]

    while (option != counter+1){                         //Loop to scan in indices from terminal in order
        scanf("%d", &option);                              //to choose which atoms to graph
        if (option > counter+1 || option < 1)
        printf("Not an available option, pick again\n");
        else if (option==counter+1);
        else {
            choice[a]=option; 
            a++;
        }
    }

    fclose(inp);
    fclose(outp1);                                        //Close input and output files
    fclose(outp2);
    fclose(outp3);
  
    if (ND==1){
        FILE *GNUoutp=fopen("SOG_NuclearDensity.gnu", "w");
        fprintf(outp6, "SOG_NuclearDensity.gnu\n");
        fprintf(GNUoutp, "set term pdf\nset output \"SOG_NuclearDensity.pdf\"\n");
        fprintf(outp6, "SOG_NuclearDensity.pdf\n");
        fprintf(GNUoutp, "set title \"Nuclear Density v. A\"\n");
        fprintf(GNUoutp, "set ylabel \"fm^-3\"\n set xlabel \"A\"\n");
        fprintf(GNUoutp, "plot \"SOG_NuclearDensity.dat\" u 1:2 ti \"\"\n");
        fclose(GNUoutp);
        system ("gnuplot < SOG_NuclearDensity.gnu");
    }

    for(n=0; n<a; n++){
        FILE *inp2=fopen("titles.tmp", "r");                //open file each loop in order to
        if (choice[n]==counter){                            //read from the beginning each time
            FILE *GNUoutp=fopen("SOG_GNUoutp_All.gnu", "w");
            if (pick==2) {
      	        fprintf(GNUoutp, "set term pdf\nset output \"SOG_AZ.pdf\"\n");
      	        fprintf(outp6, "SOG_AZ.pdf\n");
      	    }
            else {
      	        fprintf(GNUoutp, "set term pdf\nset output \"SOG.pdf\"\n");
      	        fprintf(outp6, "SOG.pdf\n");
      	    }
            fprintf(GNUoutp, "set title \"Nuclear Charge Density (SOG)\"\n");
            fprintf(GNUoutp, "set xlabel \"Radius (fm)\"\n");
            fprintf(GNUoutp, "set ylabel \"Charge Density (e*fm^-3)\"\n");
            
            if (pick==1) fprintf(GNUoutp, "set yrange [0:0.12]\nset xrange [0:10]\n");
            else fprintf(GNUoutp, "set yrange [0:0.24]\nset xrange [0:10]\n");
            
            fprintf(GNUoutp, "plot \"SOG_chargedistribution.dat\" using 1:2 with lines lw 2 ");
            fprintf(GNUoutp, "title \"\"\n");                 //GNUplot formatting
            fprintf(GNUoutp, "unset output\n");
            fclose(GNUoutp);
            system("gnuplot < SOG_GNUoutp_All.gnu");
        }
        
        else {
            FILE *GNUoutp=fopen("SOG_GNUoutp.gnu", "w");
            fprintf(outp6, "SOG_GNUoutp.gnu\n");
            FILE *GNUcharge=fopen("charge.tmp", "w");
            fprintf(outp6, "SOG_GNUoutp.gnu\n");
            for (j=0; j<loopcounter[choice[n]-1]; j++){                           //Loop that prints data for select atoms
    	        fprintf(GNUcharge, "%-.3f %-15.6g\n", results[j][0], results[j][choice[n]]);
            }

            for(i=0; i<choice[n]; i++) {                      //Loop that reads data from file and 
    	        fscanf(inp2, "%s", label);                      //saves name of atom for labeling on graph
            }

            if (pick==2) {
          	    fprintf(GNUoutp, "set term pdf\nset output \"%s_SOG_AZ.pdf\"\n", label);
          	    fprintf(outp6, "%s_SOG_AZ.pdf\n", label);
          	}

            else {
          	    fprintf(GNUoutp, "set term pdf\nset output \"%s_SOG.pdf\"\n", label);
          	    fprintf(outp6, "%s_SOG.pdf\n", label);
          	}

            fprintf(GNUoutp, "set title \"Nuclear Charge Density (SOG) for %s\"\n", label);
            fprintf(GNUoutp, "set xlabel \"Radius (fm)\"\n");
            fprintf(GNUoutp, "set ylabel \"Charge Density (e*fm^-3)\"\n");
            if (pick==1) fprintf(GNUoutp, "set yrange [0:0.12]\nset xrange [0:10]\n");
            else fprintf(GNUoutp, "set yrange [0:0.24]\nset xrange [0:10]\n");
            fprintf(GNUoutp, "plot \"charge.tmp\" using 1:2 with lines lw 4 ");
            fprintf(GNUoutp, "title \"\"\n");                 //GNUplot formatting
            fprintf(GNUoutp, "unset output\n");
            fclose(GNUoutp);
            fclose(GNUcharge);
            system("gnuplot < SOG_GNUoutp.gnu");
            system("rm charge.tmp"); 
        }
        fclose(inp2);
    }  

    FILE *inp2=fopen("titles.tmp", "r");

    for (n=1; n<counter; n++){                             //Loop to print indices and atom names
        fscanf(inp2, "%s", label);                           //to screen for selecting
        printf("%d - %s\n", n, label);
    }
    printf("%d - done\n", counter);
    printf("Which elements do you want to graph together?\n");
    fclose(inp2);

    option=counter-1;                                      //sets option !=counter so while loop
    a=0;                                                   //works, and reset a=0 for counting

    while (option != counter){                             //Loop to scan indices so graph only
        scanf("%d", &option);                                //graphs specific atoms
        if (option > counter || option < 1)
            printf("Not an available option, pick again\n");
        else if (option==counter);
        else {
            choice[a]=option+1;                                //add 1 because when graph column 1 is
            a++;                                               //radius, graph 1:choice[] is to plot
        }                                                    //multiple selected atoms on 1 graph
    }

    if(a!=0){
        int graphtype;
        printf("Now that you have picked your atoms would you like them in a 3D fenceplot or a normal 2D graph?\n");
        printf("2D=1, 3D=2, Both=3\n");

        while (1){                         //Loop to scan in indices from terminal in order
            scanf("%d", &graphtype);                              //to choose which atoms to graph
            if (graphtype > 3 || graphtype < 1){
                printf("Not an available option, pick again\n");
                printf("2D=1, 3D=2, Both=3\n");
            }
            else break;
        }

        if(graphtype==1 || graphtype==3){
                                                    //Creates graph of multiple atoms
            FILE *GNUoutp=fopen("SOG_GNUoutp_Group.gnu", "w");
            fprintf(outp6, "SOG_GNUoutp_Group.gnu\n");
    
            if (pick==2) {
        	    fprintf(GNUoutp, "set term pdf\nset output \"SOG_Group_AZ.pdf\"\n");
        	    fprintf(outp6, "SOG_Group_AZ.pdf\n");
            }
            else {
        	    fprintf(GNUoutp, "set term pdf\nset output \"SOG_Group.pdf\"\n");
        	    fprintf(outp6, "SOG_Group.pdf\n");
            }
            fprintf(GNUoutp, "set key outside\n");
            fprintf(GNUoutp, "set title \"Nuclear Charge Density (SOG)\"\n");
            fprintf(GNUoutp, "set xlabel \"Radius (fm)\"\n");
            fprintf(GNUoutp, "set ylabel \"Charge Density (e*fm^-3)\"\n");
        
            if (pick==2) fprintf(GNUoutp, "set yrange[0:0.24]\nset xrange [0:10]\n");
            else fprintf(GNUoutp, "set yrange [0:0.12]\n set xrange [0:10]\n");
        
            fprintf(GNUoutp, "plot ");
            for(n=0; n<a; n++){
                FILE *inp2=fopen("titles.tmp", "r");                //Loop to read titles.tmp in order 
                for (i=1; i<choice[n]; i++){                        //to label each line, must open and
	                fscanf(inp2, "%s", label);                        //close file to read from start of
                }                                                   //file each time
                fclose(inp2);
                fprintf(GNUoutp, "\"SOG_chargedistribution_columns.dat\" ");
                fprintf(GNUoutp, "using 1:%d with lines lw 4 lc %d ", choice[n], n+1);
                fprintf(GNUoutp, "title \"%s\" ", label);          //use previously created file to plot
                if (n==a-1) fprintf(GNUoutp, "\n");                //each individual graph by using choice[]
                else fprintf(GNUoutp, ", ");                       //to state which columnes to use
            }
            fprintf(GNUoutp, "\n");
            fprintf(GNUoutp, "unset output\n");
            fclose(GNUoutp);
            system("gnuplot < SOG_GNUoutp_Group.gnu");
        }
    
        if(graphtype==2 || graphtype==3) {
        	for(n=0; n<a; n++){
	            for(i=0; i<1000; i++){
	                if (pmax<results[i][choice[n]-1]) pmax=results[i][choice[n]-1];
	                if (rmax<results[i][0]) rmax=results[i][0];
	                if (results[i][choice[n]-1]<0.0001) i=1000;
	                else fprintf(outp4, "%d %f %f\n", n+1, results[i][0], results [i][choice[n]-1]);
	            }
	            fprintf(outp4, "\n\n");
	        }
	    
	        fclose(outp4);
	        FILE *GNUoutp=fopen("SOG_GNUoutp_Fence.gnu", "w");
	        fprintf(outp6, "SOG_GNUoutp_Fence.gnu\n");

	        fprintf(GNUoutp, "set term pdf\n");
	        if (pick==2) {
	      	    fprintf(GNUoutp, "set output \"SOG_Fence_AZ.pdf\"\n");
	      	    fprintf(outp6, "SOG_Fence_AZ.pdf\n");
	      	}
	        else {
	      	    fprintf(GNUoutp, "set output \"SOG_fence.pdf\"\n");
	      	    fprintf(outp6, "SOG_fence.pdf\n");
	      	}
	  	
            fprintf(GNUoutp, "set key outside\n");
            fprintf(GNUoutp, "set title \"Nuclear Charge Density (SOG)\"\n");
            fprintf(GNUoutp, "set xlabel \"Radius (fm)\"\n");
            fprintf(GNUoutp, "set zlabel \"Charge Density\\n\(e*fm^-3)\" offset -5\n");
            fprintf(GNUoutp, "set xrange[0:%lf] #range for radius\n", rmax);
            fprintf(GNUoutp, "set yrange[0:%d] #range for counter\n", a);
            fprintf(GNUoutp, "set zrange[0:%lf] #range for charge density\n", pmax);
            fprintf(GNUoutp, "set view 72, 355\n");
            fprintf(GNUoutp, "set xyplane 0\n");
            fprintf(GNUoutp, "set style data filledcurves x1\n");
            fprintf(GNUoutp, "splot ");
	        for(n=a-1; n>-1; n--){
	            FILE *inp2=fopen("titles.tmp", "r");
	            for(i=1; i<choice[n]; i++){
	                fscanf(inp2, "%s", label);
	            }
	            fclose(inp2);
	            fprintf(GNUoutp, "\"SOG_chargedistribution_fence.dat\" u ($1==%d?($2):1/0):($1):3 ti '%s' lw 2 lc %d", n+1, label, n+1);
	            if (n==0) fprintf(GNUoutp, "\n");
	            else fprintf(GNUoutp, ", ");
	        }
	        fclose(GNUoutp);
	        system("gnuplot < SOG_GNUoutp_Fence.gnu");
        }
    }
    fclose(outp6);
    system("rm titles.tmp");                               //remove temporary file
    printf("The program has created the following files:\n");  
    FILE *Files=fopen("Files.txt", "r");
    char file[35];  
    while(1){
        status=fscanf(Files, "%s", file);
        if (status==EOF) break;
        printf("%s\n", file);
    }
    return 0;
}