#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;
	float min_wl, max_wl, inc_wl;
	float bw, cutoff;
	float *a_hour, *zenith, *azimuth, *wl;
	float *Conv, **Irrad, **M_irrad, **Filter;
	float *filt_wl, *filt, *filt2;
	float float_in, time, x_in;
	char st[32], inst[8];
	char time_file[16], in_file[16], out_file[16], base_file[16];
	char wl_file[16];
	
	struct scan_header 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);
	}
	
	/* Enter parameters for analysis */
	
	printf("\nConvolution: \t0 for None\n\t\t1 for Spline fit\n\t\t2 for Triangle\n\t\t");
	printf("3 for Rectangle\n\t\t4 for Gaussian\n");
	scanf("%d", &conv_fnc);								// input analysis choice 
	
	if((conv_fnc != 0) && (conv_fnc != 1) && (conv_fnc != 5)){
		printf("Bandwidth (nm): ");						// enter bandwidth 
		scanf("%f", &bw);
	}
	if(conv_fnc != 5){
		printf("Min. WL, Max. WL, Inc. WL: ");			// enter wavelengths 
		scanf("%f,%f,%f", &min_wl, &max_wl, &inc_wl);
		n_wl = (int)(ceil((max_wl - min_wl)/inc_wl) + 1.0);			// calculate no. of wavelengths
		wl = vector(1, n_wl);
		for(i = 1; i <= n_wl; i++){						// calculate wavelengths
			wl[i] = min_wl + (float)(i - 1)*inc_wl;
		}
	}
	if(conv_fnc == 4){
		printf("Cut-off WL: ");							//cut-off wavelength
		scanf("%f", &cutoff);
	}
	
	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);
		
		if(conv_fnc == 1) conv_spline(n_wl, wl, hd_irrad.n_row,
			M_irrad, Conv);
		if(conv_fnc == 2) conv_tri(bw, n_wl, wl, 
			hd_irrad.n_row, M_irrad, Conv);
		if(conv_fnc == 3) conv_rect(bw, n_wl, wl, 
			hd_irrad.n_row, M_irrad, Conv);
		if(conv_fnc == 4) conv_gauss(bw, n_wl, wl, cutoff, 
			hd_irrad.n_row, M_irrad, Conv);
		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");
	store_header(fp_out, time_file, conv_fnc, bw, min_wl, max_wl, inc_wl,
		n_time, a_day, a_hour, zenith, azimuth);
	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;
}