#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"

#define PI 3.14159

main()
{
	int i, j, n;
	int n_diff_row, n_total_row, n_file, n_col;
	int ok, yy;
	float **diff_sig, **total_sig, **direct_sig;
	float *diff_wl, *diff, *diff2, diff_fit;
	float *diff_unc, *diff_unc2, diff_unc_fit;
	float **irrad, *irr_wl, *irr, *irr_Runc, *irr_Sunc;
	float *irr2, *irr_Runc2, *irr_Sunc2;
	float irr_fit, irr_Runc_fit, irr_Sunc_fit, irr_drv;
	float **resp;
	char sch_file[16], diff_file[80], total_file[80];
	char lamp_file[16], unc_file[16], data_file[16], out_file[16];
	char st[80], file[16];
	char inst[8], sig_unit[16];
	
	struct scan_header diff_hd, total_hd, resp_hd;
	struct lamp_header lamp_hd;
	struct unc_data unc;
	
	FILE *fp_sch, *fp_in, *fp_out, *fp_data;
	
	printf("Schedule File: ");						//enter schedule file
	scanf("%s", &sch_file);
	sscanf(sch_file, "%2d%s", &yy, &st);
	i = 0, j = 0;									//copy instrument
	while(sch_file[i] != '.') i++;
	while(sch_file[i] != '\0'){
		inst[j] = sch_file[i+1];
		i++;
		j++;
	}
	inst[j] = '\0';
	
	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", &diff_file) != EOF){		//read diffuse file
		printf("Diffuse\t%s\t", diff_file);
		if(strcmp(diff_file, "none") != 0){				
			n_diff_row = number_rows(diff_file);		//allocate for diffuse file
			diff_sig = matrix(1, n_diff_row, 1, 3);
			diff_wl = vector(1, n_diff_row);
			diff = vector(1, n_diff_row);
			diff2 = vector(1, n_diff_row);
			diff_unc = vector(1, n_diff_row);
			diff_unc2 = vector(1, n_diff_row);
														//calculate diffuse signals
			calc_signal(diff_file, n_diff_row, &diff_hd, diff_sig, sig_unit);
			for(i = 1; i <= n_diff_row; i++){
				diff_wl[i] = diff_sig[i][1];			//diffuse vectors
				diff[i] = diff_sig[i][2];
				diff_unc[i] = diff_sig[i][3];
			}
		}
		fscanf(fp_sch, "%s", total_file);			//read total file
		printf("Total\t%s\n", total_file);
		
		n_total_row = number_rows(total_file);		//allocate for total file
		total_sig = matrix(1, n_total_row, 1, 3);
		direct_sig = matrix(1, n_total_row, 1, 3);
														//calculate total signals
		calc_signal(total_file, n_total_row, &total_hd, total_sig, sig_unit);
		
		if(strcmp(diff_file, "none") != 0){			//calculate direct signals
			spline(diff_wl, diff, n_diff_row, 1e-35, 1e-35, diff2);
			spline(diff_wl, diff_unc, n_diff_row, 1e-35, 1e-35, diff_unc2);
			for(i = 1; i <= n_total_row; i++){
				direct_sig[i][1] = total_sig[i][1];
				splint(diff_wl, diff, diff2, n_diff_row, total_sig[i][1], &diff_fit);
				direct_sig[i][2] = total_sig[i][2] - diff_fit;
				splint(diff_wl, diff_unc, diff_unc2, n_diff_row, total_sig[i][1], &diff_unc_fit);
				direct_sig[i][3] = sqrt(SQR(total_sig[i][3]*total_sig[i][2]) + SQR(diff_unc_fit*diff_fit))/direct_sig[i][2];
			}
			free_matrix(diff_sig, 1, n_diff_row, 1, 3);
			free_vector(diff_wl, 1, n_diff_row);
			free_vector(diff, 1, n_diff_row);
			free_vector(diff2, 1, n_diff_row);
			free_vector(diff_unc, 1, n_diff_row);
			free_vector(diff_unc2, 1, n_diff_row);
		}
		else{
			for(i = 1; i <= n_total_row; i++){
				direct_sig[i][1] = total_sig[i][1];
				direct_sig[i][2] = total_sig[i][2];
				direct_sig[i][3] = total_sig[i][3];
			}
		}
		free_matrix(total_sig, 1, n_total_row, 1, 3);
		
		fscanf(fp_sch, "%s", &lamp_file);				//load lamp file
		printf("Lamp\t%s\t", lamp_file);
		fp_in = fopen(lamp_file, "r");
		
		read_lamp_header(fp_in, &lamp_hd);				//read lamp header
		
		irrad = matrix(1, lamp_hd.n_row, 1, 4);			//allocate for lamp file 
		irr_wl = vector(1, lamp_hd.n_row);
		irr = vector(1, lamp_hd.n_row);
		irr_Runc = vector(1, lamp_hd.n_row);
		irr_Sunc = vector(1, lamp_hd.n_row);
		irr2 = vector(1, lamp_hd.n_row);
		irr_Runc2 = vector(1, lamp_hd.n_row);
		irr_Sunc2 = vector(1, lamp_hd.n_row);
		
		read_data(fp_in, lamp_hd.n_row, lamp_hd.n_col, irrad);		//read lamp data
		
		fclose(fp_in);
		
		for(i = 1; i <= lamp_hd.n_row; i++){			//irradiance vectors
			irr_wl[i] = irrad[i][1];
			irr[i] = irrad[i][2];
			irr_Runc[i] = irrad[i][3];
			irr_Sunc[i] = irrad[i][4];
		}
		free_matrix(irrad, 1, lamp_hd.n_row, 1, 4);
														//spline fit irradiance
		spline(irr_wl, irr, lamp_hd.n_row, 1e-35, 1e-35, irr2);
		spline(irr_wl, irr_Runc, lamp_hd.n_row, 1e-35, 1e-35, irr_Runc2);
		spline(irr_wl, irr_Sunc, lamp_hd.n_row, 1e-35, 1e-35, irr_Sunc2);
		
		resp = matrix(1, n_total_row, 1, 4);			//allocate response 
		
		for(i = 1; i <= n_total_row; i++){				//calculate responsivity
			resp[i][1] = direct_sig[i][1];
			splint(irr_wl, irr, irr2, lamp_hd.n_row, direct_sig[i][1], &irr_fit);
			resp[i][2] = direct_sig[i][2]/irr_fit;
		}
		
		fscanf(fp_sch, "%s", &unc_file);			//load uncertainty file
		printf("Uncert\t%s\n", unc_file);
		read_unc(unc_file, &unc);					//read uncertainty data
		
		for(i = 1; i <= n_total_row; i++){
			resp[i][3] = 0;
			resp[i][4] = 0;
		}
		for(i = 1; i <= n_total_row; i++){
											//random uncertainties
			resp[i][3] += SQR((654.6/direct_sig[i][1])*0.0006*unc.curr_rand);
			resp[i][3] += SQR(direct_sig[i][3]);
			resp[i][3] = sqrt(resp[i][3]);
											//systematic uncertainties
			if(unc.radius > 0){	
			splint(irr_wl, irr_Runc, irr_Runc2, lamp_hd.n_row, direct_sig[i][1], &irr_Runc_fit);
			resp[i][4] += SQR(irr_Runc_fit);
			splint(irr_wl, irr_Sunc, irr_Sunc2, lamp_hd.n_row, direct_sig[i][1], &irr_Sunc_fit);
			resp[i][4] += SQR(irr_Sunc_fit);
			resp[i][4] += SQR(1.266474e-4 - 3.050758e-6*unc.radius - 3.947413e-4*SQR(unc.radius));
			resp[i][4] += SQR(atan(unc.radius/50.0)*90/PI*(1 - lamp_hd.avg_gonio));
			splint(irr_wl, irr, irr2, lamp_hd.n_row, direct_sig[i][1], &irr_fit);
			spldrv(irr_wl, irr, irr2, lamp_hd.n_row, direct_sig[i][1], &irr_drv);		
			resp[i][4] += SQR(irr_drv*(unc.wl_c0 + unc.wl_c1*direct_sig[i][1] + unc.wl_c2*SQR(direct_sig[i][1]))/irr_fit);	
			resp[i][4] += SQR(lamp_hd.max_gonio*unc.perp/sqrt(3.0));
			resp[i][4] += SQR(atan(unc.center/50.0)*180/PI*lamp_hd.max_gonio/sqrt(3.0));
			resp[i][4] += SQR(2.0/sqrt(3.0)*unc.dist/50.0);
			resp[i][4] = sqrt(resp[i][4]);
			}
		}										//free memory
		
		free_matrix(direct_sig, 1, n_total_row, 1, 3);
		free_vector(irr_wl, 1, lamp_hd.n_row);
		free_vector(irr, 1, lamp_hd.n_row);
		free_vector(irr2, 1, lamp_hd.n_row);
		free_vector(irr_Runc, 1, lamp_hd.n_row);
		free_vector(irr_Sunc, 1, lamp_hd.n_row);
		free_vector(irr_Runc2, 1, lamp_hd.n_row);
		free_vector(irr_Sunc2, 1, lamp_hd.n_row);
		
		if(n_file == 1){								//create summary file
			sprintf(data_file, "%2dR.%s", yy, inst);
			printf("%s\n", data_file);
			
			n_col = 23;
			if(total_hd.start_wl > -9999) n_col++;
			if(total_hd.stop_wl > -9999) n_col++;
			if(total_hd.inc_wl > -9999) n_col++;
			if(total_hd.dark_sig > -9999) n_col++;
			if(total_hd.filter_temp > -9999) n_col++;
			if(total_hd.detect_temp > -9999) n_col++;
			if(total_hd.motor_temp > -9999) n_col++;
			if(total_hd.housing_temp > -9999) n_col++;
			if(total_hd.pmt_volt > -9999) n_col++;
			if(total_hd.int_time > -9999) n_col++;
			if(total_hd.dead_time > -9999) n_col++;
			
			fp_data = fopen(data_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\tLamp_housing\tLamp_power_supply\t");
			fprintf(fp_data, "Number_of_scans\tNumber_of_columns\t");
			fprintf(fp_data, "Number_of_rows\tStart_time\tStop_time\t");
			if(total_hd.start_wl > -9999) fprintf(fp_data, "Start_wavelength\t");
			if(total_hd.stop_wl > -9999) fprintf(fp_data, "Stop_wavelength\t");
			if(total_hd.inc_wl > -9999) fprintf(fp_data, "Increment_wavelength\t");
			if(total_hd.dark_sig > -9999) fprintf(fp_data, "Dark_signal\t");
			if(total_hd.filter_temp > -9999) fprintf(fp_data, "Filter_temperature\t");
			if(total_hd.detect_temp > -9999) fprintf(fp_data, "Detector_temperature\t");
			if(total_hd.motor_temp > -9999) fprintf(fp_data, "Motor_temperature\t");
			if(total_hd.housing_temp > -9999) fprintf(fp_data, "Housing_temperature\t");
			if(total_hd.pmt_volt > -9999) fprintf(fp_data, "PMT_voltage\t");
			if(total_hd.int_time > -9999) fprintf(fp_data, "Integration_time\t");
			if(total_hd.dead_time > -9999) fprintf(fp_data, "Dead_time\t");
			fprintf(fp_data, "Comments\tDiffuse_file\tTotal_file\t");
			fprintf(fp_data, "Irradiance_file\tUncertainty_file\tAnalysis_file\n");
			fprintf(fp_data, "\t\t\t\t\t\t\t\t\t\t\t\t\t\t\t");
			fprintf(fp_data, "%s\t%s\t", total_hd.start_time_unit, total_hd.start_time_unit);
			if(total_hd.start_wl > -9999) fprintf(fp_data, "[nm]\t[nm]\t[nm]\t");
			if(total_hd.dark_sig > -9999) fprintf(fp_data, "%s\t", total_hd.dark_sig_unit);
			if(total_hd.filter_temp > -9999) fprintf(fp_data, "[C]\t");
			if(total_hd.detect_temp > -9999) fprintf(fp_data, "[C]\t");
			if(total_hd.motor_temp > -9999) fprintf(fp_data, "[C]\t");
			if(total_hd.housing_temp > -9999) fprintf(fp_data, "[C]\t");
			if(total_hd.pmt_volt > -9999) fprintf(fp_data, "[V]\t");
			if(total_hd.int_time > -9999) fprintf(fp_data, "[s]\t");
			if(total_hd.dead_time > -9999) fprintf(fp_data, "[ns]\t");
			fprintf(fp_data, "\n");
			fprintf(fp_data, "Data:\n");
		}
											// create formatted file
		i = 1;
		do{
			sprintf(out_file, "%2d%03dR%02d.%s", yy, total_hd.day, i, inst);
			ok = file_exist(out_file);
			i++;
		} while(ok == 1);
		
		fp_out = fopen(out_file, "w");			// store header in files
		printf("%s\n", out_file);
		fprintf(fp_data, "%s\t", out_file);
		fprintf(fp_out, "Header:\n");
		fprintf(fp_out, "Instrument\t%s\n", total_hd.inst);
		fprintf(fp_data, "%s\t", total_hd.inst);
		fprintf(fp_out, "Source\t%s\n", total_hd.source);
		fprintf(fp_data, "%s\t", total_hd.source);
		fprintf(fp_out, "Location\t%s\n", total_hd.locat);
		fprintf(fp_data, "%s\t", total_hd.locat);
		fprintf(fp_out, "Environment\t%s\n", total_hd.envir);
		fprintf(fp_data, "%s\t", total_hd.envir);
		fprintf(fp_out, "Year\t%d\n", total_hd.year);
		fprintf(fp_data, "%d\t", total_hd.year);
		if(strcmp(diff_file, "none") != 0){
			fprintf(fp_out, "Day\t%d\n", diff_hd.day);
			fprintf(fp_data, "%d\t", diff_hd.day);
		}
		else{
			fprintf(fp_out, "Day\t%d\n", total_hd.day);
			fprintf(fp_data, "%d\t", total_hd.day);
		}
		if(strcmp(diff_file, "none") != 0){
			fprintf(fp_out, "Date\t%s\n", diff_hd.date);
			fprintf(fp_data, "%s\t", diff_hd.date);
		}
		else{
			fprintf(fp_out, "Date\t%s\n", total_hd.date);
			fprintf(fp_data, "%s\t", total_hd.date);
		}
		fprintf(fp_out, "Operation\t%s\n", total_hd.operat);
		fprintf(fp_data, "%s\t", total_hd.operat);
		fprintf(fp_out, "Measurement_purpose\t%s\n", total_hd.meas_purp);
		fprintf(fp_data, "%s\t", total_hd.meas_purp);
		fprintf(fp_out, "Lamp_housing\t%s\n", total_hd.lamp_hs);
		fprintf(fp_data, "%s\t", total_hd.lamp_hs);
		fprintf(fp_out, "Lamp_power_supply\t%s\n", total_hd.lamp_ps);
		fprintf(fp_data, "%s\t", total_hd.lamp_ps);
		fprintf(fp_out, "Number_of_scans\t%d\n", number_files(total_file));
		fprintf(fp_data, "%d\t", number_files(total_file));
		fprintf(fp_out, "Number_of_columns\t4\n");
		fprintf(fp_data, "4\t");
		fprintf(fp_out, "Number_of_rows\t%d\n", total_hd.n_row);
		fprintf(fp_data, "%d\t", total_hd.n_row);
		if(strcmp(diff_file, "none") != 0){
			fprintf(fp_out, "Start_time\t%.5f\t%s\n", diff_hd.start_time, total_hd.start_time_unit);
			fprintf(fp_data, "%.5f\t", diff_hd.start_time);
		}
		else{
			fprintf(fp_out, "Start_time\t%.5f\t%s\n", total_hd.start_time, total_hd.start_time_unit);
			fprintf(fp_data, "%.5f\t", total_hd.start_time);
		}
		fprintf(fp_out, "Stop_time\t%.5f\t%s\n", total_hd.stop_time, total_hd.start_time_unit);
		fprintf(fp_data, "%.5f\t", total_hd.stop_time);
		
		if(total_hd.start_wl > -9999){
			fprintf(fp_out, "Start_wavelength\t%.2f\t[nm]\n", total_hd.start_wl);
			fprintf(fp_data, "%.2f\t", total_hd.start_wl);
			fprintf(fp_out, "Stop_wavelength\t%.2f\t[nm]\n", total_hd.stop_wl);
			fprintf(fp_data, "%.2f\t", total_hd.stop_wl);
			fprintf(fp_out, "Increment_wavelength\t%.2f\t[nm]\n", total_hd.inc_wl);
			fprintf(fp_data, "%.2f\t", total_hd.inc_wl);
		}
		if(total_hd.dark_sig > -9999){ 
			fprintf(fp_out, "Dark_signal\t%.5f\t%s\n", total_hd.dark_sig, total_hd.dark_sig_unit);
			fprintf(fp_data, "%.5f\t", total_hd.dark_sig);
		}
		if(total_hd.filter_temp > -9999){ 
			fprintf(fp_out, "Filter_temperature\t%.2f\t[C]\n", total_hd.filter_temp);
			fprintf(fp_data, "%.2f\t", total_hd.filter_temp);
		}
		if(total_hd.detect_temp > -9999){ 
			fprintf(fp_out, "Detector_temperature\t%.2f\t[C]\n", total_hd.detect_temp);
			fprintf(fp_data, "%.2f\t", total_hd.detect_temp);
		}
		if(total_hd.motor_temp > -9999){ 
			fprintf(fp_out, "Motor_temperature\t%.2f\t[C]\n", total_hd.motor_temp);
			fprintf(fp_data, "%.2f\t", total_hd.motor_temp);
		}
		if(total_hd.housing_temp > -9999){ 
			fprintf(fp_out, "Housing_temperature\t%.2f\t[C]\n", total_hd.housing_temp);
			fprintf(fp_data, "%.2f\t", total_hd.housing_temp);
		}
		if(total_hd.pmt_volt > -9999){ 
			fprintf(fp_out, "PMT_voltage\t%.2f\t[V]\n", total_hd.pmt_volt);
			fprintf(fp_data, "%.2f\t", total_hd.pmt_volt);
		}
		if(total_hd.int_time > -9999){ 
			fprintf(fp_out, "Integration_time\t%.4f\t[s]\n", total_hd.int_time);
			fprintf(fp_data, "%.4f\t", total_hd.int_time);
		}
		if(total_hd.dead_time > -9999){ 
			fprintf(fp_out, "Dead_time\t%.2f\t[ns]\n", total_hd.dead_time);
			fprintf(fp_data, "%.2f\t", total_hd.dead_time);
		}
			
		fprintf(fp_out, "Comments\t%s\n", total_hd.comments);
		fprintf(fp_data, "%s\t", total_hd.comments);
		fprintf(fp_out, "Diffuse_file\t%s\n", diff_file);
		fprintf(fp_data, "%s\t", diff_file);
		fprintf(fp_out, "Total_file\t%s\n", total_file);
		fprintf(fp_data, "%s\t", total_file);
		fprintf(fp_out, "Irradiance_file\t%s\n", lamp_file);
		fprintf(fp_data, "%s\t", lamp_file);
		fprintf(fp_out, "Uncertainty_file\t%s\n", unc_file);
		fprintf(fp_data, "%s\t", unc_file);
		fprintf(fp_out, "Analysis_program\tCalc_Resp.c\n");
		fprintf(fp_data, "Calc_Resp.c\n");
		fprintf(fp_out, "Labels:\n");
		fprintf(fp_out, "Wavelength\tResponsivity\tdR(R)/R\tdR(S)/R\n");
		fprintf(fp_out, "[nm]\t");
		j= 0;
		while(sig_unit[j] != ']'){
			fprintf(fp_out, "%c", sig_unit[j]);
			j++;
		}
		fprintf(fp_out, "_mW^-1_m^2_nm]\n");
		fprintf(fp_out, "Data:\n");
		for(i = 1; i <= total_hd.n_row; i++){				// store data
			fprintf(fp_out, "%.2f\t%.5f\t%.4f\t%.4f\n", resp[i][1], resp[i][2], resp[i][3], resp[i][4]);
		}
		
		fclose(fp_out);
		
		free_matrix(resp, 1, n_total_row, 1, 4);
		
		n_file++;
	}	
}