#include <stdio.h>
#include <math.h>
#include <stddef.h>
#include <stdlib.h>
#include <time.h>
/*#include "nr.h"*/
#include "nrutil.h"
/*#include "nrfitting.c"*/

#define PI 3.14159

int solar_angles(float J_day, float gmt, float lat, float lon, float *zenith, float *azimuth)

{

	double day_angle, declination, eqn_time;
	double time_zone, std_long, local_time, std_time;
	double hour_angle, altitude;
	double int_part, fr_part;
	
	day_angle = 2.0*PI*(J_day-1)/365.0;
	declination = (0.006918 - 0.399912*cos(day_angle) + 0.070257*sin(day_angle)
		- 0.006758*cos(2.0*day_angle) + 0.000907*sin(2.0*day_angle)
		- 0.002697*cos(3.0*day_angle) + 0.00148*sin(3.0*day_angle));
	eqn_time = (0.000075 + 0.001868*cos(day_angle) - 0.032077*sin(day_angle) 
		- 0.014615*cos(2.0*day_angle) - 0.04089*sin(2.0*day_angle))*229.18/60.0;
	time_zone = lon/15.0;
	fr_part = modf(time_zone, &int_part);
	if(fr_part < 0.5) time_zone = int_part;
	else time_zone = int_part + 1;
	std_long = 15.0*time_zone;
	local_time = gmt - time_zone;
	std_time = local_time + 4.0*(std_long - lon)/60.0 + eqn_time;
	hour_angle = (12.0 - std_time)*15.0*PI/180.0;
	*zenith = acos(sin(declination)*sin(lat*PI/180.0) + cos(declination)*
		cos(lat*PI/180.0)*cos(hour_angle))*180.0/PI;
	altitude = (90.0 - *zenith)*PI/180.0;
	*azimuth = acos((sin(altitude)*sin(lat*PI/180.0) - sin(declination))/
		(cos(altitude)*cos(lat*PI/180.0)))*180.0/PI;
	if(hour_angle < 0) *azimuth *= -1.0;
	
	return 1;
}

/* --------------------------------------------------------------------- */

int conv_spline(int n_wl, float *wl, int n_irrad, 
	float **m_irrad, float *conv)
{
	int i, j;
	float *irrad_wl, *irrad;
	float irrad_fit;
	float *irrad2;
	char st[6];
	
	irrad_wl = vector(1, n_irrad);
	irrad = vector(1, n_irrad);
	irrad2 = vector(1, n_irrad);
	
	for(i = 1; i <= n_irrad; i++){
		irrad_wl[i] = m_irrad[i][1];
		irrad[i] = m_irrad[i][2];
	}
	
	spline(irrad_wl, irrad, n_irrad, 1e40, 1e40, irrad2);
	
	for(i = 1; i <= n_wl; i++){
		splint(irrad_wl, irrad, irrad2, n_irrad, wl[i], &irrad_fit);
		conv[i] = irrad_fit;
	}
	
	free_vector(irrad_wl, 1, n_irrad);
	free_vector(irrad, 1, n_irrad);
	free_vector(irrad2, 1, n_irrad);
	
	return 1;
}

/* --------------------------------------------------------------------- */

int conv_tri(float bw, int n_wl, float *wl, int n_irrad, 
	float **m_irrad, float *conv)
{
	int i, j;
	double sum, prod, R_unc, S_unc;
	float tri;
	
	for(i = 1; i <= n_wl; i++){
		sum = 0, prod = 0, R_unc = 0, S_unc = 0;
		for(j = 1; j <= n_irrad; j++){
			if(m_irrad[j][1] >= (wl[i] - bw)){
				if(m_irrad[j][1] > (wl[i] + bw)) break;
				if(m_irrad[j][1] <= wl[i]){
					tri = (m_irrad[j][1] - wl[i])/bw + 1;
					sum += tri;
					prod += tri*m_irrad[j][2];
					R_unc += SQR(tri*m_irrad[j][2]*m_irrad[j][3]);
					S_unc += tri*m_irrad[j][2]*m_irrad[j][4];
				}
				if(m_irrad[j][1] > wl[i]){
					tri = (wl[i] - m_irrad[j][1])/bw + 1;
					sum += tri;
					prod += tri*m_irrad[j][2];
					R_unc += SQR(tri*m_irrad[j][2]*m_irrad[j][3]);
					S_unc += tri*m_irrad[j][2]*m_irrad[j][4];
				}
			}
		}
		conv[i] = prod/sum;
	}
	
	return 1;
}

/* --------------------------------------------------------------------- */

int conv_rect(float bw, int n_wl, float *wl, int n_irrad, 
	float **m_irrad, float *conv)
{
	int i, j;
	double sum, prod, R_unc, S_unc;
	
	for(i = 1; i <= n_wl; i++){
		sum = 0, prod = 0, R_unc = 0, S_unc = 0;
		for(j = 1; j <= n_irrad; j++){
			if(m_irrad[j][1] >= (wl[i] - bw/2.0)){
				if(m_irrad[j][1] > (wl[i] + bw/2.0)) break;
				sum += 1;
				prod += m_irrad[j][2];
				R_unc += SQR(m_irrad[j][2]*m_irrad[j][3]);
				S_unc += m_irrad[j][2]*m_irrad[j][4];
			}
		}
		conv[i] = prod/sum;
	}
	
	return 1;
}
	
/* --------------------------------------------------------------------- */

int conv_gauss(float bw, int n_wl, float *wl, float cutoff, 
	int n_irrad, float **m_irrad, float *conv)
{
	int i, j;
	double sum, prod, R_unc, S_unc;
	float gauss;
	
	for(i = 1; i <= n_wl; i++){
		sum = 0, prod = 0, R_unc = 0, S_unc = 0;
		for(j = 1; j <= n_irrad; j++){
			if(fabs(wl[i] - m_irrad[j][1]) <= 1.7*bw) gauss = exp(-1*SQR((m_irrad[j][1] - wl[i])/bw));
			else gauss = 0.00001;
			sum += gauss;
			prod += gauss*m_irrad[j][2];
			R_unc += SQR(gauss*m_irrad[j][2]*m_irrad[j][3]);
			S_unc += gauss*m_irrad[j][2]*m_irrad[j][4];
			if(m_irrad[j][1] > cutoff) break;
		}
		conv[i] = prod/sum;
	}
	
	return 1;
}
				
/* ------------------------------------------------------------------------ */

int conv_filt(float **filter, float *wl, 
	int n_filt_row, int n_filt_col, int n_wl, int n_irrad, 
	float **m_irrad, float *conv)
{
	int i, j;
	double sum, prod;
	float spl_fit, wl_fit;
	float *filt, *filt_wl, *filt2;
	
	filt_wl = vector(1, n_filt_row);
	filt = vector(1, n_filt_row);
	filt2 = vector(1, n_filt_row);
	
	for(i = 1; i <= n_wl; i++){
		for(j = 1; j <= n_filt_row; j++){
			filt_wl[j] = filter[j][1];
			filt[j] = filter[j][i+1];
			if(filt[j] <= 0.00001) filt[j] = 0.00001;
		}
		spline(filt_wl, filt, n_filt_row, 1e40, 1e40, filt2);
		
		sum = 0, prod = 0;
		for(j = 1; j <= n_irrad; j++){
			wl_fit = m_irrad[j][1];
			if((wl_fit >= filt_wl[1]) && (wl_fit <= filt_wl[n_filt_row])){
				splint(filt_wl, filt, filt2, n_filt_row, wl_fit, &spl_fit);
				sum += spl_fit;
				prod += spl_fit*m_irrad[j][2];
			}
		}
		printf("prod = %f\tsum = %f\n", prod, sum);
		conv[i] = prod/sum;
	}
	
	free_vector(filt_wl, 1, n_filt_row);
	free_vector(filt, 1, n_filt_row);
	free_vector(filt2, 1, n_filt_row);

	return 1;
}


/* ------------------------------------------------------------------------ */

int store_header(FILE *fp_out, char time_file[], int conv_fnc, float bw, 
	float min_wl, float max_wl, float inc_wl, int n_time, 
	int day[], float time[], float zenith[], float azimuth[])
{
	int i;
	char st[50];
	
	fprintf(fp_out, "Time File:\t%s\n", time_file);
	if(conv_fnc == 1) sprintf(st, "Spline Fit");
	if(conv_fnc == 2) sprintf(st, "Triangle Convolve");
	if(conv_fnc == 3) sprintf(st, "Rectangle Convolve");
	if(conv_fnc == 4) sprintf(st, "Gaussian Convolve");
	if(conv_fnc == 5) sprintf(st, "Filter Convolve");
	fprintf(fp_out, "Analysis:\t%s\n", st);
	if((conv_fnc != 1) &&(conv_fnc != 6)) fprintf(fp_out, "\t%.2f\tnm bandwidth\n", bw);
	if(conv_fnc != 5) fprintf(fp_out, "\t%.2f\tto\t%.2f\tnm by\t%.2f\tnm\n", min_wl, max_wl, inc_wl);
	fprintf(fp_out, "\t%s\n\n", st);
	
	fprintf(fp_out, "Day:\t");
	for(i = 1; i <= n_time; i++) fprintf(fp_out, "%d\t", day[i]);
	fprintf(fp_out, "\n");
	fprintf(fp_out, "Time:\t");
	for(i = 1; i <= n_time; i++) fprintf(fp_out, "%.1f\t", time[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");
	
	return 1;
}