%% SWAN loop procedure
clear; 
clc;
close all;
tic;

%% loading the data
data = readtable('Offshore Data.txt');  % wave climate data
bm = readtable('Bathymetry.txt');       % bathymetry

Hs  =   data{:,1};  % wave height [m]
Tp  =   data{:,2};  % wave period [s]
Dir =   data{:,3};  % wave direction [deg], relative to the north
w_spd = data{:,4};  % wind speed [m/s]
w_dir = data{:,5};  % wind direction [deg], relative to the north

xs  =   bm{:,1};    % distance [m]
d   =   bm{:,2};    % depth [m]

%% calculation settings
coast_angle = 105;  % angle of the coastline [deg], clockwise to the north

% grid
ng = 100;           % amount of grid points desired
xend =  xs(end);    % end of the grid
dx = (xend)/ng;     % gridsize

% savepath for output data
savepath = 'C:\Users\Menno\OneDrive\Documenten\TU Delft Civiele Techniek\Thesis\Coding\';

%% generate project_001_bot.DAT
for i = 1:(ng+1)
    s(i) = dx*(i-1);
    botdat(i) = interp1(xs,-d,(i-1)*dx);
end
writematrix(botdat', 'project_001_bot.DAT')

%% direction check
if min(Dir) < coast_angle + 90 - 70 || max(Dir) > coast_angle +90 + 70
    disp('Warning! Wave direction out of bounds! Calculation may not be accurate']
    return
end

%% SWAN loop
for ii=1:length(Hs)
try
% a new project file is written every iteration
% wave and wind variables are changed every iteration
prjfile = fopen('project_001.SWN', 'w');
A = [       "PROJECT 'project_001''A577' '1'",  
            "MODE STATIONARY ONED",        
            "SET DEPMIN=0.01 MAXMES=999 MAXERR=3 PWTAIL=4",              
            "SET LEVEL 0",             
            "SET CARTESIAN",             
            "SETUP 0",
            strcat('CGRID REGULAR 0 0 0',{' '},num2str(xend),' 0.',{' '}, num2str(ng), ' 0 &');                 
            "CIRCLE 36 0.03 0.8 30",       
            "BOUND SHAPE JONSWAP 3.3 PEAK DSPR POWER",
            strcat('INPGRID BOTTOM REGULAR 0 0 0',{' '},num2str(ng),' 0',{' '}, num2str(dx), ' 1');          
            "READ BOTTOM 1. 'project_001_bot.dat '1 0 FREE",
            "!NPGRID CURRENT REGULAR 0 0 0 100 0 10 1",            
            "!EAD CURRENT",                                                          
            strcat('BOUN SIDE W CONSTANT PAR',{' '},num2str(Hs(ii)),{' '},num2str(Tp(ii)),{' '},num2str(wrapTo360(Dir(ii)-90-coast_angle)),{' '},'2');    
            "",
            strcat('WIND',{' '},num2str(round(w_spd(ii,1),1)),{' '},num2str(round(wrapTo360(w_dir(ii,1)-90-coast_angle),1)));
            "GEN3 KOMEN",  
            "BREAK 1. 0.73",            
            "FRIC JONSWAP",           
            "TRIAD",             
            "NUM ACCUR 0.02 0.02 0.02 98 15",          
            "OUTPUT OPTIONS '%' TABLE 16 BLOCK 9 1000 SPEC 8",
            strcat("CURVE 'curve' 0 0",{' '},num2str(ng),{' '}, num2str(xend), ' 0');                
            "Table 'curve' HEADER 'project_001.tab' XP YP HSIGN RTP TM01 TMM10 DIR &",  
            "DSPR DEPTH SETUP",       
            "POIN 'S0'1 .00", 
            "SPEC 'S0' SPEC1D ABS 'project_001_point1.sp1'",
            "SPEC 'S0' SPEC2D ABS 'project_001_point1.sp2'",
            "POIN 'S1'56709 .00", 
            "SPEC 'S1' SPEC1D ABS 'project_001_point2.sp1'",
            "SPEC 'S1' SPEC2D ABS 'project_001_point2.sp2'",
            "TEST 0 0",         
            "COMPUTE",
            "STOP",
            ""];
fprintf(prjfile,'%s\n', A);
fclose(prjfile);

% perform the SWAN calculation
[status,cmdout] = system(['swanrun ' 'project_001']);
disp(['Succes! Calculation completed for iteration ', num2str(ii)])

% get results
result      =   txt2mat('project_001.tab');
HsEns(ii,:) =   result(:,3);
TpEns(ii,:) =   result(:,4);
DirEns(ii,:)=   result(:,7);

% error catcher, in case the process fails display the iteration number
catch ME
    disp(['Something went wrong! Error occured at iteration ', num2str(ii)])
end
end

%% new data
Hs_on       =   nonzeros(HsEns(:,end));
Tp_on       =   nonzeros(TpEns(:,end));
Dir_on      =   wrapTo360(nonzeros(DirEns(:,end)+90+coast_angle));


%% save data
% Create a new table with the splines
newtable = table(Hs_on, Tp_on, Dir_on, w_spd, w_dir, 'VariableNames', {'Wave height','Wave period','Wave direction','Wind speed','Wind direction'} );

% Write data to text file
filename = [savepath, 'Onshore Data.txt'];
writetable(newtable, filename)

load handel
sound(y,Fs)
toc;