#include <stdio.h>
#include <math.h>
#include <stddef.h>
#include <stdlib.h>
#include <string.h>
#include "nrutil.h"
#include "nrutil.c"
#include "nrfitting.c"
#include "UVStruct.c"
#include "UVCent.c"
#include "UVFile.c"

void main(void)

{
	int i, j, k, m, i_pk;
	int n_col, n_file, n_row, ok;
	int n_line_col, n_line_row;
	int n1, n2, i_min, i_max;
	float **data, **bk_data;
	float line_wl, bw, bw_mult;
	float dum, min, max, min_unc, max_unc;
	float peak, fwhm, fwhm_unc, fraction, centroid, centroid_unc;
	char st[64];
	char sch_file[20], in_file[20], out_file[20], line_file[20];
	
	struct scan_header L;
	
	FILE *fp_sch, *fp_in, *fp_line, *fp_out;
	
	printf("Schedule File: ");				// input files
	scanf("%s", &sch_file);	
	printf("Output File: ");
	scanf("%s", &out_file);
	
	fp_out = fopen(out_file, "w");			// header for output
	fprintf(fp_out, "\t\t\tPeak\t\t\tCentroid\t\t\t\t\tBandwidth\n");
	fprintf(fp_out, "File\tSource\tLine\tSignal\tIntercept\tFit\t0.5\t0.25\t0.1\tBW\tFWHM\tFWHM\n\n");
	
	fp_sch = fopen(sch_file, "r");
	
	fscanf(fp_sch, "%s %s %d", &st, &st, &n_col);	// read schedule header
	fscanf(fp_sch, "%s %d", &st, &n_row);
	do{
		fscanf(fp_sch, "%s", &st);
	} while(strcmp(st, "Data:") != 0);
	
	while(fscanf(fp_sch, "%s", &in_file) != EOF){	// read formatted and line files
		fscanf(fp_sch, "%s", &line_file);
		
		fp_in = fopen(in_file, "r");				// open formatted file
		read_scan_header(fp_in, &L);				// read header
		data = matrix(1, L.n_row, 1, L.n_col);		// allocate data matrix	
		bk_data = matrix(1, L.n_row, 1, L.n_col);	// allocate background matrix
		read_data(fp_in, L.n_row, L.n_col, data);	// read data
		fclose(fp_in);								// close formatted file
		
		for(i = 1; i <= L.n_row; i++){				// copy data to bk_data
			bk_data[i][1] = data[i][1];
			bk_data[i][2] = data[i][2];
			bk_data[i][3] = data[i][3];
		}
		
		fp_line = fopen(line_file, "r");					// open line file
		fscanf(fp_line, "%s %s %d", &st, &st, &n_line_col);	// read header
		fscanf(fp_line, "%s %d", &st, &n_line_row);
		fscanf(fp_line, "%s %f %s", &st, &bw, &st);
		do{
			fscanf(fp_line, "%s", &st);
		}while(strcmp(st, "Data:") != 0);
		for(m = 1; m <= n_line_row; m++){
			fscanf(fp_line, "%f %f", &line_wl, &bw_mult);			// read wavelength & mult.
			if((line_wl > L.start_wl) && (line_wl < L.stop_wl)){	// determine if line in scan
				printf("%s\t%f\t%f\n", L.source, line_wl, bw_mult);
				bksub_calc(L.n_row, bw, line_wl, bw_mult, data, bk_data);	// subtract background
				i_pk = find_peak(L.n_row, bw, line_wl, bk_data);			// find peak
				if(i_pk != 0){
					fprintf(fp_out, "%s\t%s\t", in_file, L.source);				// store file and line
					fprintf(fp_out, "%.3f\t%.3f\t", line_wl, bk_data[i_pk][1]);	// store wl and peak
					
					peak_intercept_calc(i_pk, bw, bk_data, &peak);				// calc peak from
					/*printf("peak\t%f\n", peak);*/								// intercept
					fprintf(fp_out, "%.3f\t", peak);							// store peak
					
					peak_fit_calc(i_pk, bk_data, &peak);						// calc peak from
					/*printf("peak\t%f\n", peak);*/								// spline fit
					fprintf(fp_out, "%.3f\t", peak);							// store peak
					
					fraction = 0.5;												// centroid
					find_limits_fr(i_pk, L.n_row, fraction, bk_data, &i_min, &i_max);
					centroid_calc(i_min, i_max, bk_data, &centroid, &centroid_unc);
					/*printf("centroid\t%f\t+-%f\n", centroid, centroid_unc);*/
					fprintf(fp_out, "%.3f\t", centroid);
					
					fraction = 0.25;											// centroid
					find_limits_fr(i_pk, L.n_row, fraction, bk_data, &i_min, &i_max);
					centroid_calc(i_min, i_max, bk_data, &centroid, &centroid_unc);
					/*printf("centroid\t%f\t+-%f\n", centroid, centroid_unc);*/
					fprintf(fp_out, "%.3f\t", centroid);
					
					fraction = 0.1;												// centroid
					find_limits_fr(i_pk, L.n_row, fraction, bk_data, &i_min, &i_max);
					centroid_calc(i_min, i_max, bk_data, &centroid, &centroid_unc);
					/*printf("centroid\t%f\t+-%f\n", centroid, centroid_unc);*/
					fprintf(fp_out, "%.3f\t", centroid);
					
					find_limits_bw(i_pk, L.n_row, bw*bw_mult, bk_data, &i_min, &i_max);	// centroid
					centroid_calc(i_min, i_max, bk_data, &centroid, &centroid_unc);
					/*printf("centroid\t%f\t+-%f\n", centroid, centroid_unc);*/
					fprintf(fp_out, "%.3f\t", centroid);
					
					fwhm_calc(i_pk, bk_data, &min, &max, &min_unc, &max_unc);		// centroid
					centroid = (max + min)/2.0;
					centroid_unc = sqrt(SQR(min_unc) + SQR(max_unc))/2.0;
					/*printf("centroid\t%f\t+-%f\n", centroid, centroid_unc);*/
					fprintf(fp_out, "%.3f\t", centroid);
					
					fwhm = max - min;											// bandwidth
					fwhm_unc = sqrt(SQR(min_unc) + SQR(max_unc));
					/*printf("bandwidth\t%f\t+-%f\n", fwhm, fwhm_unc);*/
					fprintf(fp_out, "%.3f\t", fwhm);
					
					fprintf(fp_out, "\n");
				}	
			}
		}
		fclose(fp_line);
		
	}
	fclose(fp_sch);	
	fclose(fp_out);
		
}