//Code for creating fort.22 wind file in NWS=5 format
//Convention is Nautical for PARAM, SWAN, and ADCIRC (0 deg=North, 90 deg=East, 180 deg=South, 270 deg=West)
//PARAM and ADCIRC use Meteorlogical convention (wind from North is WVY= -value)
//SWAN uses Oceanographis convention (wind from North is WVY = +value)
//
//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>

using namespace std;

int main(){

int IWTIME;		//time of the wind field
int IREFYR = 19;	//start time (year,month,day,hr,min)
int IREFMO = 11;
int IREFDAY = 13;
int IREFHR = 00;
int IREFMIN = 00;
int NDSETSE;
int NSPOOLGE = 3600;		//the number of time steps at which information is written to fort.63
int dims = 2;			//dimension of the wind dataset
double pi = 4*atan(1);
double umax;
double dir,u;
int theta_res;
double wind_converge=0.0001;
int hr_count = 0;
int day_count = 0;
double xwind,ywind,WVX,WVY,PRN=0;
double WVXa[NDSETSE],WVYa[NDSETSE],WVXb[NDSETSE],WVYb[NDSETSE],WVXAD[NDSETSE],WVYAD[NDSETSE];
int index=0;
int wind_type;
bool check;


cout<<endl;
cout<<"Enter wind field desired (0=const. wind, 1=theta wind, 2=sin wind): ";
cin>>wind_type;
cout<<endl;

cout<<"Enter total number of datasets to be written to fort.63 (ADCIRC NDSETSE value): ";
cin>>NDSETSE;

	if(wind_type==1){

		check=false;

		if(360%NDSETSE == 0){

			check=true;
			cout<<endl;

		}

		while(check==false){

			cout<<"Inter valid NDSETSE (evenly divisible by 360): ";
			cin>>NDSETSE;

			if(360%NDSETSE == 0){

				check=true;

			}

		}

	}


//double WVXAD[NDSETSE];
//double WVYAD[NDSETSE];

cout<<"Enter max wind speed (m/s): ";
cin>>umax;
cout<<endl;
						
ofstream outfile1("PARAM_theta_U_30.22");
ofstream outfile3("SWAN_theta_U_30.22"); 
ofstream outfile2("ADCIRC_theta_U_30.22");



//*****************************theta wind********************************

if(wind_type==1){

ofstream outfile1("PARAM_theta_U_30.22");
ofstream outfile2("SWAN_theta_U_30.22");
ofstream outfile3("ADCIRC_theta_U_30.22"); 

	int theta_res;

	cout<<"Enter angular resolution for wind field rotation (deg): ";
	cin>>theta_res;
	cout<<endl;

//********iterating to solve for WVX and WVY given max wind speed and angular resolution

	for(int theta=0; theta<=360-theta_res; theta=theta+theta_res){

		index=index+1;

		for(double WY=wind_converge; WY<=umax; WY=WY+wind_converge){

			double theta_check;

			WVX=sqrt(umax*umax - WY*WY);
			theta_check=atan(WVX/WY)*180/pi;

			if(abs(theta_check-theta) < wind_converge*30){
			
				WVY=WY;

				WVXa[index]=WVX;
				WVYa[index]=WVY;

cout<<"converged for theta: "<<theta<<" WVXa["<<index<<"]="<<WVXa[index]<<" WVYa["<<index<<"]="<<WVYa[index]<<endl;
				break;

			}

		}

	}

//Rearranging the indexes for all 4 quadrents

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

		int ii;

		if(i>=1 && i<=10){		//Q1 (0-90 deg)

			WVXb[i]=WVXa[i];
			WVYb[i]=WVYa[i];
			WVXAD[i]=-WVXa[i];
			WVYAD[i]=-WVYa[i];

		}

		if(i>10 && i<=20){		//Q4 (90-180 deg)

			ii=(20-i);

			WVXb[i]=WVXa[ii];
			WVYb[i]=-WVYa[ii];
			WVXAD[i]=-WVXa[ii];
			WVYAD[i]=WVYa[ii];
		}

		if(i>18 && i<=28){		//Q3 (180-270 deg)

			ii=i-18;

			WVXb[i]=-WVXa[ii];
			WVYb[i]=-WVYa[ii];
			WVXAD[i]=WVXa[ii];
			WVYAD[i]=WVYa[ii];
		}

		if(i>28 && i<37){		//Q2 (270-360 deg)

			ii=(38-i);

			WVXb[i]=-WVXa[ii];
			WVYb[i]=WVYa[ii];
			WVXAD[i]=WVXa[ii];
			WVYAD[i]=-WVYa[ii];
		}
		//cout<<"WVXb["<<i<<"] = "<<WVXb[i]<<"  WVYb["<<i<<"] = "<<WVYb[i]<<endl;
		//cout<<"WVXAD["<<i<<"] = "<<WVXAD[i]<<"  WVYAD["<<i<<"] = "<<WVYAD[i]<<endl;
		//cout<<endl;
	}

	cout<<"Writing fort.22"<<endl;

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

		cout<<"timestep "<<i<<" wind dir "<<(i-1)*theta_res<<"(deg N)"<<endl;

		for(int n=1; n<=126772; n++){//can split the domain with multiple loops here for different wind

			//umax=30;	//function for varying u with time: sin(pi*i/NDSETSE)*umax;
			PRN=10.3;	//surface atmospheric pressure in meters of water(can vary with time)

				//PARAM

				outfile1<<n<<" "<<WVXAD[i]<<" "<<WVYAD[i]<<" "<<" "<<PRN<<endl;

				//ADCIRC

				outfile3<<n<<" "<<WVXAD[i]<<" "<<WVYAD[i]<<" "<<" "<<PRN<<endl;

		
			if(n==1){

				cout<<"PARAM    WVX: "<<WVXAD[i]<<"  WVY: "<<WVYAD[i]<<endl;
				cout<<"ADCIRC   WVX: "<<WVXAD[i]<<"  WVY: "<<WVYAD[i]<<endl;

			}

		}	// close for all nodes
	
	}	// close for NDSETSE



	for(int n=1; n<=NDSETSE; n++){		//creating file for swan wind field input 

		cout<<"SWAN     WVX: "<< WVXb[n]<<"  WVY: "<<WVYb[n]<<endl;

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

			outfile2<<WVXb[n]<<endl;

		}

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

			outfile2<<WVYb[n]<<endl;

		}

	}	//close swan wind field


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

}	//close theta wind option

//*************************Constant Wind Option***************************

if(wind_type==0){

ofstream outfile1("PARAM_const_U_30.22");
ofstream outfile2("SWAN_const_U_30.22");
ofstream outfile3("ADCIRC_const_U_30.22"); 

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

		cout<<"writing fort.22 (constant wind) at timestep "<<i<<endl;

		for(int n=1; n<=126772; n++){//can split the domain with multiple loops here for different wind

			PRN=10.3;	//surface atmospheric pressure in meters of water(can vary with time)
			xwind=0;
			ywind=umax;


			if(n==1){
				cout<<"PARAM    U="<<-umax<<"(m/s)"<<endl;
			}
				//PARAM
				outfile1<<n<<" "<<xwind<<" "<<-ywind<<" "<<PRN<<endl;

			if(n==1){
				cout<<"ADCIRC   U="<<-umax<<"(m/s)"<<endl;
			}
				//ADCIRC
				outfile3<<n<<" "<<xwind<<" "<<-ywind<<" "<<PRN<<endl;
		
		}
	
	}	// close for NDSETSE


	for(int n=1; n<=NDSETSE; n++){		//creating file for swan wind field input 

		cout<<"SWAN     U="<<umax<<"(m/s)   at timestep "<<n<<endl;

		xwind=0;
		ywind=umax;


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

			outfile2<<xwind<<endl;

		}

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

			outfile2<<ywind<<endl;

		}

	}


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

}	//close wind type (constant)



//*************************Sin Ramping Wind Option***************************


if(wind_type==2){

u=0;

ofstream outfile1("PARAM_sin_U_30.22");
ofstream outfile2("SWAN_sin_U_30.22");
ofstream outfile3("ADCIRC_sin_U_30.22"); 

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

		//can split the domain with multiple loops here for different wind

		cout<<"writing fort.22  (sin wind) at timestep "<<i<<endl;		

		for(int n=1; n<=126772; n++){

			u=sin(pi*i/NDSETSE)*umax;

			PRN=10.3;	//surface atmospheric pressure in meters of water(can vary with time)


			if(n==1){
				cout<<"PARAM     fort.22    U="<<-u<<"(m/s)"<<endl;
			}

			outfile1<<n<<" "<<0<<" "<<-u<<" "<<" "<<PRN<<endl;

			if(n==1){
				cout<<"ADCIRC    fort.22    U="<<-u<<"(m/s)"<<endl;
			}

			outfile3<<n<<" "<<0<<" "<<-u<<" "<<" "<<PRN<<endl;
				
		}
	
	}	// close for NDSETSE


	u=0;

	//SWAN 

	for(int n=1; n<=NDSETSE; n++){		

		u=sin(pi*n/NDSETSE)*umax;

		cout<<"SWAN      fort.22    U="<<u<<"(m/s)"<<endl;


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

			outfile2<<0<<endl;

		}

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

			outfile2<<u<<endl;

		}

	}

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

}	//close sin wind option


return 0;
}