/*program to calculate the charge density
using FB parameters spherical Bessel functions
http://en.wikipedia.org/wiki/Bessel_function.
*/

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

int main(){

    FILE *count=fopen("FB_data.dat", "r");
    int atomcount, s;
    char c;
                                                          //Skips first line of file
    while (c != '\n') c=fgetc(count);                       //that only contains column titles
  
    atomcount=0;
    do{
        s=fgetc(count);                                       //counts number of atoms (lines) in
        c=fgetc(count);                                       //file in order to seed array that
        if (s == '\n') atomcount++;                           //contains results
    } while (c != EOF);
  
    fclose(count);

    float x, rn[17], r, RR, term;
    char label[10];
    int counter, Z, A;
    int n, i, j, a;
    int status, pick;
    int loopcounter[atomcount], marker;
    float  results[1000][atomcount];
    float pmax=0, rmax=0, az, rfm;
    float sum, integral, r2, Ratom, Re, densityR, densityRe;
 
    FILE *inp=fopen("FB_data.dat", "r");           //Open input and output files
    FILE *outp1=fopen("FB_chargedistribution.dat", "w");
    FILE *outp2=fopen("FB_chargedistribution_columns.dat", "w");
    FILE *outp3=fopen("titles.tmp", "w");
    FILE *outp4=fopen("FB_readinwriteout.txt", "w");
    FILE *outp5=fopen("FB_chargedistribution_fence.dat", "w");
    FILE *outp6=fopen("FB_nucleardensity.dat", "w");

    if (inp==NULL){
        printf("Can't Find File\n");                          //if file doesn't exist
        return 0;                                             //exit program
    }
  
    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");
    }
  
    while (c != '\n') c=fgetc(inp);

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

    while (1){
        status=fscanf(inp, "%s %d %d", label, &A, &Z);
        if (status==EOF) break;                            //while file can be read, read in atomic label
        fprintf(outp4, "%s %d %d ", label, A, Z);          //(A-atom), A, and Z.  Break at EOF
        r=0.000000001;     
        marker=0;
        integral=0.0;
        sum=0.0;
      
        fprintf(outp1, "%s\n", label);                     //headers for output files
        fprintf(outp2, "----- %s ----", label);
        fprintf(outp3, "%s\n", label);                     //print atom names to file for labeling in graphs
        printf("%d - %s\n", counter, label);               //Print index of atoms in file, used for selecting
                                                         //what is to be graphed

        for (n=0; n<17; n++){                              //Loops over the twelve R's in data file
	        fscanf(inp, "%f", &rn[n]);                           //and places them in array for equation
	        fprintf(outp4, "%2.5g ", rn[n]);  
        }
      
        fscanf(inp, "%f", &RR);
        fprintf(outp4, " %0.2f\n", RR);

	    if (pick==2) az=(float)A/Z;
	    else az=1;

        for (n=0; n<1000; n++){                             //Loop for each radius, 0.01 intervals
	        term=0.0;
	        i=1;                                            //Reset both A and p to zero for each new radius 
	        for (j=0; j<17; j++){                           //loop for each R and Q, making sure that p=sum of
	            x=i*PI*r/RR;                                  //p's for each R and Q used.
	            term=term+rn[j]*sin(x)/x;
	            i++;
	        }  
	        integral=r*r*r*r*term*0.01;
            sum+=integral;
	        marker++;                                  
	        results[n][0]=r;                                //placing r and density into results array
	        results[n][counter]=term*az;                       //setting charge density per unit charge
	        fprintf(outp1, "%f %g\n", results[n][0], results[n][counter]);
	        r=r+0.01;     
	        if (term<0.0) n=1000;                             //adds 0.01 to r and begin another loop
                                                              //if term<0, breaks loop
        }
      
        loopcounter[counter-1]=marker;
        counter++;                                        //places next density results in next columns of array
        fprintf(outp1, "\n\n");
        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(outp6, "%d %lf %lf\n", A, densityRe, r2);
      
    }

    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[]
    int option;

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

    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);
        }
        fprintf(outp2, "\n");
    }

    fclose(inp);
    fclose(outp1);                                        //Close input and output files
    fclose(outp2);
    fclose(outp3);
    fclose(outp4);
    fclose(outp6);
  
    for(n=0; n<a; n++){
        FILE *inp2=fopen("titles.tmp", "r");                //open file each loop so in order to
        if (choice[n]==counter){                            //read from the beginning each time
            FILE *GNUoutp=fopen("FB_GNUoutp_All.gnu", "w");
            
            if (pick==2) fprintf(GNUoutp, "set term pdf\n set output \"FB_AZ.pdf\"\n");
            else fprintf(GNUoutp, "set term pdf\nset output \"FB.pdf\"\n");
            
            fprintf(GNUoutp, "set key outside\n");
            fprintf(GNUoutp, "set title \"Fourier Bessel Charge Distribution\"\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.1]\nset xrange [0:10]\n");
            else fprintf(GNUoutp, "set yrange [0:0.24]\nset xrange [0:10]\n");
            
            fprintf(GNUoutp, "plot \"FB_chargedistribution.dat\" using 1:2 w l lw 2 ");
            fprintf(GNUoutp, "title \"\"\n");                 //GNUplot formatting
            fprintf(GNUoutp, "unset output\n");
            fclose(GNUoutp);
            system("gnuplot < FB_GNUoutp_All.gnu");
        }
        else {
            FILE *GNUoutp=fopen("FB_GNUoutp.gnu", "w");
            FILE *GNUcharge=fopen("charge.tmp", "w");
            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=1; 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\n set output \"%s_FB_AZ.pdf\"\n", label);
            else fprintf(GNUoutp, "set term pdf\nset output \"%s_FB.pdf\"\n", label);
            
            fprintf(GNUoutp, "set key outside\n");
            fprintf(GNUoutp, "set title \"Fourier Bessel Charge Distribution 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.1]\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 w l lw 4 ");
            fprintf(GNUoutp, "title \"\"\n");                 //GNUplot formatting
            fprintf(GNUoutp, "unset output\n");
            fclose(GNUoutp);
            fclose(GNUcharge);
            system("gnuplot < FB_GNUoutp.gnu");
        }
        fclose(inp2);
        system("rm charge.tmp");							//removes temporary files from director
    }  

    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
    }

	int graphtype;

    if(a==0) graphtype = 10;                                      //exits if no atoms selected
    else {
        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("FB_GNUoutp_Group.gnu", "w");
  
        if (pick==2) fprintf(GNUoutp, "set term pdf\n set output \"FB_Group_AZ.pdf\"\n");
        else fprintf(GNUoutp, "set term pdf\nset output \"FB_Group.pdf\"\n");
    
        fprintf(GNUoutp, "set key outside\n");
        fprintf(GNUoutp, "set title \"Fourier Bessel Charge Distribution\"\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.1]\nset xrange [0:10]\n");
        else fprintf(GNUoutp, "set yrange [0:0.24]\nset xrange [0:10]\n");
    
        fprintf(GNUoutp, "set style fill solid 1\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, "\"FB_chargedistribution_columns.dat\" ");
            fprintf(GNUoutp, "using 1:%d w l 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, "unset output\n");
        fclose(GNUoutp);
        system ("gnuplot < FB_GNUoutp_Group.gnu");
    }
  
    if (graphtype==2 || graphtype==3) {

        for(n=0; n<a; n++){
	        printf("%d\n", choice[n]-1);
	        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.00001) i=1000;
	            else {
	                fprintf(outp5, "%d %f %f\n", n+1, results[i][0], results [i][choice[n]-1]);
	            }
	        }
	        fprintf(outp5, "\n\n");
        }
        fclose(outp5);

	    FILE *GNUoutp=fopen("FB_GNUoutp_Fence.gnu", "w");

	    fprintf(GNUoutp, "set term pdf\n");
	
	    if (pick==2) fprintf(GNUoutp, "set output \"FB_Fence_AZ.pdf\"\n");
	    else fprintf(GNUoutp, "set output \"FB_fence.pdf\"\n");
	
	    fprintf(GNUoutp, "set key outside\n");
	    fprintf(GNUoutp, "set title \"Nuclear Charge Density (FB)\"\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 filledcurve 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, "\"FB_chargedistribution_fence.dat\" u ($1==%d?($2):1/0):($1):3 ti '%s' lw 1 lc %d", n+1, label, n+1);
	        if (n==0) fprintf(GNUoutp, "\n");
	        else fprintf(GNUoutp, ", ");
	
	    }
	    fclose(GNUoutp);
	    system("gnuplot < FB_GNUoutp_Fence.gnu");
	}
	
    system("rm titles.tmp");                               //remove temporary file
    FILE *GNUoutp=fopen("FB_NuclearDensity.gnu", "w");
    fprintf(GNUoutp, "set term pdf\nset output \"FB_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 \"FB_nucleardensity.dat\" u 1:2 ti \"\"\n");
    fclose(GNUoutp);
    system ("gnuplot < FB_NuclearDensity.gnu");
    return 0;
}