// Code for reading in the fetch and depth data (raw files from pre_proc.cpp)
// and modifying according to the averaging schemes (IDW depth and cosine weighted fetch)
// based on the parameters entered here. 
// This code eliminates having to re-run pre_proc.cpp for updating averaging parameters.  
// The final result is 3 modified files (depth, fetch, and bottom slope) that are used with 
// param_wave.cpp to compute wave heights.
//
// By: Samuel Carter Boyd

#include <cmath>  
#include <iostream>
#include <iomanip>
#include <fstream>
#include <string>
#include <stdlib.h>
#include <cstdlib>
#include <math.h> 
#include <vector>
#include <sstream>
#include <algorithm>
#include <ctime>
#include <stdio.h>
#include <string.h>

using namespace std;


int line=0;
int x=0;
int x1=0;
int x2=0;
int node_tot,elm_tot;
char content[10],temp[10];
string s; 
int dir_tot,node,k3,k1,k2,direction,num_depths;

int Angres;
//const int Angres=2;

double Eff_top,Eff_bot,Eff_fetch,depth_temp,depth_tot;
double pi = 4*atan(1);

bool check;
bool check2=true;

double p=1;
double p_temp;
double Up;
double Eff_res;
double depth_cur,depth_prev,denominator,Eff_depth;
int res;

int i,n,di_low,di_high,c,c2;
double Hw,Hwabove,Hwbelow,u,di_temp,di_input,dir,y3,PRN, HWEMO,HSPM,HCEM,HTMA,WSPM,TSPM,TCEM,TTMA,Tw;

double WEMO_shoal,WEMO_fric,WEMO_cum,H_break;
double kSPM,LSPM,ksSPM,kfSPM,SPM_shoal,SPM_fric,SPM_cum;
double kCEM,LCEM,ksCEM,kfCEM,CEM_shoal,CEM_fric,CEM_cum;
double kTMA,LTMA,ksTMA,kfTMA,TMA_shoal,TMA_fric,TMA_cum;

double sigma;
double tol;
double eps;
double knew,kold,k0;
double k;
double L;
double ks,kf,Tm,fm;


double g=9.7918;	//gravity at meanlat (28.1091 deg)

int NDSETSE,NP,DTDP,NSPOOLGE,IRTYPE,TIME,IT;

int Ftop,Ftop2,Favg=0;
double Dtop,Dtop2,Davg,dconst,Mtop,Mtop2,Mavg=0;


int main(){

	ifstream infile;
	infile.open("fort.14");			//name of grid file

	if(infile.is_open()){

		cout<<endl;
		cout<<"Reading domain (fort.14 open)"<<endl;

	}

	else{

		cout<<"!!!ERROR domain not read (fort.14 not open)!!!"<<endl;
	}

	
		while(infile >> content){		//reading element and node total from fort.14		
	
			x=x+1;

			if(x==2){	

				elm_tot = atoi(content);
			
			}

			if(x==3){

				node_tot = atoi(content);
				break;

			}

		}

	x=0;

	infile.close();

	cout<<endl;
	cout<<"Element Total: "<<elm_tot<<"          Node Total: "<<node_tot<<endl;
	cout<<endl;


cout<<"Enter Angular Resolution used to generate distance and depth data (deg): ";
cin>>Angres; 
cout<<endl;

cout<<"Enter Step Resolution used to generate distance and depth data (m): ";
cin>>res; 


double fetch[node_tot][360/Angres];
double F[node_tot][360/Angres];
double depth[node_tot][360/Angres];
double d[node_tot][360/Angres];
double slope[node_tot][360/Angres];	
//double p[node_tot][360/Angres];	

double WVX[node_tot];
double WVY[node_tot];
double WDIR[node_tot];
double WVEL[node_tot];

	ifstream infile1;
	infile1.open("Distance_data_2_deg.txt");		//reading in fetch data

	float clock_t,time_req;
	time_req=clock();			//initializing the clock for timing tasks

	if(infile1.is_open()){

		cout<<"Reading Fetch Data"<<endl;
		cout<<endl;
	
		while(getline(infile1,s)){

			line=line+1;
			
			if(line>3 && line<126773){	//ignoring header lines

				if(s.length()>7){	//if the length of line > 7 characters

					stringstream ss;
					ss<<s;			//read the line as string (s)
					ss>>content;		//read number by number from s as (content)
					x=0;			//indexing variable reset to 0

					while(!ss.eof()){	//while values in the line (string ss)

						ss>>temp;	//read each value
						x=x+1;		//indexing variable 		

						if(x==1){	//if first number, store value as node number

							node=atof(content);
							//cout<<"NODE: "<<node<<endl;

						}

							//if any of the next 360/Angres numbers, store as fetch
						if(x<1+(360/Angres)){

							fetch[node][x]=atof(temp);
							//cout<<fetch[node][x]<<endl;

						}

					}	//close while !eof for ss

				}	//close if s.length>7
	
			}	//close if header >3 and <126773

		}	//close while getline for "Distance_data4.txt"

	}	//close if "Distance_data4.txt" is open

	else{

		cout<<"!!! ERROR Distance_data_2_deg.txt not read !!!"<<endl;

	}

// AVG fetch


	//ofstream outfile4("fetch_dist.csv");
/*
	for(int j=1; j<=node_tot; j++){	

		Ftop=0;	

		for(int i=1; i<=360/Angres; i++){

			Ftop=Ftop+fetch[j][i];

			if(fetch[j][i]!=0){
				//outfile4<<fetch[j][i]<<",";
			}

		}

		Ftop=Ftop/(360/Angres);
		Ftop2=Ftop2+Ftop;
		c=c+1;

	}

	//outfile4.close();

	Favg=Ftop2/c;

	cout<<"Average Fetch in domain: "<<Favg<<" m"<<endl;
*/



	cout<<"Enter a positive power parameter for depth weights"<<endl;
	cout<<"(100 = no depth weighting, 0 = variable weighting): ";
	cin>>p_temp;

// VARIABLE WEIGHTING SHOULD NOT BE USED UNTIL A BETTER P-VALUE RELATIONSHIP IS DETERMINED

	check=false;

	if(p_temp>0){

		check=true;
		cout<<endl;

	}

	if(p_temp==0){		//variable weighting 0 number

		check=true;
		cout<<endl;

		cout<<"WARNING: variable p-value should not be used until"<<endl;
		cout<<"         a better p-value relationship is determined"<<endl;
		cout<<endl;
		cout<<"Enter Avg. wind speed expected in domain (m/s)"<<endl;
		cout<<"(this will be used for variable IDW p-value): ";
		cin>>Up;
		cout<<endl;

	}

	while(check==false){

		cout<<"Inter valid power parameter (greater than or equal to zero): ";
		cin>>p_temp;

		if(p_temp>0){

			check=true;
			cout<<endl;

		}

	}



	cout<<"Enter an angular resolution for effective fetch weighting (deg): ";
	cin>>Eff_res;
	cout<<"Effective angular averaging will consider "<<floor(Eff_res/Angres)<<" rays on either side"<<endl;
	cout<<endl;


//****************Reading in Depth Data

	line=0;

	ifstream infile2;
	infile2.open("Depth_data_2_deg.txt");		//reading in depth data


//	if(p==100){
//
//		cout<<"should skip if p=100"<<endl;

//	}

	if(infile2.is_open()){

		cout<<"Reading Depth Data and Computing Averages with p: "<<p_temp<<endl;
	
		while(getline(infile2,s)){

			line=line+1;

			if(line==2){	//second line is where node total and direction total are stored

				stringstream ss;
				ss<<s;
				ss>>content;
				node_tot=atof(content);
				//cout<<"Node_tot: "<<node_tot<<endl;

				ss>>temp;
				dir_tot=atof(temp);
				//cout<<"Dir_tot: "<<dir_tot<<endl;

			}

//x2 indexes each line of Depth_data from 0 to num_depths and resets after every direction
//x1 indexes the number of values per line after the first value

			if(line>2){	//for all other lines

				x2=x2+1;		//indexing variable 2 progress
				x1=0;			//indexing variable 1 reset when there is a new line
				stringstream ss;
				ss<<s;
				ss>>content;

				while(!ss.eof()){

					x2=0;		//indexing varriable 2 reset when there is a new value
					ss>>temp;	//read in each value of ss
					x1=x1+1;	//progress indexing variable 1
					check=false;	//reset boolean check value to false	

					if(x1==1){	//first value is direction (1:4*Ares)

						direction=atof(temp);
						check=true;	//if the direction value is read, check is reset 							         //to true because you know a new node has been reached

					}

					if(check=true && depth_tot!=0){		//if a new node has been reached 

						direction=direction-1;		//reset direction to the previous

						if(direction==0){		//if==0 reset to the last value

							direction=360/Angres;

						}

					//depth averaging scheme

						d[node][direction]=depth_tot/denominator;

						//cout<<"p["<<node<<"]["<<direction<<"]: "<<p<<endl;

				//cout<<"Depth_avg["<<node<<"]["<<direction<<"]: "<<d[node][direction]<<endl;

						depth_tot=0;	//resetting depth_tot to 0 once value is stored
						denominator=0;

					}

					if(x1==2){	//second value is number of depths(needed for averaging)

						num_depths=atof(temp);

					}

				}	//close while !ss.eof

				if(x2==0){		//reading in node number
				
					node=atof(content);
					//cout<<"N: "<<node<<endl;
					

				}

				if(x2==1){		//first depth value special case
					
					depth_temp = atof(content);

//****************************************************************************************
//variable p-value, calculated every new fetch angle, constant along that single fetch direction
//Equation came from exponential fitting the relationships for development curves
//****************************************************************************************

					p=p_temp;

					if(p_temp==0){

					p=log((10924*exp(-0.064*Up))/depth_temp) / log(fetch[node][direction]);

					}
					
					if(p < 0){	//occurs when depth is set to -99999

						if(p_temp==0){

					p=log((10924*exp(-0.064*Up))/depth_temp) / log(fetch[node][direction]);

						}

						if(p_temp!=0){

							p=p_temp;
						}

					//cout<<"p<0 at ["<<node<<"]["<<direction<<"]"<<endl;

					}

					depth_tot = (2*depth_temp)/(pow(res,p));
					denominator = 2/(pow(res,p));
//cout<<"depth["<<x2<<"]: "<<depth_temp<<" depth_tot: "<<depth_tot<<" denominator: "<<denominator<<endl;
					
				}

			if(p==0){
				
				//if(check2==true){
					//cout<<"Raw depth values being used (no IDW) at node: "<<node<<endl;
					//check2=false;
				//}

				if(x2==1 && direction==1){

					dconst=atof(content);
					//cout<<"depth["<<node<<"]: "<<depth_temp<<endl;

				}

				d[node][direction]=dconst;
				//cout<<"d["<<node<<"]["<<direction<<"]: "<<d[node][direction]<<endl;

			}

				
					
				if(x2>2){			//skip 2nd line, built into special case, read 3rd line+ in the depths and sum together for average
					
					depth_temp = atof(content);

					depth_tot = depth_tot+(depth_temp/(pow(((x2-1)*res),p)));
					denominator = denominator+(1/(pow(((x2-1)*res),p)));
//cout<<"depth["<<x2<<"]: "<<depth_temp<<" depth_tot: "<<depth_tot<<" denominator: "<<denominator<<endl;

				}


			}	//close if line > 2

		}	//close while getline

	}	//close if infile2 is_open

	else{

		cout<<"!!! ERROR Depth_data_2_deg.txt not read !!!"<<endl;

	}




//avg bottom slope

	ifstream infile4;
	infile4.open("Slope_upwind_2_deg.txt");		//reading in bottom slope data

	x=-1;
	n=0;

	if(infile4.is_open()){

		cout<<"Reading bottom slope"<<endl;
		cout<<endl;

		while(infile4 >> content){

			x=x+1;

			if(x%(dir_tot+1)==0){		

				n=n+1;				//node counter increment
				i=-1;				//direction counter reset at new node 

			}

			i=i+1;					//direction counter increment
			slope[n][i]=atof(content);		//fetch values

			//if(n>6 && n<9){

				//cout<<"slope["<<n<<"]["<<i<<"]: "<<slope[n][i]<<" x: "<<x<<endl;

			//}
	
		}

	}

	else{

		cout<<"!!! ERROR Slope_upwind_2_deg.txt not read !!!"<<endl;

	}

	infile4.close();



//*****Avg fetch and avg depth calculations

	//ofstream outfile4("fetch_dist.csv");
	//ofstream outfile2("IDW_depth_check.txt");


//*****Avg fetch, depth, and slope calculations


	for(int j=1; j<=node_tot; j++){	

		Ftop=0;	
		Dtop=0;
		Mtop=0;
		c2=0;			//number of non-zero values per node (0-180) if c2=0 then node is dry

		for(int i=1; i<=360/Angres; i++){


			if(d[j][i]>0){

				Ftop=Ftop+fetch[j][i];
				Dtop=Dtop+d[j][i];

				if(slope[j][i] < 1){
			
					Mtop=Mtop+slope[j][i];
				}
	
				c2=c2+1;

				if(j>6 && j<10){

					//cout<<"slope["<<j<<"]["<<i<<"]: "<<slope[j][i]<<" Mtop: "<<Mtop<<endl;
					//cout<<"fetch["<<j<<"]["<<i<<"]: "<<fetch[j][i]<<endl;
					//cout<<"    d["<<j<<"]["<<i<<"]: "<<d[j][i]<<endl;

				}

			}

		}
		

		if(c2!=0){		//wet node	

			c=c+1;		//number of non-zero value nodes (wet nodes considered for avg)

			//cout<<"!!!non- zero node"<<endl;

			Ftop=Ftop/c2;
			Dtop=Dtop/c2;
			Mtop=Mtop/c2;
			Ftop2=Ftop2+Ftop;
			Dtop2=Dtop2+Dtop;
			Mtop2=Mtop2+Mtop;

		}

	}


	Favg=Ftop2/c;
	Davg=Dtop2/c;
	Mavg=Mtop2/c;

	cout<<"Average Fetch in domain: "<<Favg<<" m"<<endl;
	cout<<"Average Depth in domain: "<<Davg<<" m"<<endl;
	cout<<"Average slope in domain: "<<Mavg<<" = 1/"<<floor(1/Mavg)<<endl;
	cout<<endl;




//******Effective fetch calculation*******

//	ofstream outfile1("Eff_fetch_1_40.txt");

	cout<<"Calculating Effective Fetch with res: "<<Eff_res<<" deg"<<endl;


	for(int j=1; j<=node_tot; j++){			//for all nodes

		//cout<<"Node(j): "<<j<<endl;

//		outfile1<<j<<"  ";


		for(int i=1; i<=360/Angres; i++){	//for all angles

			k2=0;				//reset all counting and temp variables
			k1=0;
			Eff_top=0;
			Eff_bot=0;

			//cout<<"Angle(i): "<<i<<endl;

			if(Eff_res==0){			//case for no effective resolution (use raw fetches)

				F[j][i]=fetch[j][i];
		
//				outfile1<<F[j][i]<<" ";

			}

			else{				//case for any effective resolution

				F[j][i]=0;

//rounding down to the nearest integer division for direction given the Effective fetch resolution specified
			for(int n=-floor(Eff_res/Angres); n<=floor(Eff_res/Angres); n++){


				//cout<<"k0: "<<i+n<<endl;
				check=false;				//resetting special case check value

	//case for start nodes (1:4) when <0 is encountered

				if(i+n < 1 && i+n+360/Angres>k1){//direction within Eff_res deg left or right of 								//the current node. for values that go below 1

					k1=i+n+(360/Angres);	//reset index to (33:36)
					//cout<<"k1: "<<k1<<endl;
					
					//cosine weighted average numerator and denominator, as summations

					Eff_top = fetch[j][k1]*cos(n*Angres*pi/180) + Eff_top;		
					Eff_bot = cos(n*Angres*pi/180) + Eff_bot;

					check=true;		//checks for this special case

				}

	//case for start nodes (33:36) when >36 is encountered

	//for values that go above max direction

				if(i+n > 360/Angres && i+n-360/Angres>k2){

					k2=i+n-(360/Angres);	//reset index to (1:4)
					//cout<<"k2: "<<k2<<endl;

					Eff_top = fetch[j][k2]*cos(n*Angres*pi/180) + Eff_top;
					Eff_bot = cos(n*Angres*pi/180) + Eff_bot;

					check=true;		//checks for this special case

				}

	//case for all other nodes(normal +/- Angres deg on each side of start node				

				if(check==false){		//if neither special case, than regular case

					k3=i+n;
					//cout<<"k : "<<k<<endl;

					Eff_top = fetch[j][k3]*cos(n*Angres*pi/180) + Eff_top;
					Eff_bot = cos(n*Angres*pi/180) + Eff_bot;

				}

			}	//close for -4<n<4 loop

			}	//close else statement

			if(Eff_res!=0){		

				Eff_fetch = Eff_top/Eff_bot;	//devide sum top/sum bot for Effective fetch

				F[j][i]=Eff_fetch;
		
//				outfile1<<F[j][i]<<" ";

			}

		//cout<<"fetch: "<<fetch[j][i]<<"  Eff_fetch: "<<Eff_fetch<<endl;

		}	//close for 1<i<36 loop

//		outfile1<<endl;

	} 	//close for all nodes loop


//storing depth, fetch, and fetch MODified avgs in output files

	ofstream outfile1("IDW_depth_MOD.txt");
	ofstream outfile2("Eff_fetch_MOD.txt");

	outfile2<<node_tot<<"  "<<Angres<<endl;		//header info for wave4_07

	for(int l=1; l<=node_tot; l++){

		outfile1<<l<<"  ";
		outfile2<<l<<"  ";

		for(int j=1; j<=360/Angres; j++){

			outfile1<<d[l][j]<<" ";			
			outfile2<<F[l][j]<<" ";			

		}

		outfile1<<endl;
		outfile2<<endl;

	} 

	outfile1.close();
	outfile2.close();


	
	time_req=clock()- time_req;

	cout<<"***Total time taken: "<<(float)time_req/CLOCKS_PER_SEC<<" seconds***"<<endl;

return 0;
}