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

int main()
{

	int i, j, k;
	int n_time, conv_fnc, n_wl, n_row, n_col;
	int int_in;
	int *a_year, *a_day, *match;
	float *a_hour, *zenith, *azimuth, *wl;
	float *Conv, **Irrad, **M_irrad, **Filter;
	float *filt_wl, *filt, *filt2;
	float float_in, time, x_in;
	float sum, prod, wl_fit, spl_fit;
	char st[32], inst[8];
	char time_file[16], in_file[16], out_file[16], base_file[16], filt_file[16];
	char wl_file[16];
	
	struct scan_header hd_filt, hd_irrad;
	
	FILE *fp_in, *fp_file, *fp_out;
	
	/* Input time file and load instruments and times */
	
	printf("Times File: ");
	scanf("%s", &time_file);
								
	a_day = ivector(1, 100);							// allocate memory 
	a_year = ivector(1, 100);
	a_hour = vector(1, 100);
	
	fp_in = fopen(time_file, "r");						// open file 
	fscanf(fp_in, "%s%s%s", &st, &st, &st);
	n_time = 1;
	while(fscanf(fp_in, "%d", &int_in) != EOF){			// read in times 
		a_year[n_time] = int_in;
		fscanf(fp_in, "%d%f", &int_in, &float_in);
		a_day[n_time] = int_in;
		a_hour[n_time] = float_in;
		n_time++;
	}
	n_time--;											// no. of times 
	fclose(fp_in);
	
	printf("Instrument: ");
	scanf("%s", &inst);
	
	zenith = vector(1, n_time);							// allocate zenith and azimuth 
	azimuth = vector(1, n_time);
	
	for(i = 1; i <= n_time; i++){						// calculate angles 
		solar_angles(a_day[i], a_hour[i], 40.125, 105.237, zenith + i, azimuth + i);
	}
	
	/* Load filter functions  */
	
	printf("Filter file: ");
	scanf("%s", &filt_file);
	printf("Wavelength file: ");
	scanf("%s", &wl_file);
	
	fp_in = fopen(filt_file, "r");
	read_scan_header(fp_in, &hd_filt);
	Filter = matrix(1, hd_filt.n_row, 1, hd_filt.n_col); 
	read_data(fp_in, hd_filt.n_row, hd_filt.n_col, Filter);
	fclose(fp_in);
	
	fp_in = fopen(wl_file, "r");
	do{
		fscanf(fp_in, "%s", &st);
	} while(strcmp(st, "Number_of_rows") != 0);
	fscanf(fp_in, "%d", &n_wl);
	wl = vector(1, n_wl);
	do{
		fscanf(fp_in, "%s", &st);
	} while(strcmp(st, "Data:") != 0);
	for(i = 1; i <= n_wl; i++){
		fscanf(fp_in, "%s%f", &st, wl + i);
	}
	fclose(fp_in);
	printf("%d\n", n_wl);
	
	match = ivector(1, n_wl);
	
	printf("Output File: ");							// enter output file 
	scanf("%s", &base_file);
	
	Conv = vector(1, n_wl);								// allocate for convolution 
	Irrad = matrix(1, n_wl, 1, n_time);					// allocate for analysis files 
	
	/* Loop through files, loading and analyzing data */
	
	for(i = 1; i <= n_time; i++){						// loop through times 
		sprintf(in_file, "%2d%3d%3.0f.%s",a_year[i], a_day[i], a_hour[i]*10, inst);	// load AES data 
		printf("%s\t",in_file);
		fp_in = fopen(in_file, "r");
		read_scan_header(fp_in, &hd_irrad);
		M_irrad = matrix(1, hd_irrad.n_row, 1, hd_irrad.n_col);
		read_data(fp_in, hd_irrad.n_row, hd_irrad.n_col, M_irrad);
		fclose(fp_in);
		
		filt_wl = vector(1, hd_filt.n_row);
		filt = vector(1, hd_filt.n_row);
		filt2 = vector(1, hd_filt.n_row);
		
		for(j = 1; j <= n_wl; j++){
			for(k = 1; k <= hd_filt.n_row; k++){
				filt_wl[k] = Filter[k][1];
				filt[k] = Filter[k][j+1];
				if(filt[k] <= 0.00001) filt[j] = 0.00001;
			}
			spline(filt_wl, filt, hd_filt.n_row, 1e40, 1e40, filt2);
			
			sum = 0, prod = 0;
			for(k = 1; k <= hd_irrad.n_row; k++){
				wl_fit = M_irrad[k][1];
				if((wl_fit >= filt_wl[1]) && (wl_fit <= filt_wl[hd_filt.n_row])){
					splint(filt_wl, filt, filt2, hd_filt.n_row, wl_fit, &spl_fit);
					sum += spl_fit;
					prod += spl_fit*M_irrad[k][2];
				}
			}
			Conv[j] = prod/sum;
		}
		
		free_vector(filt_wl, 1, hd_filt.n_row);
		free_vector(filt, 1, hd_filt.n_row);
		free_vector(filt2, 1, hd_filt.n_row);
		
		free_matrix(M_irrad, 1, hd_irrad.n_row, 1, hd_irrad.n_col);
		/*for(j = 1; j <= n_wl; j++){
			printf("%d\t%.2f\t%f\n", j, wl[j], Conv[j]);
		}*/
		
		/* Compile results */
		
		for(j = 1; j <= n_wl; j++){
			Irrad[j][i] = Conv[j];
		}
		
	}
	
	/* Store Data */
	
	sprintf(out_file, "%s.irr", base_file);
	fp_out = fopen(out_file, "w");
	fprintf(fp_out, "Day:\t");
	for(i = 1; i <= n_time; i++) fprintf(fp_out, "%d\t", a_day[i]);
	fprintf(fp_out, "\n");
	fprintf(fp_out, "Time:\t");
	for(i = 1; i <= n_time; i++) fprintf(fp_out, "%.1f\t", a_hour[i]);
	fprintf(fp_out, "\n");
	fprintf(fp_out, "Zenith Angle [deg]:\t");
	for(i = 1; i <= n_time; i++) fprintf(fp_out, "%.4f\t", zenith[i]);
	fprintf(fp_out, "\n");
	fprintf(fp_out, "Azimuth Angle [deg]:\t");
	for(i = 1; i <= n_time; i++) fprintf(fp_out, "%.4f\t", azimuth[i]);
	fprintf(fp_out, "\n\n");
	fprintf(fp_out, "Wavelength\t");
	for(i = 1; i <= n_time; i++){
		fprintf(fp_out, "Irradiance\t");
	}
	fprintf(fp_out, "\n[nm]\t");
	for(i = 1; i <= n_time; i++){
		fprintf(fp_out, "[mW/m^2 nm]\t");
	}
	fprintf(fp_out, "\n\n");
	for(i = 1; i <= n_wl; i++){
		fprintf(fp_out, "%.2f\t", wl[i]);
		for(j = 1; j <= n_time; j++){
			fprintf(fp_out, "%.4f\t", Irrad[i][j]);
		}
		fprintf(fp_out, "\n");
	}
	fclose(fp_out);
	
	return 1;
}