% apply t_tide to Rosario Sea Spider data
% source code from https://www-old.eoas.ubc.ca/~rich/#T_Tide
%
% J. Thomson, Mar 2025

clear all

load('SeaSpider_Sig250_Rosario.mat')

%% remove points out of water (instead of just making NaN in other processing)

mindepth = 80;

t( waterdepth < mindepth ) = [];
u( waterdepth < mindepth, : ) = [];
v( waterdepth < mindepth, : ) = [];
waterdepth( waterdepth < mindepth ) = [];

%% tidal evelataion

%h = detrend(waterdepth);
%[hHARMONICS,hPREDICT]=t_tide(h,'interval',1/6,'lattitude',48); %percent var predicted/var original= 96.0 %

%% hubheight velocity (2D)

hubbin = 70; % MLLW - R - offset : 90 - 13.5 - 3.  
uv = complex(u(:,hubbin), v(:,hubbin));
%[uvHARMONICS,uvPREDICT]=t_tide(uv,'interval',1/6,'lattitude',48,'start',t(1),'shallow','M10','error','linear','synthesis',1); %percent var predicted/var original= 96.2 %
[uvHARMONICS,uvPREDICT]=t_tide(uv,'interval',1/6,'lattitude',48,'start',t(1)); %percent var predicted/var original= 96.2 %


%% predict full year 

dt = 1; % spacing (hours)

t_annual = t(1)+ [1:dt:(365*24)]'; % one hour spacing
uvPREDICT_annual = t_predic(t_annual, uvHARMONICS, 'latitude',48);

save SeaSpider_Sig250_Rosario.mat uvHARMONICS uvPREDICT_annual -append 

%% time series
figure(1), clf
plot(t,real(uv),'c',t,imag(uv),'b'), hold on
plot(t,real(uvPREDICT),'m',t,imag(uvPREDICT),'r')
datetick

%% u,v bin scatter

figure(2), clf
subplot(1,2,1)
binscatter(real(uvPREDICT_annual), imag(uvPREDICT_annual),[30 30])
set(gca,'fontsize',14,'fontweight','demi')
axis([-3 3 -3 3]), xlabel('u [m/s]'), ylabel('v [m/s]'), grid
title('Annual prediction from harmonics')
cb1 = colorbar; cb1.Ticks = linspace(0,593,6); cb1.TickLabels = num2str([0:0.2:1]');

subplot(1,2,2)
binscatter(real(uv), imag(uv), [30 30])
set(gca,'fontsize',14,'fontweight','demi')
axis([-3 3 -3 3]), xlabel('u [m/s]'), ylabel('v [m/s]'), grid
title('103 day observations')
cb2 = colorbar; cb2.Ticks = linspace(0,700,6); cb2.TickLabels = num2str([0:0.2:1]');

print -dpng Rosario_SeaSpider_Sig250_harmonicUVscatter.png


%% speed exceedance

figure(3), clf
eX = eprob(abs(uv));
plot(eX.sort,eX.eprob * 100,'b-','LineWidth',2), hold on
eX = eprob(abs(uvPREDICT));
plot(eX.sort,eX.eprob * 100,'g-','LineWidth',2), hold on
eX = eprob(abs(uvPREDICT_annual));
plot(eX.sort,eX.eprob * 100,'m-','LineWidth',2), hold on
legend('103 day observations','103 day reconstruction','1 yr reconstruction')
ylabel('Exceedance Probability (%)')
xlabel('Scalar speed at hub depth [m/s]')
set(gca,'fontsize',14,'fontweight','demi')

print -dpng Rosario_SeaSpider_Sig250_exceed.png


%% AEP again

hubspeed = abs(uvPREDICT_annual);

% turbine params
cutin = 0.7; % m/s
rated = 2.25; % m/s
efficiency = 0.39;
R = 13.5; % rotor radius (m)
hubdepth = 3.5 + R; % hubheight (m) below surface

% power regions
region1 = find( hubspeed < cutin );
region2 = find( hubspeed > cutin & hubspeed < rated );
region3 = find( hubspeed > rated );

power(region1) =  0;
power(region2) = 0.5 * 1030 * 3.14 * R^2 * hubspeed(region2).^3 * efficiency;
power(region3) = 0.5 * 1030 * 3.14 * R^2 * rated.^3 * efficiency;
AEP = trapz( power ) * dt;
AEP_GWh = AEP./1e9;  % 2.2 GWh, within 10% of model
AEP = 2*AEP % two rotors