/*  ew.c  by CMS  12.10.02                                           */
/*  v1.0                                                             */
/*  Compile with "gcc ew.c -o ew"                                    */
/*  Run with "ew" and enter EW of blue & red lines.                  */
/*  Program will calculate doublet ratio along curve of growth       */
/*  and output C, F(C), F(2C), F(2C)/F(C), b and N(Na I).            */

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

/* Define the data file.                                             */
/* Define the "interpolate" function to interpolate between two      */
/* points (upperc,upperdr) and (lowerc, lowerdr) to find at what     */
/* x value the line goes through the value y.                        */

 #define DATA "doublet.data"
 #define interpolate(x1,y1,x2,y2,y) (((y-y1)*(x2-x1)/(y2-y1))+x1)

/* Define variables and start program */

 int main(char *argv[])
 {
    double dr, drclose, upperdr, lowerdr;
    double upperc, lowerc, input, wred, wblue;
    double upperfc, upperf2c, lowerfc, lowerf2c;
    double cinterp, fcinterp, f2cinterp, z, b, N, T;
    double *fc, *f2c, *ccc, *f2cfc;
    int s, drcloseindex, i, imax;
    FILE *datafile;

/* Allocate memory */

    fc = (double *)malloc(111111*sizeof(double));
    f2c = (double *)malloc(111111*sizeof(double));
    f2cfc = (double *)malloc(111111*sizeof(double));
    ccc = (double *)malloc(111111*sizeof(double));
    
/* Ask user for doublet ratio to interpolate to */

    printf("Doublet Ratio Method Calculator \n"); 
    printf("------------------------------- \n");
    printf("Enter W(5890):  ");
    scanf("%lf",&wblue);
    printf("Enter W(5896):  ");
    scanf("%lf",&wred);
    printf("Enter redshift (z):  ");
    scanf("%lf",&z);
    input = (wblue/wred);
    /* printf("input = %lf\n",input); */
    dr = input*((5895.9243*5895.9243)/(5889.9512*5889.9512));
    /* printf("doublet ratio = %lf\n",dr); */
/* Open file to read and read in all data, from 1 to imax. */

    datafile = fopen(DATA, "r");  
    i = 0;
    while(fscanf(datafile, "%lf %lf %lf %lf", 
		 &ccc[i], &f2cfc[i], &fc[i], &f2c[i]) == 4)
      i++;    
    imax = i;
    fclose(datafile);

/* Look for closest values to dr */
    drcloseindex = 0;
    drclose = f2cfc[0];

    for(i = 0; i < imax; i++) {    
    	if(fabs(f2cfc[i] - dr) < fabs(drclose - dr)) {
		drclose = f2cfc[i];                            
	  	drcloseindex = i;
	}
    }

/* Set upper values to closest one above given doublet ratio. */

      upperc = ccc[drcloseindex];
      upperdr = f2cfc[drcloseindex];
      upperfc = fc[drcloseindex];
      upperf2c = f2c[drcloseindex];      

/* Set lower limits as next lower values. */

      if(abs(f2cfc[drcloseindex-1]-dr) < abs(f2cfc[drcloseindex+1]-dr)) {
	lowerc = ccc[drcloseindex-1];
	lowerdr = f2cfc[drcloseindex-1];
	lowerfc = fc[drcloseindex-1];
	lowerf2c = f2c[drcloseindex-1];
      }
      else {
        lowerc = ccc[drcloseindex+1];
        lowerdr = f2cfc[drcloseindex+1];
	lowerfc = fc[drcloseindex+1];
	lowerf2c = f2c[drcloseindex+1];
      }

/* Interpolate closest values of C, F(C), & F(2C) using above defn. */
      cinterp = interpolate(upperc, upperdr, lowerc, lowerdr, dr);
      fcinterp = interpolate(upperfc, upperdr, lowerfc, lowerdr, dr);
      f2cinterp = interpolate(upperf2c, upperdr, lowerf2c, lowerdr, dr);
 
/* Print out resulting interpolated points */
      printf("Interpolation gives:\n");
      printf("C = %lf\t\t\tF(2C)/F(C) = %lf\n", cinterp, dr);
      printf("F(C) = %lf\t\t\tF(2C) = %lf\n", fcinterp, f2cinterp);

/* Calculate the b-value (cm/s) for the given equivalent widths. */

      if (input > 1 && input < 1.98 ) {
	b = (wred/5895.9243)*(1+z)*(3E10)/(2*fcinterp);
	T = (b/1.29E4)*(b/1.29E4)*23;
	N = 66.8*cinterp*b/(0.6311*5889.9512E-8);
	printf("\nb = %g cm/s \t\tT = %g K\n", b, T);
	printf("N(Na I) = %g cm^-2\n", N);
      } else {
	printf("Sorry, no can do!\n");
      }
/* Release memory. */
      free(fc);
      free(f2c);
      free(f2cfc);
      free(ccc);
      
/* Exit program. */
    return 0;
 }