#include <stdio.h>
#include <math.h>
#include <stddef.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#include "nrutil.h"
#include "nrutil.c"
#include "nrfitting.c"
#include "UVStruct.c"
#include "UVFIle.c"
#include "UVCommon.c"

main()
{
	int i, j, n;
	int n_row, n_col, n_file, n_stray;
	int ok, year, day, hour, tenth;
	int **label;
	float stray_wl, stray_sig, stray_unc, stray_sum, stray_sum2;
	float **sig, **resp, **irrad;
	float *resp_wl, *resp_val, *resp_val2;
	float *resp_Runc, *resp_Runc2, *resp_Sunc, *resp_Sunc2;
	float resp_fit, resp_Runc_fit, resp_Sunc_fit;
	float time, signal, aux;
	double integ, fract;
	char sch_file[16], resp_file[16], data_file[16], out_file[16];
	char st[80], file[16];
	char inst[8], sig_unit[16];
	
	struct scan_header hd_sig, hd_resp;
	struct label_unit unit;
	
	FILE *fp_sch, *fp_in, *fp_resp, *fp_out, *fp_data;
	
	printf("Schedule File: ");						//enter schedule file
	scanf("%s", &sch_file);
	printf("Stray-light Subtraction Limit (0 if none): ");
	scanf("%f", &stray_wl);
	
	hd_sig.start_wl = -9999;						//initialize optional headers
	hd_sig.stop_wl = -9999;
	hd_sig.inc_wl = -9999;
	hd_sig.dark_sig = -9999;
	hd_sig.filter_temp = -9999;
	hd_sig.detect_temp = -9999;
	hd_sig.motor_temp = -9999;
	hd_sig.housing_temp = -9999;
	hd_sig.pmt_volt = -9999;
	hd_sig.int_time = -9999;
	hd_sig.dead_time = -9999;
	hd_sig.n_chan = -9999;
			
	n_file = 1;
	fp_sch = fopen(sch_file, "r");					//open schedule file
	do{							
		fscanf(fp_sch, "%s", &st);
	} while(strcmp(st, "Data:") != 0);
	while(fscanf(fp_sch, "%s", &data_file) != EOF){			//read data file
	
		printf("Data\t%s\t", data_file);	
		fp_in = fopen(data_file, "r");	
		read_scan_header(fp_in, &hd_sig);					//read data header
		if(strcmp(hd_sig.operat, "Discrete") == 0){
			sig = matrix(1, hd_sig.n_row, 1, hd_sig.n_col);
		}
		label = imatrix(1, hd_sig.n_col, 1, 64);			//allocate for labels
		read_labels(fp_in, hd_sig.n_col, label, &unit);		//read labels
		/*read_data(fp_in, hd_sig.n_row, hd_sig.n_col, sig);	//read signals
		fclose(fp_in);*/
		
		fscanf(fp_sch, "%s", &resp_file);					//read resp file
		printf("Resp\t%s\t", resp_file);
		fp_resp = fopen(resp_file, "r");
		read_scan_header(fp_resp, &hd_resp);					//read resp header
		resp = matrix(1, hd_resp.n_row, 1, 4);
		read_data(fp_resp, hd_resp.n_row, hd_resp.n_col, resp);	//read resp
		fclose(fp_resp);				
		
		sscanf(data_file, "%2d%3dF%2d.%s", &year, &day, &n, &inst);	//create output file
		if(strcmp(hd_sig.operat, "Discrete") == 0){					//discrete files
			fract = modf(hd_sig.start_time, &integ);
			hour = (int)integ;
			fract *= 10.0;
			fract = modf(fract, &integ);
			if(fract >= 0.5) integ += 1.0;
			tenth = (int)integ;
			sprintf(out_file, "%2d%03d%02d%d.%s", year, day, hour, tenth, inst);
		}
		else{														//continuous files
			sprintf(out_file, "%2d%03d.%s", year, day, inst);
		}
		printf("Output\t%s\n", out_file);
		
		if(n_file == 1){								//create summary file
			sprintf(file, "%2dE.%s", year, inst);
			printf("%s\n", file);
			
			n_col = 19;
			if(hd_sig.start_wl > -9999) n_col++;
			if(hd_sig.stop_wl > -9999) n_col++;
			if(hd_sig.inc_wl > -9999) n_col++;
			if(hd_sig.dark_sig > -9999) n_col++;
			if(hd_sig.filter_temp > -9999) n_col++;
			if(hd_sig.detect_temp > -9999) n_col++;
			if(hd_sig.motor_temp > -9999) n_col++;
			if(hd_sig.housing_temp > -9999) n_col++;
			if(hd_sig.pmt_volt > -9999) n_col++;
			if(hd_sig.int_time > -9999) n_col++;
			if(hd_sig.dead_time > -9999) n_col++;
			if(hd_sig.n_chan > -9999) n_col++;
			if(stray_wl != 0) n_col += 2;
			
			fp_data = fopen(file, "w");	// compilation of header data
			fprintf(fp_data, "Header:\n");
			fprintf(fp_data, "Number_of_columns\t%d\n", n_col);
			fprintf(fp_data, "Labels:\n");
			fprintf(fp_data, "File\tInstrument\tSource\tLocation\tEnvironment\t");
			fprintf(fp_data, "Year\tDay\tDate\tOperation\t");
			fprintf(fp_data, "Measurement_purpose\t");
			fprintf(fp_data, "Number_of_scans\tNumber_of_columns\t");
			fprintf(fp_data, "Number_of_rows\t");
			fprintf(fp_data, "Start_time\tStop_time\t");
			if(hd_sig.start_wl > -9999) fprintf(fp_data, "Start_wavelength\t");
			if(hd_sig.stop_wl > -9999) fprintf(fp_data, "Stop_wavelength\t");
			if(hd_sig.inc_wl > -9999) fprintf(fp_data, "Increment_wavelength\t");
			if(hd_sig.n_chan > -9999) fprintf(fp_data, "Number_of_channels\t");
			if(hd_sig.dark_sig > -9999) fprintf(fp_data, "Dark_signal\t");
			if(hd_sig.filter_temp > -9999) fprintf(fp_data, "Filter_temperature\t");
			if(hd_sig.detect_temp > -9999) fprintf(fp_data, "Detector_temperature\t");
			if(hd_sig.motor_temp > -9999) fprintf(fp_data, "Motor_temperature\t");
			if(hd_sig.housing_temp > -9999) fprintf(fp_data, "Housing_temperature\t");
			if(hd_sig.pmt_volt > -9999) fprintf(fp_data, "PMT_voltage\t");
			if(hd_sig.int_time > -9999) fprintf(fp_data, "Integration_time\t");
			if(hd_sig.dead_time > -9999) fprintf(fp_data, "Dead_time\t");
			if(stray_wl != 0){
				fprintf(fp_data, "Stray_light_wavelength\tStray_light_signal\t");
			}
			fprintf(fp_data, "Comments\tFormatted_file\tResponsivity_file\t");
			fprintf(fp_data, "Analysis_file\n");
			fprintf(fp_data, "\t\t\t\t\t\t\t\t\t\t\t\t\t");
			fprintf(fp_data, "%s\t%s\t", hd_sig.start_time_unit, hd_sig.start_time_unit);
			if(hd_sig.start_wl > -9999) fprintf(fp_data, "[nm]\t[nm]\t[nm]\t");
			if(hd_sig.n_chan > -9999) fprintf(fp_data, "\t");
			if(hd_sig.dark_sig > -9999) fprintf(fp_data, "%s\t", hd_sig.dark_sig_unit);
			if(hd_sig.filter_temp > -9999) fprintf(fp_data, "[C]\t");
			if(hd_sig.detect_temp > -9999) fprintf(fp_data, "[C]\t");
			if(hd_sig.motor_temp > -9999) fprintf(fp_data, "[C]\t");
			if(hd_sig.housing_temp > -9999) fprintf(fp_data, "[C]\t");
			if(hd_sig.pmt_volt > -9999) fprintf(fp_data, "[V]\t");
			if(hd_sig.int_time > -9999) fprintf(fp_data, "[s]\t");
			if(hd_sig.dead_time > -9999) fprintf(fp_data, "[ns]\t");
			if(stray_wl != 0){
				fprintf(fp_data, "[nm]\t%s\t", unit.sig_unit);
			}
			fprintf(fp_data, "\n");
			fprintf(fp_data, "Data:\n");
		}
		
		if(strcmp(hd_sig.operat, "Discrete") == 0){
			resp_wl = vector(1, hd_resp.n_row);					//allocate for responsivity
			resp_val = vector(1, hd_resp.n_row);
			resp_Runc = vector(1, hd_resp.n_row);
			resp_Sunc = vector(1, hd_resp.n_row);
			resp_val2 = vector(1, hd_resp.n_row);
			resp_Runc2 = vector(1, hd_resp.n_row);
			resp_Sunc2 = vector(1, hd_resp.n_row);
			
			for(i = 1; i <= hd_resp.n_row; i++){
				resp_wl[i] = resp[i][1];
				resp_val[i] = resp[i][2];
				resp_Runc[i] = resp[i][3];
				resp_Sunc[i] = resp[i][4];
			}
																//spline fit responsivity
			spline(resp_wl, resp_val, hd_resp.n_row, 1e-35, 1e-35, resp_val2);
			spline(resp_wl, resp_Runc, hd_resp.n_row, 1e-35, 1e-35, resp_Runc2);
			spline(resp_wl, resp_Sunc, hd_resp.n_row, 1e-35, 1e-35, resp_Sunc2);
			
			read_data(fp_in, hd_sig.n_row, hd_sig.n_col, sig);	//read signals
			fclose(fp_in);
			
			irrad = matrix(1, hd_sig.n_row, 1, hd_sig.n_col + 1);
			
			stray_sig = 0;
			stray_unc = 0;
			if(stray_wl != 0){
				n_stray = 0;
				stray_sum = 0;
				stray_sum2 = 0;
				for(i = 1; i <= hd_sig.n_row; i++){
					if(sig[i][1] <= stray_wl){
						stray_sum += sig[i][2];
						stray_sum2 += sig[i][2]*sig[i][2];
						n_stray++;
					}
					else{
						break;
					}
				}
				stray_sig = stray_sum/(float)n_stray;
				stray_unc = sqrt((stray_sum2 - 2*stray_sig*stray_sum + n_stray*stray_sig*stray_sig)/((float)((n_stray-1)*n_stray)))/stray_sig;
			}
			for(i = 1; i <= hd_sig.n_row; i++){					//calculate irradiance
				irrad[i][1] = sig[i][1];
				splint(resp_wl, resp_val, resp_val2, hd_resp.n_row, sig[i][1], &resp_fit);
				splint(resp_wl, resp_Runc, resp_Runc2, hd_resp.n_row, sig[i][1], &resp_Runc_fit);
				splint(resp_wl, resp_Sunc, resp_Sunc2, hd_resp.n_row, sig[i][1], &resp_Sunc_fit);
				irrad[i][2] = (sig[i][2] - stray_sig)/resp_fit;
				sig[i][3] = sqrt(SQR(sig[i][3]*sig[i][2]) + SQR(stray_unc*stray_sig))/(sig[i][2] - stray_sig);
				irrad[i][3] = sqrt(SQR(sig[i][3]) + SQR(resp_Runc_fit));
				irrad[i][4] = resp_Sunc_fit;
			}
		}
		
		fp_out = fopen(out_file, "w");			// store header in files
		fprintf(fp_data, "%s\t", out_file);
		fprintf(fp_out, "Header:\n");
		fprintf(fp_out, "Instrument\t%s\n", hd_sig.inst);
		fprintf(fp_data, "%s\t", hd_sig.inst);
		fprintf(fp_out, "Source\t%s\n", hd_sig.source);
		fprintf(fp_data, "%s\t", hd_sig.source);
		fprintf(fp_out, "Location\t%s\n", hd_sig.locat);
		fprintf(fp_data, "%s\t", hd_sig.locat);
		fprintf(fp_out, "Environment\t%s\n", hd_sig.envir);
		fprintf(fp_data, "%s\t", hd_sig.envir);
		fprintf(fp_out, "Year\t%d\n", hd_sig.year);
		fprintf(fp_data, "%d\t", hd_sig.year);
		fprintf(fp_out, "Day\t%d\n", hd_sig.day);
		fprintf(fp_data, "%d\t", hd_sig.day);
		fprintf(fp_out, "Date\t%s\n", hd_sig.date);
		fprintf(fp_data, "%s\t", hd_sig.date);
		fprintf(fp_out, "Operation\t%s\n", hd_sig.operat);
		fprintf(fp_data, "%s\t", hd_sig.operat);
		fprintf(fp_out, "Measurement_purpose\t%s\n", hd_sig.meas_purp);
		fprintf(fp_data, "%s\t", hd_sig.meas_purp);
		fprintf(fp_out, "Number_of_scans\t1\n");
		fprintf(fp_data, "1\t");
		if(strcmp(hd_sig.operat, "Discrete") == 0){
			fprintf(fp_out, "Number_of_columns\t%d\n", hd_sig.n_col + 1);
			fprintf(fp_data, "%d\t", hd_sig.n_col + 1);
		}
		else{
			fprintf(fp_out, "Number_of_columns\t%d\n", hd_sig.n_col);
			fprintf(fp_data, "%d\t", hd_sig.n_col);
		}
		fprintf(fp_out, "Number_of_rows\t%d\n", hd_sig.n_row);
		fprintf(fp_data, "%d\t", hd_sig.n_row);
		fprintf(fp_out, "Start_time\t%.5f\t%s\n", hd_sig.start_time, hd_sig.start_time_unit);
		fprintf(fp_data, "%.5f\t", hd_sig.start_time);
		fprintf(fp_out, "Stop_time\t%.5f\t%s\n", hd_sig.stop_time, hd_sig.start_time_unit);
		fprintf(fp_data, "%.5f\t", hd_sig.stop_time);
		
		if(hd_sig.start_wl > -9999){
			fprintf(fp_out, "Start_wavelength\t%.2f\t[nm]\n", hd_sig.start_wl);
			fprintf(fp_data, "%.2f\t", hd_sig.start_wl);
			fprintf(fp_out, "Stop_wavelength\t%.2f\t[nm]\n", hd_sig.stop_wl);
			fprintf(fp_data, "%.2f\t", hd_sig.stop_wl);
			fprintf(fp_out, "Increment_wavelength\t%.2f\t[nm]\n", hd_sig.inc_wl);
			fprintf(fp_data, "%.2f\t", hd_sig.inc_wl);
		}
		if(hd_sig.n_chan > -9999){ 
			fprintf(fp_out, "Number_of_channels\t%d\n", hd_sig.n_chan);
			fprintf(fp_data, "%d\t", hd_sig.n_chan);
		}
		if(hd_sig.dark_sig > -9999){ 
			fprintf(fp_out, "Dark_signal\t%.5f\t%s\n", hd_sig.dark_sig, hd_sig.dark_sig_unit);
			fprintf(fp_data, "%.5f\t", hd_sig.dark_sig);
		}
		if(hd_sig.filter_temp > -9999){ 
			fprintf(fp_out, "Filter_temperature\t%.2f\t[C]\n", hd_sig.filter_temp);
			fprintf(fp_data, "%.2f\t", hd_sig.filter_temp);
		}
		if(hd_sig.detect_temp > -9999){ 
			fprintf(fp_out, "Detector_temperature\t%.2f\t[C]\n", hd_sig.detect_temp);
			fprintf(fp_data, "%.2f\t", hd_sig.detect_temp);
		}
		if(hd_sig.motor_temp > -9999){ 
			fprintf(fp_out, "Motor_temperature\t%.2f\t[C]\n", hd_sig.motor_temp);
			fprintf(fp_data, "%.2f\t", hd_sig.motor_temp);
		}
		if(hd_sig.housing_temp > -9999){ 
			fprintf(fp_out, "Housing_temperature\t%.2f\t[C]\n", hd_sig.housing_temp);
			fprintf(fp_data, "%.2f\t", hd_sig.housing_temp);
		}
		if(hd_sig.pmt_volt > -9999){ 
			fprintf(fp_out, "PMT_voltage\t%.2f\t[V]\n", hd_sig.pmt_volt);
			fprintf(fp_data, "%.2f\t", hd_sig.pmt_volt);
		}
		if(hd_sig.int_time > -9999){ 
			fprintf(fp_out, "Integration_time\t%.4f\t[s]\n", hd_sig.int_time);
			fprintf(fp_data, "%.4f\t", hd_sig.int_time);
		}
		if(hd_sig.dead_time > -9999){ 
			fprintf(fp_out, "Dead_time\t%.2f\t[ns]\n", hd_sig.dead_time);
			fprintf(fp_data, "%.2f\t", hd_sig.dead_time);
		}
		if(stray_wl != 0){
			fprintf(fp_out, "Stray_light_wavelength\t%.2f\t[nm]\n", stray_wl);
			fprintf(fp_data, "%.2f\t", stray_wl);	
			fprintf(fp_out, "Stray_light_signal\t%.5f\t%s\n", stray_sig, unit.sig_unit);
			fprintf(fp_data, "%.5f\t", stray_sig);
		}
		fprintf(fp_out, "Comments\t%s\n", hd_sig.comments);
		fprintf(fp_data, "%s\t", hd_sig.comments);
		fprintf(fp_out, "Formatted_file\t%s\n", data_file);
		fprintf(fp_data, "%s\t", data_file);
		fprintf(fp_out, "Responsivity_file\t%s\n", resp_file);
		fprintf(fp_data, "%s\t", resp_file);
		fprintf(fp_out, "Analysis_program\tCalc_Irrad.c\n");
		fprintf(fp_data, "Calc_Irrad.c\n");
		fprintf(fp_out, "Labels:\n");
		if(strcmp(hd_sig.operat, "Discrete") == 0){
			fprintf(fp_out, "Wavelength\tIrradiance\tdE(R)/E\tdE(S)/E\n");
			fprintf(fp_out, "[nm]\t[mW_m^-2_nm^-1]\n");
			fprintf(fp_out, "Data:\n");
			for(i = 1; i <= hd_sig.n_row; i++){				// store data
				fprintf(fp_out, "%.2f\t%.5f\t%.4f\t%.4f\n", irrad[i][1], irrad[i][2], irrad[i][3], irrad[i][4]);
			}
			free_matrix(sig, 1, hd_sig.n_row, 1, hd_sig.n_col);
			free_matrix(irrad, 1, hd_sig.n_row, 1, hd_sig.n_col + 1);
			free_vector(resp_wl, 1, hd_resp.n_row);					
			free_vector(resp_val,1, hd_resp.n_row);
			free_vector(resp_Runc, 1, hd_resp.n_row);
			free_vector(resp_Sunc, 1, hd_resp.n_row);
			free_vector(resp_val2, 1, hd_resp.n_row);
			free_vector(resp_Runc2, 1, hd_resp.n_row);
			free_vector(resp_Sunc2, 1, hd_resp.n_row);
		}
		else{
			fprintf(fp_out, "Time\t");
			for(i = 1; i <= hd_resp.n_row; i++){
				fprintf(fp_out, "Irradiance_%.2f[nm]\t", resp[i][1]);
			}
			for(i = (hd_sig.n_chan + 2); i <= hd_sig.n_col; i++){
				j = 1;
				while(label[i][j] != '\0'){
					fprintf(fp_out, "%c", label[i][j]);
					j++;
				}
				fprintf(fp_out, "\t");
			}
			fprintf(fp_out, "\n");
			fprintf(fp_out, "[UTC]\t");
			for(i = 1; i <= hd_resp.n_row; i++){
				fprintf(fp_out, "[mW_m^-2_nm^-1]\t");
			}
			for(i = (hd_sig.n_chan + 2); i <= hd_sig.n_col; i++){
				j = 1;
				while(label[i][j] != '\0'){
					st[j-1] = label[i][j];
					j++;
				}
				st[j-1] = '\0';
				if(strcmp(st, "Dark_signal") == 0) fprintf(fp_out, "%s\t", unit.dark_sig_unit);
				if(strcmp(st, "Dark_J") == 0) fprintf(fp_out, "%s\t", unit.darkj_sig_unit);
				if(strcmp(st, "Dark_T") == 0) fprintf(fp_out, "%s\t", unit.darkt_sig_unit);
				if(strcmp(st, "Filter_temperature") == 0) fprintf(fp_out, "%s\t", unit.filt_temp_unit);
				if(strcmp(st, "Integration_time") == 0) fprintf(fp_out, "%s\t", unit.int_time_unit);
			}
			fprintf(fp_out, "\n");
			fprintf(fp_out, "Data:\n");
			do{
				fscanf(fp_in, "%s", &st);
			} while(strcmp(st, "Data:") != 0);
			while(fscanf(fp_in, "%f", &time) != EOF){
				fprintf(fp_out, "%.5f\t", time);
			
				for(i = 1; i <= hd_sig.n_chan; i++){					//calculate irradiance
					fscanf(fp_in, "%f", &signal);
					fprintf(fp_out, "%.5f\t", signal/resp[i][2]);
				}
				for(i = (hd_sig.n_chan + 2); i <= hd_sig.n_col; i++){
					fscanf(fp_in, "%f", &aux);
					fprintf(fp_out, "%.5f\t", aux);
				}
				fprintf(fp_out, "\n");
			}
			fclose(fp_in);
		}
		
		free_matrix(resp, 1, hd_resp.n_row, 1, 4);
		free_imatrix(label, 1, hd_sig.n_col, 1, 64);
		n_file++;
		fclose(fp_out);
	}
	fclose(fp_data);
	fclose(fp_sch);	
}