// This program pre-processes the fetch and depth data needed for the parametric wave code (param_wave.cpp).
// It reads in a grid file (fort.14) and calculats fetch rays (cosine averaged), depth values (I.D.W. averaged), 
// and bottom slopes for all wet nodes in the domain at angular increments (2 deg reccomended)
//
// By: Samuel Carter Boyd
// Date: 04/13/2020


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


using namespace std;

char content[10];
double a;
double distlon,distlat;
double lat0,lon0,dlat1,dlat2,dlat3,dlat4,dlon1,dlon2,dlon3,dlon4;
int wet_tot,dry_tot,x,i,c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,ii,step1,step2,step3,step4,node1,node_tot,elm_tot;

double res;	//resolution (meters) for spatial step
int Angres;	//res (deg) for angular step
bool check;
double p,p_temp,Eff_res;	//IDW power value and cosine weighted angular averaging value

double vx[4] = {};
double vy[4] = {};
double dx;		//deg lat
double dy;		//deg lon
double xloc[4],yloc[4] = {};
double D,delx,dely,theta;
double crossprod1,crossprod2,crossprod3;
int contELM1,contELM2,contELM3,contELM4;
double w1,w2,w3;


vector<double> dephold1,dephold2,dephold3,dephold4;

double meanwetlat,meanwetlon,wetlonsum,wetlatsum = 0;
double pi = 4*atan(1);
double Rearth;
double Req=6378137;	
double Rpo=6356752;
double dhav1,dhav2,dhav3,dhav4;
int emin1,emax1,emin2,emax2,emin3,emax3,emin4,emax4=0;


int line,x1,x2,c,Ftop,Ftop2,Favg,k,k1,k2,num_depths,direction,dir_tot=0;
string s;
char temp[10];
double g=9.7918;
double Eff_top,Eff_bot,Eff_fetch,depth_temp,depth_tot,depth_cur,depth_prev,Eff_depth,denominator,Up;
int Ftop,Ftop2,Favg=0;
double Dtop,Dtop2,Davg,dconst,Mtop,Mtop2,Mavg=0;
	

int main(){

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

//************************************************************************************************
//These are the names of the final output files created for the grid's depth and fetch information
//************************************************************************************************

	ofstream outfile6("Distance_data_2_deg.txt");
	ofstream outfile7("Depth_data_2_deg.txt");
	ofstream outfile8("Slope_upwind_2_deg.txt");
	ofstream outfile9("p_values_2_deg.txt");
	ofstream outfile2("IDW_depths_2_deg.txt");
	ofstream outfile1("Eff_fetch_2_deg.txt");
	ofstream outfile5("Raw_depth_2_deg.txt");


	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;
	}

						//reading element and node total from fort.14
	
		while(infile >> content){			
	
			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;
			
					//defining arrays once dimensions have been read
double elv[node_tot+1];		//all nodes elevation
double lon[node_tot+1];		//all nodes longitude
double lat[node_tot+1];		//all nodes node latitude
double node[node_tot+1];	//all nodes node number
double wetnode[node_tot+1];	//wet node number
double wetelv[node_tot+1];	//wet node elevation
double wetlon[node_tot+1];	//wet node longitude
double wetlat[node_tot+1];	//wet node latitude
int ELM[elm_tot];		//triangular elements
int NE[elm_tot][3];		//neighboring nodes in each element

	infile.open("fort.14");
	
					//reading lat,lon,elv,node,ELM,and NE data from fort.14
						
		while(infile >> content){

			x=x+1;

			if(x>=4 && x<= 4 + (node_tot*4)){	//4 is from node number(1),lat(2),lon(3),elv(4)
					
				i = i+1;
				a = atof(content);

				if(((i+0)%4) == 0){
					c1=c1+1;
					elv[c1] = a;
				}

				if(((i+1)%4) == 0){
					c2=c2+1;
					lat[c2] = a;
				}

				if(((i+2)%4) == 0){
					c3=c3+1;
					lon[c3] = a;
				}

				if(((i+3)%4) == 0){
					c4=c4+1;					
					node[c4] = a;
				}
			}

//Element information is being read once code reaches end of regular nodal info
 
			if(x >= (4 + (node_tot*4)) && x < (4 + (node_tot*4))+(elm_tot*5)){

				ii = ii+1;
				a = atoi(content);

				if(((ii+4)%5) == 0){
					c7 = c7+1;
					ELM[c7] = a;
				}

				if(((ii+2)%5) == 0){
					c8 = c8+1;
					NE[c8][1] = a;
				}

				if(((ii+1)%5) == 0){
					c9 = c9+1;
					NE[c9][2] = a;
				}

				if(((ii+0)%5) == 0){
					c10 = c10+1;
					NE[c10][3] = a;
				}

			}

		}	

	infile.close();

//calculating average wet lat and lon in domain
				
		for(int j=1; j<=node_tot; j++){
				
			if(elv[j] > 0){

				c5=c5+1;	

				wetlatsum = wetlatsum+lat[j];
				wetlonsum = wetlonsum+lon[j];
					
			}	

			if(elv[j] < 0){

				c6=c6+1;
				
			}

			
		}	//close node_tot loop


		wet_tot = c5;
		dry_tot = c6;

		meanwetlat = wetlatsum/wet_tot;
		meanwetlon = wetlonsum/wet_tot;
		cout<<"Mean longitude (wet): "<<meanwetlon<<"        mean latitude wet(wet): "<<meanwetlat<<endl;
		cout<<endl;

//interpolating the radius of earth between equator and poles wrt avg wet lattitude of the domain

	Rearth = sqrt((pow((Req*Req*cos(meanwetlat*pi/180)),2.0)+pow((Rpo*Rpo*sin(meanwetlat*pi/180)),2.0))/
		 (pow((Req*cos(meanwetlat*pi/180)),2.0)+pow((Rpo*sin(meanwetlat*pi/180)),2.0)));

	cout<<"Radius of Earth at this location: "<<Rearth/1000<<" km"<<endl;

//haversine forumla distance(m) for 1 deg lat and 1 deg lon based on avg wetlat and wetlon of domain

	distlon = 2*Rearth*asin(sqrt(cos((meanwetlat-.5)*pi/180)*cos((meanwetlat+.5)*pi/180)*
				      sin(.5*pi/180)*sin(.5*pi/180)));
				    	
	distlat = 2*Rearth*asin(sqrt(sin(.5*pi/180)*sin(.5*pi/180)));	

	cout<<"1 deg lattitude = "<<distlat<<" m"<<endl;
	cout<<"1 deg longitude = "<<distlon<<" m"<<endl;
	cout<<endl;

	double d12;
	double d13;
	double d23;
	double dmin;
	double dhold=100000;

//calculating minimum distance (dhold) between wet nodes for spatial step reccomendation

	for(int e=1; e<=elm_tot; e++){

		if(elv[NE[e][1]]>0 && elv[NE[e][2]]>0 && elv[NE[e][3]]>0){

		d12=distlon*sqrt(pow(lon[NE[e][1]]-lon[NE[e][2]],2)+(pow(lat[NE[e][1]]-lat[NE[e][2]],2)));
		d13=distlon*sqrt(pow(lon[NE[e][1]]-lon[NE[e][3]],2)+(pow(lat[NE[e][1]]-lat[NE[e][3]],2)));
		d23=distlon*sqrt(pow(lon[NE[e][2]]-lon[NE[e][3]],2)+(pow(lat[NE[e][2]]-lat[NE[e][3]],2)));

			//cout<<"d12: "<<d12<<" d13: "<<d13<<" d23: "<<d23<<endl;

			if(d12<d13 && d12<d23){
				dmin=d12;	
			}
			if(d13<d12 && d13<d23){
				dmin=d13;	
			}
			if(d23<d12 && d23<d13){
				dmin=d23;	
			}

			//cout<<"dmin: "<<dmin<<endl;

			if(dmin<dhold){
				
				dhold=dmin;
	//cout<<"e: "<<e<<" NE[1]: "<<NE[e][1]<<" NE[2]: "<<NE[e][2]<<" NE[3]: "<<NE[e][3]<<" dmin: "<<dhold<<endl;
			}

		}

		dmin=0;

	}

double steprec=floor(3*dhold);	//recomended spatial step resolution

//User input for step and angular resolution

	cout<<"Enter spatial step resolution (recommended = "<<steprec<<"m): ";
	cin>>res;

	dx = res/distlon;	//50m step in lat/lon degrees specified by the step resolution
	dy = res/distlat;

	//cout<<"dx("<<res<<"m): "<<dx<<" deg lon"<<endl;
	//cout<<"dy("<<res<<"m): "<<dy<<" deg lat"<<endl;

	cout<<"Enter angular resolution (2,5,10,18,or 30 deg): ";
	cin>>Angres;

	check=false;

	if(Angres==2 || Angres==5 || Angres==10 || Angres==18 || Angres==30){

		check=true;
		cout<<endl;

	}

	while(check==false){

		cout<<"Inter valid angular resolution (2,5,10,18, or 30 deg): ";
		cin>>Angres;

		if(Angres==2 || Angres==5 || Angres==10 || Angres==18 || Angres==30){

			check=true;
			cout<<endl;

		}

	}
	

//initializing solution array once angular resolution is known	

	double d[node_tot+1][360/Angres];
	double Ares = (360/Angres)/4;		//Num of angles per quadrent
	int Ares2 = Ares;			//cast to integer for indexing array
	double fetch[node_tot+1][360/Angres];
	double F[node_tot+1][360/Angres];
	double depth[node_tot+1][360/Angres];

	double slopehold[node_tot][360/Angres];
	double startelv;

	for (int j=1; j<=node_tot; j++){		//initializing slope array with zeros
		for(int m=1; m<=360/Angres; m++){
			slopehold[j][m]=0;
		}
	}

double Pelv1[1+(360/(4*Angres))],Pelv2[1+(360/(4*Angres))],Pelv3[1+(360/(4*Angres))],Pelv4[1+(360/(4*Angres))];
double XP1[1+(360/(4*Angres))],XP2[1+(360/(4*Angres))],XP3[1+(360/(4*Angres))],XP4[1+(360/(4*Angres))];
double YP1[1+(360/(4*Angres))],YP2[1+(360/(4*Angres))],YP3[1+(360/(4*Angres))],YP4[1+(360/(4*Angres))];
double Phold1[1+(360/(4*Angres))],Phold2[1+(360/(4*Angres))],Phold3[1+(360/(4*Angres))],Phold4[1+(360/(4*Angres))];
double dpast1[1+(360/(4*Angres))],dpast2[1+(360/(4*Angres))],dpast3[1+(360/(4*Angres))],dpast4[1+(360/(4*Angres))];


	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;
	
	//creating raw depth output (just the depth at each node)


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

		outfile5<<i<<"  "<<elv[i]<<endl;

	} 

	outfile5.close();


//***************************************Begin Fetch and Depth Code****************************************

	cout<<"**********Calculating Fetch and Depth**********"<<endl; 

	outfile6<<"Fetch Distances at spatial res: "<<res<<"[m] and angular res: "<<Angres<<"[deg]"<<endl;
	outfile6<<"Grid used to generate file: 'fort.14'"<<endl;

	outfile7<<"Water Depths along fetch rays from 'fort.14'"<<endl;
	outfile7<<node_tot<<" "<<360/Angres<<endl;

	//outfile8<<"Upwind bottom slope from node used for breaking criteria"<<endl;
	//outfile8<<"Grid used to generate file: 'fort.14'"<<endl;
	

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

		cout<<"node: "<<j<<endl;
		outfile6<<j<<"  ";
		outfile8<<j<<"  ";

		lon0 = lon[j];
		lat0 = lat[j];

		if(elv[j] < 0){				//if node is dry, write 0 slopes for output file

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

				outfile8<<"0 ";

			}
		}

		if(elv[j] > 0){				//if node is wet

			startelv=elv[j];		//depth at the start node

//*************************************Q1********************************************************************

			for(int m=0; m<Ares; m++){   //Angular resolution including 0 excluding Ares [0-90) deg

				outfile7<<j<<" ";
				outfile7<<m+1<<" ";
				dephold1.clear();	//erases elements from dephold vector
	
				step1 = 0;
				Pelv1[m] = elv[j];
				Phold1[m] = elv[j];
				dpast1[m]=0;

							//min and max from elm search optimization analysis

				emin1=2*j-floor(.0163*j+1500);
				emax1=2*j+1500;

				while(Pelv1[m] > 0){	//while (XP,YP) is wet
				
					dephold1.push_back(Pelv1[m]);	//stores every Pelv value
						
					step1=step1+1;		//increment XP and YP 1 step further

					XP1[m] = lon0 + (dx*cos((pi/180)*(90-(Angres*m))))*step1;
					YP1[m] = lat0 + (dy*sin((pi/180)*(90-(Angres*m))))*step1;

//if first step already taken, compare the 1st neighboring node of the previous containing element 
//so that the restricted element search is wrt the containing element instead of the start node j


//Pelv1 at step1==2 corresponds to the first step's interpolated depth. Calculate slope from start to here
//BUT slope values are stored in slopehold array in flipped indexes. Q1<-->Q3 and Q2<-->Q4 because we want
//the index to correspond to upwind slope for determining wave breaking

					if(step1==2){

						slopehold[j][m+1]=(startelv-Pelv1[m])/res;

//if slope is less than 1/1000 set to 0, this inludes negative slopes (from shallow to deeper water)

						if(slopehold[j][m+1] < 1/1000){ 

							slopehold[j][m+1]=0;

						}

						outfile8<<slopehold[j][m+1]<<" ";

					}

	

					if(step1>1){

						emin1=2*NE[contELM1][1]-floor(.0163*NE[contELM1][1]+1500);
						emax1=2*NE[contELM1][1]+1500;

					}

					if(emin1<1){
							
						emin1=1;

					}

					if(emax1>=elm_tot){
							
						emax1=elm_tot;

					}

					contELM1 = false;	//resent containing element

//****Element Search*****

//restricted element search according to emin and emax values

					for(int e=emin1; e<=emax1; e++){

						for(int i=1; i<=3; i++){ 	//vectors of neighboring nodes

							xloc[i] = lon[NE[e][i]];
							yloc[i] = lat[NE[e][i]];
							delx = xloc[i] - XP1[m];
							dely = yloc[i] - YP1[m];
							D = pow((pow(delx,2.0)+pow(dely,2.0)),0.5);
							theta = atan2(dely,delx);    			
							vx[i] = D*cos(theta);
							vy[i] = D*sin(theta);

						}
								//cross products of neighboring node vectors
							
						crossprod1 = vx[1]*vy[2] - vy[1]*vx[2];
						crossprod2 = vx[2]*vy[3] - vy[2]*vx[3];
						crossprod3 = vx[3]*vy[1] - vy[3]*vx[1];

								//if cross products are all >= 0, containing ELM

						if(crossprod1>=0 && crossprod2>=0 && crossprod3>=0){
			
							contELM1 = ELM[e];

//cout<<"R.SearchFound.   contELM1: "<<contELM1<<"   emin:  "<<emin1<<"   emax: "<<emax1<<endl;

//****Barycentric elevation interpolation****

				w1 =((yloc[2]-yloc[3])*(XP1[m]-xloc[3])+(xloc[3]-xloc[2])*(YP1[m]-yloc[3]))/
				    ((yloc[2]-yloc[3])*(xloc[1]-xloc[3])+(xloc[3]-xloc[2])*(yloc[1]-yloc[3]));

				w2 =((yloc[3]-yloc[1])*(XP1[m]-xloc[3])+(xloc[1]-xloc[3])*(YP1[m]-yloc[3]))/
				    ((yloc[2]-yloc[3])*(xloc[1]-xloc[3])+(xloc[3]-xloc[2])*(yloc[1]-yloc[3]));

				w3 = 1-w1-w2;

				Pelv1[m] = (w1*elv[NE[e][1]]+w2*elv[NE[e][2]]+w3*elv[NE[e][3]])/(w1+w2+w3);
		
							if(Pelv1[m] > 0){
										
								Phold1[m] = Pelv1[m];	//The last wet point

							}

							break;	//out of for elements loop

						}	//close if(crossprod > 0)

					}	//close restricted element search for loop

				if(contELM1 == false){

//cout<<"Long Search Q1************** steps: "<<step1<<endl;

//if containing element is not found, check all elements to see if restricted seach missed it

					for(int e=1; e<elm_tot; e++){

						for(int i=1; i<=3; i++){ 	//vectors of neighboring nodes

							xloc[i] = lon[NE[e][i]];
							yloc[i] = lat[NE[e][i]];
							delx = xloc[i] - XP1[m];
							dely = yloc[i] - YP1[m];
							D = pow((pow(delx,2.0)+pow(dely,2.0)),0.5);
							theta = atan2(dely,delx);    			
							vx[i] = D*cos(theta);
							vy[i] = D*sin(theta);

						}
								//cross products of neighboring node vectors
							
						crossprod1 = vx[1]*vy[2] - vy[1]*vx[2];
						crossprod2 = vx[2]*vy[3] - vy[2]*vx[3];
						crossprod3 = vx[3]*vy[1] - vy[3]*vx[1];

								//if cross products are all >= 0, containing ELM

						if(crossprod1>=0 && crossprod2>=0 && crossprod3>=0){
			
							contELM1 = ELM[e];

//****Barycentric elevation interpolation****

				w1 =((yloc[2]-yloc[3])*(XP1[m]-xloc[3])+(xloc[3]-xloc[2])*(YP1[m]-yloc[3]))/
				    ((yloc[2]-yloc[3])*(xloc[1]-xloc[3])+(xloc[3]-xloc[2])*(yloc[1]-yloc[3]));

				w2 =((yloc[3]-yloc[1])*(XP1[m]-xloc[3])+(xloc[1]-xloc[3])*(YP1[m]-yloc[3]))/
				    ((yloc[2]-yloc[3])*(xloc[1]-xloc[3])+(xloc[3]-xloc[2])*(yloc[1]-yloc[3]));

				w3 = 1-w1-w2;

				Pelv1[m] = (w1*elv[NE[e][1]]+w2*elv[NE[e][2]]+w3*elv[NE[e][3]])/(w1+w2+w3);
		
							if(Pelv1[m] > 0){
										
								Phold1[m] = Pelv1[m];	//The last wet point

							}

							break;	//out of for elements loop

						}	//close if(crossprod > 0)

					}	//close for all elements loop

//cout<<"!!!!contELM1: "<<contELM1<<"   emin:  "<<emin1<<"   emax: "<<emax1<<endl;

				}	//close if contELM1 == false

//if containing element is still false, then the element is out of bounds

					if(contELM1 == false){

						Pelv1[m]=0;	

					}

				}	//close while Pelv (XP,YP) > 0 loop

				outfile7<<dephold1.size()<<endl;

				for(int h=0; h<dephold1.size(); h++){
				
					outfile7<<dephold1[h]<<endl;

				}

					//estimating distance past the last wet point				
	
				dpast1[m] = Phold1[m]/((Phold1[m] + abs(Pelv1[m]))/res); 

				dlon1 = abs(XP1[m] - lon0);	//distance from start to land in lat/lon	
				dlat1 = abs(YP1[m] - lat0);

	//if NOT the first step and does NOT go out of bounds, this is the usual condition
	//distance from start to land-res = dist from start to step before land + dist past that last point 

				dhav1 = 2*Rearth*asin(sqrt(sin(((dlat1)/2)*pi/180)*
							  sin(((dlat1)/2)*pi/180)+
						   (cos(lat0*pi/180)*cos(YP1[m]*pi/180)*
							  sin(((dlon1)/2)*pi/180)*
							  sin(((dlon1)/2)*pi/180))))-res+dpast1[m];


	//if first step hits land, dhav = dpast

				if(step1 == 1){			

					dhav1=dpast1[m];

	//if first step is out of bounds, set dhav = 0

					if(contELM1 == false){	
						
						dhav1=0;

					}

				}

	//if NOT first step and out of bounds,calculate dhav = d(out of bounds)-res/2, this is a round estimate

				if(step1!=1 && contELM1 ==false){  

					dhav1 = 2*Rearth*asin(sqrt(sin(((dlat1)/2)*pi/180)*
								  sin(((dlat1)/2)*pi/180)+
							   (cos(lat0*pi/180)*cos(YP1[m]*pi/180)*
								  sin(((dlon1)/2)*pi/180)*
								  sin(((dlon1)/2)*pi/180))))-(res/2);

				}
			
				outfile6<<dhav1<<" ";

			}	//close for m


//*************************************Q4********************************************************************


			for(int m=0; m<Ares; m++){		//including 0 excluding Ares [90-180) deg

				outfile7<<j<<" ";
				outfile7<<m+Angres<<" ";
				dephold4.clear();

				step4 = 0;
				Pelv4[m] = elv[j];
				Phold4[m] = elv[j];
				dpast4[m]=0;

				emin4=2*j-floor(.0163*j+1500);
				emax4=2*j+1500;

				while(Pelv4[m] > 0){			//while (XP,YP) is wet

					dephold4.push_back(Pelv4[m]);	//stores every Pelv value

					step4=step4+1;

					XP4[m] = lon0 + (dx*cos((pi/180)*(Angres*m)))*step4;
					YP4[m] = lat0 - (dy*sin((pi/180)*(Angres*m)))*step4;

//if first step already taken, compare the 1st neighboring node of the previous containing element 
//so that the restricted element search is wrt the containing element instead of the start node j


//Pelv4 at step4==2 corresponds to the first step's interpolated depth. Calculate slope from start to here
//BUT slope values are stored in slopehold array in flipped indexes. Q1<-->Q3 and Q2<-->Q4 because we want
//the index to correspond to upwind slope for determining wave breaking

					if(step4==2){

						slopehold[j][m+Ares2+1]=(startelv-Pelv4[m])/res;

//if slope is less than 1/1000 set to 0, this inludes negative slopes (from shallow to deeper water)

						if(slopehold[j][m+Ares2+1] < 1/1000){ 

							slopehold[j][m+Ares2+1]=0;

						}

						outfile8<<slopehold[j][m+Ares2+1]<<" ";

					}


					if(step4>1){

						emin4=2*NE[contELM4][1]-floor(.0163*NE[contELM4][1]+1500);
						emax4=2*NE[contELM4][1]+1500;

					}

					if(emin4<1){
							
						emin4=1;

					}

					if(emax4>=elm_tot){
							
						emax4=elm_tot;

					}

					contELM4 = false;	//resent containing element

//****Element Search****

					for(int e=emin4; e<=emax4; e++){

						for(int i=1; i<=3; i++){ 

							xloc[i] = lon[NE[e][i]];
							yloc[i] = lat[NE[e][i]];
							delx = xloc[i] - XP4[m];
							dely = yloc[i] - YP4[m];
							D = pow((pow(delx,2.0)+pow(dely,2.0)),0.5);
							theta = atan2(dely,delx);    			
							vx[i] = D*cos(theta);
							vy[i] = D*sin(theta);

						}

						crossprod1 = vx[1]*vy[2] - vy[1]*vx[2];
						crossprod2 = vx[2]*vy[3] - vy[2]*vx[3];
						crossprod3 = vx[3]*vy[1] - vy[3]*vx[1];

						if(crossprod1>=0 && crossprod2>=0 && crossprod3>=0){
			
							contELM4 = ELM[e];

//****arycentric elevation interpolation****

				w1 =((yloc[2]-yloc[3])*(XP4[m]-xloc[3])+(xloc[3]-xloc[2])*(YP4[m]-yloc[3]))/
				    ((yloc[2]-yloc[3])*(xloc[1]-xloc[3])+(xloc[3]-xloc[2])*(yloc[1]-yloc[3]));

				w2 =((yloc[3]-yloc[1])*(XP4[m]-xloc[3])+(xloc[1]-xloc[3])*(YP4[m]-yloc[3]))/
				    ((yloc[2]-yloc[3])*(xloc[1]-xloc[3])+(xloc[3]-xloc[2])*(yloc[1]-yloc[3]));

				w3 = 1-w1-w2;

				Pelv4[m] = (w1*elv[NE[e][1]]+w2*elv[NE[e][2]]+w3*elv[NE[e][3]])/(w1+w2+w3);
		

							if(Pelv4[m] > 0){
										
								Phold4[m] = Pelv4[m];	//The last wet point 
										

							}

						break;	//out of for elements loop

						}	//close if(crossprod > 0)

					}	//close restricted elemt search for loop

//if no containing element, search all nodes

				if(contELM4 == false){

//cout<<"Long Search Q4************** steps: "<<step4<<endl;

					for(int e=1; e<elm_tot; e++){

						for(int i=1; i<=3; i++){ 

							xloc[i] = lon[NE[e][i]];
							yloc[i] = lat[NE[e][i]];
							delx = xloc[i] - XP4[m];
							dely = yloc[i] - YP4[m];
							D = pow((pow(delx,2.0)+pow(dely,2.0)),0.5);
							theta = atan2(dely,delx);    			
							vx[i] = D*cos(theta);
							vy[i] = D*sin(theta);

						}

						crossprod1 = vx[1]*vy[2] - vy[1]*vx[2];
						crossprod2 = vx[2]*vy[3] - vy[2]*vx[3];
						crossprod3 = vx[3]*vy[1] - vy[3]*vx[1];

						if(crossprod1>=0 && crossprod2>=0 && crossprod3>=0){
			
							contELM4 = ELM[e];

//****Barycentric elevation interpolation****

				w1 =((yloc[2]-yloc[3])*(XP4[m]-xloc[3])+(xloc[3]-xloc[2])*(YP4[m]-yloc[3]))/
				    ((yloc[2]-yloc[3])*(xloc[1]-xloc[3])+(xloc[3]-xloc[2])*(yloc[1]-yloc[3]));

				w2 =((yloc[3]-yloc[1])*(XP4[m]-xloc[3])+(xloc[1]-xloc[3])*(YP4[m]-yloc[3]))/
				    ((yloc[2]-yloc[3])*(xloc[1]-xloc[3])+(xloc[3]-xloc[2])*(yloc[1]-yloc[3]));

				w3 = 1-w1-w2;

				Pelv4[m] = (w1*elv[NE[e][1]]+w2*elv[NE[e][2]]+w3*elv[NE[e][3]])/(w1+w2+w3);
		

							if(Pelv4[m] > 0){
										
								Phold4[m] = Pelv4[m];	//The last wet point 
										

							}

						break;	//out of for elements loop

						}	//close if(crossprod > 0)

					}	//close restricted elemt search for loop
//cout<<"!!!!contELM4: "<<contELM4<<"   emin:  "<<emin4<<"   emax: "<<emax4<<endl;

				}	//close if contELM4 == false


//if contELM4 is still false, then out of bounds

					if(contELM4 == false){

						Pelv4[m]=0;	

					}

				}	//close while Pelv (XP,YP) > 0 loop


				outfile7<<dephold4.size()<<endl;

				for(int h=0; h<dephold4.size(); h++){
				
					outfile7<<dephold4[h]<<endl;

				}

				dpast4[m] = Phold4[m]/((Phold4[m] + abs(Pelv4[m]))/res);

				dlon4 = abs(XP4[m] - lon0);			
				dlat4 = abs(YP4[m] - lat0);


				dhav4 = 2*Rearth*asin(sqrt(sin(((dlat4)/2)*pi/180)*
							  sin(((dlat4)/2)*pi/180)+
						   (cos(lat0*pi/180)*cos(YP4[m]*pi/180)*
							  sin(((dlon4)/2)*pi/180)*
							  sin(((dlon4)/2)*pi/180))))-res+dpast4[m];

				if(step4 == 1){	

					dhav4=dpast4[m];
			
					if(contELM4 == false){

						dhav4=0;

					}

				}

				if(step4!=1 && contELM4 ==false){
		

					dhav4 = 2*Rearth*asin(sqrt(sin(((dlat4)/2)*pi/180)*
								  sin(((dlat4)/2)*pi/180)+
							   (cos(lat0*pi/180)*cos(YP4[m]*pi/180)*
								  sin(((dlon4)/2)*pi/180)*
								  sin(((dlon4)/2)*pi/180))))-(res/2);

				}

				outfile6<<dhav4<<" ";

			}	//close for m


//*************************************Q3********************************************************************
	

			for(int m=0; m<Ares; m++){		//including 0 excluding Ares (0-80 deg)

				outfile7<<j<<" ";
				outfile7<<m+Angres+Angres-1<<" ";
				dephold3.clear();

				step3 = 0;
				Pelv3[m] = elv[j];
				Phold3[m] = elv[j];
				dpast3[m]=0;

				emin3=2*j-floor(.0163*j+1500);
				emax3=2*j+1500;

				while(Pelv3[m] > 0){			//while (XP,YP) is wet

					dephold3.push_back(Pelv3[m]);	//stores every Pelv value

					step3=step3+1;

					XP3[m] = lon0 - (dx*cos((pi/180)*(90-(Angres*m))))*step3;
					YP3[m] = lat0 - (dy*sin((pi/180)*(90-(Angres*m))))*step3;

//if first step already taken, compare the 1st neighboring node of the previous containing element 
//so that the restricted element search is wrt the containing element instead of the start node j

//Pelv3 at step3==2 corresponds to the first step's interpolated depth. Calculate slope from start to here
//BUT slope values are stored in slopehold array in flipped indexes. Q1<-->Q3 and Q2<-->Q4 because we want
//the index to correspond to upwind slope for determining wave breaking

					if(step3==2){

						slopehold[j][m+Ares2+Ares2+1]=(startelv-Pelv3[m])/res;

//if slope is less than 1/1000 set to 0, this inludes negative slopes (from shallow to deeper water)

						if(slopehold[j][m+Ares2+Ares2+1] < 1/1000){ 

							slopehold[j][m+Ares2+Ares2+1]=0;

						}

						outfile8<<slopehold[j][m+Ares2+Ares2+1]<<" ";

					}


					if(step3>1){

						emin3=2*NE[contELM3][1]-floor(.0163*NE[contELM3][1]+1500);
						emax3=2*NE[contELM3][1]+1500;

					}

					if(emin3<1){
							
						emin3=1;

					}

					if(emax3>=elm_tot){
							
						emax3=elm_tot;

					}

					contELM3 = false;	//resent containing element

//****Element Search****

					for(int e=emin3; e<=emax3; e++){

						for(int i=1; i<=3; i++){ 

							xloc[i] = lon[NE[e][i]];
							yloc[i] = lat[NE[e][i]];
							delx = xloc[i] - XP3[m];
							dely = yloc[i] - YP3[m];
							D = pow((pow(delx,2.0)+pow(dely,2.0)),0.5);
							theta = atan2(dely,delx);    			
							vx[i] = D*cos(theta);
							vy[i] = D*sin(theta);

						}

						crossprod1 = vx[1]*vy[2] - vy[1]*vx[2];
						crossprod2 = vx[2]*vy[3] - vy[2]*vx[3];
						crossprod3 = vx[3]*vy[1] - vy[3]*vx[1];

						if(crossprod1>=0 && crossprod2>=0 && crossprod3>=0){
			
							contELM3 = ELM[e];

//****Barycentric elevation interpolation****

				w1 =((yloc[2]-yloc[3])*(XP3[m]-xloc[3])+(xloc[3]-xloc[2])*(YP3[m]-yloc[3]))/
				    ((yloc[2]-yloc[3])*(xloc[1]-xloc[3])+(xloc[3]-xloc[2])*(yloc[1]-yloc[3]));

				w2 =((yloc[3]-yloc[1])*(XP3[m]-xloc[3])+(xloc[1]-xloc[3])*(YP3[m]-yloc[3]))/
				    ((yloc[2]-yloc[3])*(xloc[1]-xloc[3])+(xloc[3]-xloc[2])*(yloc[1]-yloc[3]));

				w3 = 1-w1-w2;

				Pelv3[m] = (w1*elv[NE[e][1]]+w2*elv[NE[e][2]]+w3*elv[NE[e][3]])/(w1+w2+w3);

		
							if(Pelv3[m] > 0){
										
								Phold3[m] = Pelv3[m];	//The last wet point 
										
 
							}

						break;	//out of for elements loop

						}	//close if(crossprod > 0)

					}	//close restricted element search for loop


				if(contELM3 == false){
//cout<<"Long Search Q3************** steps: "<<step3<<endl;

					for(int e=1; e<elm_tot; e++){

						for(int i=1; i<=3; i++){ 

							xloc[i] = lon[NE[e][i]];
							yloc[i] = lat[NE[e][i]];
							delx = xloc[i] - XP3[m];
							dely = yloc[i] - YP3[m];
							D = pow((pow(delx,2.0)+pow(dely,2.0)),0.5);
							theta = atan2(dely,delx);    			
							vx[i] = D*cos(theta);
							vy[i] = D*sin(theta);

						}

						crossprod1 = vx[1]*vy[2] - vy[1]*vx[2];
						crossprod2 = vx[2]*vy[3] - vy[2]*vx[3];
						crossprod3 = vx[3]*vy[1] - vy[3]*vx[1];

						if(crossprod1>=0 && crossprod2>=0 && crossprod3>=0){
			
							contELM3 = ELM[e];

//****Barycentric elevation interpolation****

				w1 =((yloc[2]-yloc[3])*(XP3[m]-xloc[3])+(xloc[3]-xloc[2])*(YP3[m]-yloc[3]))/
				    ((yloc[2]-yloc[3])*(xloc[1]-xloc[3])+(xloc[3]-xloc[2])*(yloc[1]-yloc[3]));

				w2 =((yloc[3]-yloc[1])*(XP3[m]-xloc[3])+(xloc[1]-xloc[3])*(YP3[m]-yloc[3]))/
				    ((yloc[2]-yloc[3])*(xloc[1]-xloc[3])+(xloc[3]-xloc[2])*(yloc[1]-yloc[3]));

				w3 = 1-w1-w2;

				Pelv3[m] = (w1*elv[NE[e][1]]+w2*elv[NE[e][2]]+w3*elv[NE[e][3]])/(w1+w2+w3);

		
							if(Pelv3[m] > 0){
										
								Phold3[m] = Pelv3[m];	//The last wet point 
										
 
							}

						break;	//out of for elements loop

						}	//close if(crossprod > 0)

					}	//close restricted element search for loop

//cout<<"!!!!contELM3: "<<contELM3<<"   emin:  "<<emin3<<"   emax: "<<emax3<<endl;

				}	//close if contELM3 == false;


//if still no containing element, out of bounds

					if(contELM3 == false){

						Pelv3[m]=0;	
						//break;			//break to...

					}

				}	//close while Pelv (XP,YP) > 0 loop

				outfile7<<dephold3.size()<<endl;

				for(int h=0; h<dephold3.size(); h++){
				
					outfile7<<dephold3[h]<<endl;

				}

				dpast3[m] = Phold3[m]/((Phold3[m] + abs(Pelv3[m]))/res);

				dlon3 = abs(XP3[m] - lon0);			
				dlat3 = abs(YP3[m] - lat0);


				dhav3 = 2*Rearth*asin(sqrt(sin(((dlat3)/2)*pi/180)*
							  sin(((dlat3)/2)*pi/180)+
						   (cos(lat0*pi/180)*cos(YP3[m]*pi/180)*
							  sin(((dlon3)/2)*pi/180)*
							  sin(((dlon3)/2)*pi/180))))-res+dpast3[m];

				if(step3 == 1){	
							
					//dold3[m] = dpast3[m];

					dhav3=dpast3[m];
			
					if(contELM3 == false){

						//dold3[m] = 0;

						dhav3=0;

					}

				}

				if(step3!=1 && contELM3 ==false){
		

					dhav3 = 2*Rearth*asin(sqrt(sin(((dlat3)/2)*pi/180)*
								  sin(((dlat3)/2)*pi/180)+
							   (cos(lat0*pi/180)*cos(YP3[m]*pi/180)*
								  sin(((dlon3)/2)*pi/180)*
								  sin(((dlon3)/2)*pi/180))))-(res/2);

				}

				outfile6<<dhav3<<" ";

			}	//close for m


//*************************************Q2********************************************************************


			for(int m=0; m<Ares; m++){		//including 0 excluding Ares (0-80 deg)	

				outfile7<<j<<" ";
				outfile7<<m+Angres+Angres+Angres-2<<" ";
				dephold2.clear();

				step2 = 0;
				Pelv2[m] = elv[j];
				Phold2[m] = elv[j];
				dpast2[m]=0;

				emin2=2*j-floor(.0163*j+1500);
				emax2=2*j+1500;

				while(Pelv2[m] > 0){			//while (XP,YP) is wet		

					dephold2.push_back(Pelv2[m]);	//stores every Pelv value

					step2=step2+1;

					XP2[m] = lon0 - (dx*cos((pi/180)*(Angres*m)))*step2;
					YP2[m] = lat0 + (dy*sin((pi/180)*(Angres*m)))*step2;

//if first step already taken, compare the 1st neighboring node of the previous containing element 
//so that the restricted element search is wrt the containing element instead of the start node j


//Pelv1 at step1==2 corresponds to the first step's interpolated depth. Calculate slope from start to here
//BUT slope values are stored in slopehold array in flipped indexes. Q1<-->Q3 and Q2<-->Q4 because we want
//the index to correspond to upwind slope for determining wave breaking

					if(step2==2){

						slopehold[j][m+Ares2+Ares2+Ares2+1]=(startelv-Pelv2[m])/res;

//if slope is less than 1/1000 set to 0, this inludes negative slopes (from shallow to deeper water)

						if(slopehold[j][m+Ares2+Ares2+Ares2+1] < 1/1000){ 

							slopehold[j][m+Ares2+Ares2+Ares2+1]=0;

						}

						outfile8<<slopehold[j][m+Ares2+Ares2+Ares2+1]<<" ";

					}


					if(step2>1){

						emin2=2*NE[contELM2][1]-floor(.0163*NE[contELM2][1]+1500);
						emax2=2*NE[contELM2][1]+1500;

					}

					if(emin2<1){
							
						emin2=1;

					}

					if(emax2>=elm_tot){
							
						emax2=elm_tot;

					}

					contELM2 = false;	//resent containing element


//****Element Search****
					for(int e=emin2; e<=emax2; e++){

						for(int i=1; i<=3; i++){ 

							xloc[i] = lon[NE[e][i]];
							yloc[i] = lat[NE[e][i]];
							delx = xloc[i] - XP2[m];
							dely = yloc[i] - YP2[m];
							D = pow((pow(delx,2.0)+pow(dely,2.0)),0.5);
							theta = atan2(dely,delx);    			
							vx[i] = D*cos(theta);
							vy[i] = D*sin(theta);

						}

						crossprod1 = vx[1]*vy[2] - vy[1]*vx[2];
						crossprod2 = vx[2]*vy[3] - vy[2]*vx[3];
						crossprod3 = vx[3]*vy[1] - vy[3]*vx[1];

						if(crossprod1>=0 && crossprod2>=0 && crossprod3>=0){
			
							contELM2 = ELM[e];


//****Barycentric elevation interpolation****

				w1 =((yloc[2]-yloc[3])*(XP2[m]-xloc[3])+(xloc[3]-xloc[2])*(YP2[m]-yloc[3]))/
				    ((yloc[2]-yloc[3])*(xloc[1]-xloc[3])+(xloc[3]-xloc[2])*(yloc[1]-yloc[3]));

				w2 =((yloc[3]-yloc[1])*(XP2[m]-xloc[3])+(xloc[1]-xloc[3])*(YP2[m]-yloc[3]))/
				    ((yloc[2]-yloc[3])*(xloc[1]-xloc[3])+(xloc[3]-xloc[2])*(yloc[1]-yloc[3]));

				w3 = 1-w1-w2;

				Pelv2[m] = (w1*elv[NE[e][1]]+w2*elv[NE[e][2]]+w3*elv[NE[e][3]])/(w1+w2+w3);

							if(Pelv2[m] > 0){
										
								Phold2[m] = Pelv2[m];	//The last wet point 
		 										
							}

						break;	//out of for elements loop

						}	//close if(crossprod > 0)

					}	//close restricted element search for loop

				if(contELM2 == false){

//cout<<"Long Search Q2************** steps: "<<step2<<endl;

					for(int e=1; e<elm_tot; e++){

						for(int i=1; i<=3; i++){ 

							xloc[i] = lon[NE[e][i]];
							yloc[i] = lat[NE[e][i]];
							delx = xloc[i] - XP2[m];
							dely = yloc[i] - YP2[m];
							D = pow((pow(delx,2.0)+pow(dely,2.0)),0.5);
							theta = atan2(dely,delx);    			
							vx[i] = D*cos(theta);
							vy[i] = D*sin(theta);

						}

						crossprod1 = vx[1]*vy[2] - vy[1]*vx[2];
						crossprod2 = vx[2]*vy[3] - vy[2]*vx[3];
						crossprod3 = vx[3]*vy[1] - vy[3]*vx[1];

						if(crossprod1>=0 && crossprod2>=0 && crossprod3>=0){
			
							contELM2 = ELM[e];


//****Barycentric elevation interpolation****

				w1 =((yloc[2]-yloc[3])*(XP2[m]-xloc[3])+(xloc[3]-xloc[2])*(YP2[m]-yloc[3]))/
				    ((yloc[2]-yloc[3])*(xloc[1]-xloc[3])+(xloc[3]-xloc[2])*(yloc[1]-yloc[3]));

				w2 =((yloc[3]-yloc[1])*(XP2[m]-xloc[3])+(xloc[1]-xloc[3])*(YP2[m]-yloc[3]))/
				    ((yloc[2]-yloc[3])*(xloc[1]-xloc[3])+(xloc[3]-xloc[2])*(yloc[1]-yloc[3]));

				w3 = 1-w1-w2;

				Pelv2[m] = (w1*elv[NE[e][1]]+w2*elv[NE[e][2]]+w3*elv[NE[e][3]])/(w1+w2+w3);

							if(Pelv2[m] > 0){
										
								Phold2[m] = Pelv2[m];	//The last wet point 
		 										
							}

						break;	//out of for elements loop

						}	//close if(crossprod > 0)

					}	//close all element search for loop

//cout<<"!!!!contELM2: "<<contELM2<<"   emin:  "<<emin2<<"   emax: "<<emax2<<endl;

				}	//close if contELM2 == false

//if containing element still false, out of bounds

					if(contELM2 == false){

						Pelv2[m]=0;	

					}
		
				}	//close while Pelv (XP,YP) > 0 loop

				outfile7<<dephold2.size()<<endl;

				for(int h=0; h<dephold2.size(); h++){
				
					outfile7<<dephold2[h]<<endl;

				}

				dpast2[m] = Phold2[m]/((Phold2[m] + abs(Pelv2[m]))/res);

				dlon2 = abs(XP2[m] - lon0);			
				dlat2 = abs(YP2[m] - lat0);


				dhav2 = 2*Rearth*asin(sqrt(sin(((dlat2)/2)*pi/180)*
							  sin(((dlat2)/2)*pi/180)+
						   (cos(lat0*pi/180)*cos(YP2[m]*pi/180)*
							  sin(((dlon2)/2)*pi/180)*
							  sin(((dlon2)/2)*pi/180))))-res+dpast2[m];

				if(step2 == 1){	

					dhav2=dpast2[m];
			
					if(contELM2 == false){
						
						dhav2=0;

					}

				}

				if(step2!=1 && contELM2 ==false){
		

					dhav2 = 2*Rearth*asin(sqrt(sin(((dlat2)/2)*pi/180)*
								  sin(((dlat2)/2)*pi/180)+
							   (cos(lat0*pi/180)*cos(YP2[m]*pi/180)*
								  sin(((dlon2)/2)*pi/180)*
								  sin(((dlon2)/2)*pi/180))))-(res/2);

				}

				outfile6<<dhav2<<" ";

			}	//close for m

		}	//close if wet 

		outfile6<<endl;
		outfile8<<endl;

	}		//close for all nodes loop


outfile6.close();
outfile7.close();
outfile8.close();

	time_req=clock()- time_req;

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




//**************Begin Averaging Schemes***********************************************************************


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

	float time_req2;
	time_req2=clock();			//initializing the clock for timing tasks
	c=0;

	if(infile1.is_open()){

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

			line=line+1;
			
			if(line>3 && line<node_tot+3){	//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

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

						}


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

							fetch[node1][x]=atof(temp);		
							
						}

					}	//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;

	}

	time_req2=clock()- time_req2;

	cout<<"***Time taken: "<<(float)time_req2/CLOCKS_PER_SEC<<" seconds***"<<endl;

//****************Reading in Depth Data and computing IDW depth

	line=0;

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


	//outfile9<<"p-values determined from variable calculation with U="<<Up<<endl;
	outfile9<<"IDW depth averages with p="<<p<<endl;

	float time_req3;
	time_req3=clock();
	c=0;

	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

//format is:    node   dir.index   #steps
//		     depth at step1
//		     depth at step2
//			   ...
//		     depth at stepN (node before land)
//	repeat for next dir. index, then once all directions are performed, repeat for next node
		   

			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[node1][direction]=depth_tot/denominator;
						


						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
				
					node1=atof(content);

					outfile9<<endl;
					outfile9<<node1<<"  ";
					//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
//THIS SHOULD NOT BE USED UNTIL A BETTER P-VALUE RELATIONSHIP IS DETERMINED
//****************************************************************************************

					p=p_temp;

					if(p_temp==0){

				p=log((10924*exp(-0.064*Up))/depth_temp) / log(fetch[node1][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[node1][direction]);

						}

						if(p_temp!=0){

							p=p_temp;
						}

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

					}

					outfile9<<p<<" ";

					depth_tot = (2*depth_temp)/(pow(res,p));
					denominator = 2/(pow(res,p));
					//cout<<"depth["<<x2<<"]: "<<depth_temp<<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<<endl;

				}


			}	//close if line > 2

		}	//close while getline

	}	//close if infile2 is_open

	else{

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

	}

outfile9.close();


	time_req3=clock()- time_req3;

	cout<<"***Time taken: "<<(float)time_req3/CLOCKS_PER_SEC<<" seconds***"<<endl;


//storing IDW depth avgs in an output file

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

		outfile2<<i<<"  ";

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

			outfile2<<d[i][j]<<" ";

		}

		outfile2<<endl;

	} 

	outfile2.close();



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



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

	outfile1<<node_tot<<"  "<<Angres<<endl;

	float time_req5;
	time_req5=clock();


	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;

			for(int n=-floor(Eff_res/Angres); n<=floor(Eff_res/Angres); n++){

//rounding down to the nearest integer division for direction given the Effective fetch resolution specified

				//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 
					//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

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

	//for values that go above max direction

					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 +/- Eff_res deg on each side of start node				

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

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

					Eff_top = fetch[j][k]*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


	outfile1.close();



//*****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(slopehold[j][i] < 1){
			
					Mtop=Mtop+slopehold[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;



	time_req5=clock()- time_req5;

	cout<<"***Time taken: "<<(float)time_req5/CLOCKS_PER_SEC<<" seconds***"<<endl;



return 0;
}