#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, n_wl, n_row, n_col;
	int int_in;
	int *a_year, *a_day, *match;
	float *a_hour, *zenith, *azimuth, *wl;
	float *Conv, **Irrad;
	float *filt_wl, *filt, *filt2;
	float float_in, time;
	float wl_inc, time_inc;
	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_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("Wavelength file: ");
	scanf("%s", &wl_file);
	printf("Wavelength Increment (nm): ");
	scanf("%f", &wl_inc);
	printf("Time Increment (s): ");
	scanf("%f", &time_inc);
	
	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.%s", a_year[i], a_day[i], inst);		// load SERC data 
		printf("%s\t",in_file);
		
		for(j = 1; j <= n_wl; j++){
			match[j] = 0;
		}
		
		fp_in = fopen(in_file, "r");
		do{
			fscanf(fp_in, "%s", &st);
		} while(strcmp(st, "Number_of_columns") != 0);
		fscanf(fp_in, "%d", &n_col);
		do{
			fscanf(fp_in, "%s", &st);
		} while(strcmp(st, "Data:") != 0);
		while(fscanf(fp_in, "%f", &time) != EOF){
			for(j = 1; j <= n_wl; j++){
				fscanf(fp_in, "%f", &float_in);
				if(time >= (a_hour[i] + (wl[j] - 290.0)*time_inc/3600.0/wl_inc)){
					if(match[j] == 0){
						Conv[j] = float_in;
						match[j] = 1;
					}
				}
			}
			for(j = n_wl+2; j <= n_col; j++){	
				fscanf(fp_in, "%f", &float_in);
			}
		}
		fclose(fp_in);
		
		/*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;
}