/*Multiple Model Charge Distribution
HO, MHO, 2pF, 3pF, 2pG, 3pG
 */

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

int main(){
    char d;                                           //counts number of atoms
    int atomcount;     								//in the data file
    FILE *count=fopen ("tabletr.txt","r");			//used to seed arrays
    atomcount=0;
    if (count==NULL) perror ("Error opening file");		
    else{
        do {
            d = fgetc (count);
            if (d == '\n') atomcount++;
        } while (d != EOF);
        fclose (count);
    }

    char atom[10], label[15], ref[10];
    int Z, A, model, status, pick, i=1, counter, work, v;
    float c, z, w, r2, rtest=0.0, rmax, az;
    double ind, sum, result, t=0, tstep=.001, p_0;
    float f, x, Re, Rec, densityRe, densityRec, r2c;		//r2c is calculate r2 value
    float  integral, integrate;
    struct inputdata{
        float p_0;
        float c;
        float z;
        float w;
        int model;
        char label[15];
        float r2;
        char ref[6];
    } results[atomcount];

    FILE *inp=fopen("tabletr.txt","r");					//Data file
    FILE *outp=fopen("titles.tmp", "w");					//written file with atom names for labeling
    FILE *outp2=fopen("NuclearCharge_nucleardensity.dat", "w");	//data file for n v A
    int r2issue = 0;

    if (inp==NULL){
        printf("(tabletr.txt)\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 yes (type 1) or no (type 2)\n");
    }  
                                                    	//Skips first line of file
    while (d != '\n') d=fgetc(inp);                      	//that only contains column titles
    counter=1;
  
    fprintf(outp2, "A n Calculated n r2 Calculated r2\n");

    while (1){
	    sum=0.0;
	    x=0.0;
        integral=0.0;
        status=fscanf(inp,"%s %d %d %d %f %f %f %f %s\n",atom,&Z,&A,&model,&c,&z,&w,&r2,ref);
        if (status==EOF) break;								//scans in data and breaks at EOF 

        if (model == 1){
            sum=0;
            for (t=0; t<10+tstep; t+=tstep) {
	            ind=(1+z*pow(t/c,2))*exp(-pow(t/c,2))*t*t*tstep;
	            sum+=ind;
            }
            strcpy(label,atom);
            strcat(label,"-HO");
            printf("%d - %s - HO - c=%.3lf - z=%.3lf\n", counter, atom, c, z);
            fprintf(outp, "%s-HO\n", atom);
        }
        
        if (model == 2){
            sum=0;
            for (t=0; t<10+tstep; t+=tstep) {
	            ind=(1+z*pow(t/c,2))*exp(-pow(t/c,2))*t*t*tstep;
	            sum+=ind;
            }
            strcpy(label,atom);
            strcat(label,"-MHO");
            printf("%d - %s - MHO - c=%.3lf - z=%.3lf\n", counter, atom, c, z);
            fprintf(outp, "%s-MHO\n", atom);
        }
        
        if (model == 3) {
            sum=0;
            for (t=0; t<10+tstep; t+=tstep) {
	            ind=(1.0/(1+exp((t-c)/z)))*t*t*tstep;
	            sum+=ind; 
            }
            strcpy(label,atom);
            strcat(label,"-2pF");
            printf("%d - %s - 2pF - c=%.3lf - z=%.3lf\n", counter, atom, c, z);
            fprintf(outp, "%s-2pF\n", atom);
        }
        
        if (model == 4) {
            sum=0;
            for (t=0; t<10+tstep; t+=tstep) {
	            ind=tstep*t*t*(1+w*t*t/pow(c,2))/(1+exp((t-c)/z));    
	            sum+=ind;
            }
            strcpy(label,atom);
            strcat(label,"-3pF");
            printf("%d - %s - 3pF - c=%.3lf - z=%.3lf - w=%.3lf\n", counter, atom, c, z, w);
            fprintf(outp, "%s-3pF\n", atom);
        }
        
        if (model == 5) {
            sum=0;
            for (t=0; t<10+tstep; t+=tstep) {
	            ind=tstep*t*t*(1+w*t*t/pow(c,2))/(1+exp((t*t-c*c)/pow(z,2)));    
	            sum+=ind;
            }
            strcpy(label,atom);
            strcat(label,"-3pG");
            printf("%d - %s - 3pG - c=%.3lf - z=%.3lf - w=%.3lf\n", counter, atom, c, z, w);
            fprintf(outp, "%s-3pG\n", atom);
        }
        
        if (model == 6) {
            sum=0;
            for (t=0; t<10+tstep; t+=tstep) {
	            ind=tstep*t*t*1.0/(1+exp((t*t-c*c)/pow(z,2)));   
	            sum+=ind;
            }  
            strcpy(label,atom);
            strcat(label,"-2pG");
            printf("%d - %s - 2pG =c=%.3lf - z=%.3lf\n", counter, atom, c, z);      
            fprintf(outp, "%s-2pG\n", atom);
        }
        
        result=4.0*M_PI*sum;
        if (pick==2) az=(float)A/Z;
	    else az=1;	
        p_0=(double)Z/result*az;  
        results[counter].p_0=p_0;
        results[counter].c=c;
        results[counter].z=z;
        results[counter].w=w;
        results[counter].model=model;
        results[counter].r2=r2;
        strcpy(results[counter].label,label);
    
   
        for (v=0; v<1000; v++){
        	if (model == 1) f=p_0*(1+z*(x/c)*(x/c))*exp(-(x/c)*(x/c));
        	if (model == 2) f=p_0*(1+z*(x/c)*(x/c))*exp(-(x/c)*(x/c));	
        	if (model == 3) f=p_0/(1+exp((x-c)/z));
        	if (model == 4) f=p_0*(1+w*x*x/(c*c))/(1+exp((x-c)/z));
        	if (model == 5) f=p_0*(1+w*x*x/(c*c))/(1+exp((x*x-c*c)/(z*z)));
        	if (model == 6) f=p_0/(1+exp((x*x-c*c)/(z*z)));
        	
        	integrate=x*x*x*x*f*0.01;
        	integral+=integrate;
        	x=x+0.01;
        }
    
        r2c = pow(4*PI/Z*integral, 0.5);
        Re = pow(5*r2*r2/3, 0.5);
        Rec = pow(5*r2c*r2c/3, 0.5);
	    densityRe = (3.*A)/(4.*PI*Re*Re*Re);
	    densityRec= (3.*A)/(4*PI*Rec*Rec*Rec);
     	fprintf(outp2, "%d %lf %lf %lf %lf\n", A, densityRe, densityRec, r2, r2c);
        
    }

    printf("%d - done\n", counter);
    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);

    fclose(outp);
    fclose(outp2);
    fclose(inp);

    FILE *GNUoutp=fopen("GEN_NuclearDensity.gnu", "w");
    fprintf(GNUoutp, "set term pdf\nset output \"NuclearCharge_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 \"NuclearCharge_nucleardensity.dat\" u 1:3 ti \"\", ");
    fclose(GNUoutp);
    system ("gnuplot < GEN_NuclearDensity.gnu");
    int choice[counter];
    int a=0, option, n;

    while (option != counter){                           //Loop to scan in indices from terminal in order
        scanf("%d", &option);                              //to choose which atoms to graph
        if (option > counter || option < 1)
            printf("Not an available option, pick again\n");
        else if (option==counter);
        else {
            choice[a]=option; 
            a++;
        }
    }
    for (n=0; n<a; n++){
        FILE *GNUoutp=fopen("GNUoutp.gnu", "w");
        FILE *inp2=fopen("titles.tmp", "r");
    
        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
        }
        
        p_0=results[choice[n]].p_0;
        c=results[choice[n]].c;
        z=results[choice[n]].z;
        w=results[choice[n]].w;
        model=results[choice[n]].model;
        //strcpy(label,results[choice[n]].label);

        if (pick==2) fprintf(GNUoutp, "set term pdf\nset output \"NuclearCharge_%s_AZ.pdf\"\n", label);
        else fprintf(GNUoutp, "set term pdf\nset output \"NuclearCharge_%s.pdf\"\n", label);
        fprintf(GNUoutp, "set title \"Nuclear Charge Density of %s\"\n", label);
        fprintf(GNUoutp, "set xrange [0:10]\n");
        
        if (pick==2) fprintf(GNUoutp, "set yrange [0:0.24]\n");
        else fprintf(GNUoutp, "set yrange [0:0.12]\n");
        
        fprintf(GNUoutp, "set xlabel \"r (fm)\"\n");
        fprintf(GNUoutp, "set ylabel \"Density (nucleons/fm^3)\"\n");
        fprintf(GNUoutp, "p_0=%lf\nc=%lf\nz=%lf\nw=%lf\n", p_0,c,z,w);  
        
        if (model == 1) fprintf(GNUoutp, "f(x)=p_0*(1+z*(x/c)*(x/c))*exp(-(x/c)*(x/c))\n");
        if (model == 2) fprintf(GNUoutp,"f(x)=p_0*(1+z*(x/c)*(x/c))*exp(-(x/c)*(x/c))\n");
        if (model == 3) fprintf(GNUoutp,"f(x)=p_0/(1+exp((x-c)/z))\n");
        if (model == 4) fprintf(GNUoutp,"f(x)=p_0*(1+w*x*x/(c*c))/(1+exp((x-c)/z))\n");
        if (model == 5) fprintf(GNUoutp,"f(x)=p_0*(1+w*x*x/(c*c))/(1+exp((x*x-c*c)/(z*z)))\n");
        if (model == 6) fprintf(GNUoutp,"f(x)=p_0/(1+exp((x*x-c*c)/(z*z)))\n");
        
        fprintf(GNUoutp, "plot f(x) title \'%s\' with lines lw 2\n", label);
        fprintf(GNUoutp, "unset output\n");
        fclose(inp2);
        fclose(GNUoutp);
        system ("gnuplot < GNUoutp.gnu");
    }

    FILE *inp2=fopen("titles.tmp", "r");
    printf("Which elements do you want to graph together?\n");
    for (n=1; n<counter; n++){                             //Loop to print indices and atom names
        fscanf(inp2, "%s", label);                           //to screen for selecting
        c=results[n].c;
        z=results[n].z;
        w=results[n].w;
        printf("%d - %s - c=%.3lf - z=%.3lf - w=%.3lf\n", n, label, c, z, w);
    }
    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;                                  //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) return 0;                                      //exits if no atoms selected
  
    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("GNUoutp_Group.gnu", "w");
  
        if (pick==2) fprintf(GNUoutp, "set term pdf\nset output \"NuclearCharge_Group_AZ.pdf\"\n");
        else fprintf(GNUoutp, "set term pdf\nset output \"NuclearCharge_Group.pdf\"\n");
        
        fprintf(GNUoutp, "set key outside\n");
        fprintf(GNUoutp, "set title \"Nuclear Charge Density\"\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:9]\n");
        else fprintf(GNUoutp, "set yrange [0:0.12]\nset xrange [0:9]\n");
        
        for(n=0; n<a; n++){
	        p_0=results[choice[n]].p_0;
	        c=results[choice[n]].c;
	        z=results[choice[n]].z;
	        w=results[choice[n]].w;
	        model=results[choice[n]].model;
	        fprintf(GNUoutp, "p_0%d=%lf\nc%d=%lf\nz%d=%lf\nw%d=%lf\n", n, p_0,n, c, n, z,n, w);
	        
	        if (model == 1) fprintf(GNUoutp,"f%d(x)=p_0%d*(1+z%d*(x/c%d)*(x/c%d))*exp(-(x/c%d)*(x/c%d))\n", n,n,n,n,n,n,n);
	        if (model == 2) fprintf(GNUoutp,"f%d(x)=p_0%d*(1+z%d*(x/c%d)*(x/c%d))*exp(-(x/c%d)*(x/c%d))\n", n,n,n,n,n,n,n);
	        if (model == 3) fprintf(GNUoutp,"f%d(x)=p_0%d/(1+exp((x-c%d)/z%d))\n", n,n,n,n);
	        if (model == 4) fprintf(GNUoutp,"f%d(x)=p_0%d*(1+w%d*x*x/(c%d*c%d))/(1+exp((x-c%d)/z%d))\n", n,n,n,n,n,n,n);
	        if (model == 5) fprintf(GNUoutp,"f%d(x)=p_0%d*(1+w%d*x*x/(c%d*c%d))/(1+exp((x*x-c%d*c%d)/(z%d*z%d)))\n", n,n,n,n,n,n,n,n,n);
	        if (model == 6) fprintf(GNUoutp,"f%d(x)=p_0%d/(1+exp((x*x-c%d*c%d)/(z%d*z%d)))\n", n,n,n,n,n,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=0; 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, "f%d(x) title \'%s\' with lines lw 2 lc %d", n, label, n+1);
            
            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 < GNUoutp_Group.gnu");
    }  
    if(graphtype==2 || graphtype==3){
        for(n=0; n<a; n++){
            if (rtest<results[choice[n]].r2){
	            rmax=2.1*results[choice[n]].r2;
	            rtest=results[choice[n]].r2;
	        }
        }
        
      	FILE *GNUoutp=fopen("GNUoutp_Fence.gnu", "w");
      	
      	if (pick==2) fprintf(GNUoutp, "set term pdf\nset output \"NuclearCharge_Fence_AZ.pdf\"\n");
        else fprintf(GNUoutp, "set term pdf\nset output \"NuclearCharge_Fence.pdf\"\n");
        
        fprintf(GNUoutp, "set title \"Nuclear Charge Density\"\n");
        fprintf(GNUoutp, "set key outside\n");
        fprintf(GNUoutp, "set xlabel \"Radius (fm)\"\n");
        fprintf(GNUoutp, "set zlabel \"Charge Density\\n\(e*fm^-3)\" offset -5\n");
        fprintf(GNUoutp, "set yrange [%d:0] #range for counter\n", a);
        fprintf(GNUoutp, "set xrange [0:10] #range for radius\n");
        fprintf(GNUoutp, "set view 72,355\n");
        fprintf(GNUoutp, "set style function impulses\n");
        fprintf(GNUoutp, "set samp 200\nset iso 500\n");
        fprintf(GNUoutp, "yf(y, a) = (abs(y-a)<.02) ? 1. : NaN\n");
        fprintf(GNUoutp, "set xyplane 0\n");
        for(n=0; n<a; n++){
            p_0=results[choice[n]].p_0;
            c=results[choice[n]].c;
            z=results[choice[n]].z;
            w=results[choice[n]].w;
            model=results[choice[n]].model;
            fprintf(GNUoutp, "p_0%d=%lf\nc%d=%lf\nz%d=%lf\nw%d=%lf\n", n, p_0,n, c, n, z,n, w);

            if (model == 1) fprintf(GNUoutp, "f%d(x)=p_0%d*(1+z%d*(x/c%d)*(x/c%d))*exp(-(x/c%d)*(x/c%d))\n", n,n,n,n,n,n,n);
            if (model == 2) fprintf(GNUoutp,"f%d(x)=p_0%d*(1+z%d*(x/c%d)*(x/c%d))*exp(-(x/c%d)*(x/c%d))\n", n,n,n,n,n,n,n);
            if (model == 3) fprintf(GNUoutp,"f%d(x)=p_0%d/(1+exp((x-c%d)/z%d))\n", n,n,n,n);
            if (model == 4) fprintf(GNUoutp,"f%d(x)=p_0%d*(1+w%d*x*x/(c%d*c%d))/(1+exp((x-c%d)/z%d))\n", n,n,n,n,n,n,n);
            if (model == 5) fprintf(GNUoutp,"f%d(x)=p_0%d*(1+w%d*x*x/(c%d*c%d))/(1+exp((x*x-c%d*c%d)/(z%d*z%d)))\n", n,n,n,n,n,n,n,n,n);
            if (model == 6) fprintf(GNUoutp,"f%d(x)=p_0%d/(1+exp((x*x-c%d*c%d)/(z%d*z%d)))\n", n,n,n,n,n,n);

        }
        
        fprintf(GNUoutp, "splot ");
        for(n=a-1; n>-1; n--){
            FILE *inp2=fopen("titles.tmp", "r");                //Loop to read titles.tmp in order 
            for (i=0; 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, "yf(y,%d)*f%d(x) title \'%s\' lc %d", a-n, n, label, n+1);
            if (n==0) 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 < GNUoutp_Fence.gnu");
    }
    system("rm titles.tmp");
    return 0;
}