% Analysis of Rosario Strait ADCP data collected Oct 2024-Jan 2025
% using surface buoy with downlooking Aquadopp
% Velocity and power statistics for OPALCO tidal project
% J. Thomson, Feb 2025

clear all, close all

fpath = '~/Desktop/RosarioData/Aquadopp_Rosario_2025/'

prefix  =  'OPALCO02';

filebase = [ fpath prefix ];

sen = load( [ filebase '.sen' ]);

% timestamps
month = sen(:, 1);
day = sen(:, 2);
year = sen(:, 3);
hour = sen(:, 4);
minute = sen(:, 5);
second = sen(:, 6);
t = datenum( year, month, day, hour, minute, second);

% ancillary data
voltage = sen( :, 9);
soundspd = sen( :, 10); % m/s
heading = sen( :, 11); % deg M
pitch = sen( :, 12); % deg
roll = sen( :, 13); % deg
pres = sen( :, 14); % dbar *** GAGE VALUE, NEEDS ATMOSPHERIC CORRECTION ***
temp = sen( :, 15); % deg C

% doppler data as arrays of time x cell: velocity (m/s), amplitude (cts), correlation (%)
u = load( [ filebase '.v1' ]);  
v = load( [ filebase '.v2' ]);  
w = load( [ filebase '.v3' ]);  
a1 = load( [ filebase '.a1' ]);  
a2 = load( [ filebase '.a2' ]);  
a3 = load( [ filebase '.a3' ]);  

% initialize 
coor = 'ENU';


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

%% cell depth (d) below surface

[pts cells] = size(u);
res = 1; % m (see .hdr file)
blanking = 0.4; % m (see .hdr file)
d = blanking + [1:cells]*res;  % depth of cell below surface (see .hdr file)

waterdepth = 86; % m (survey depth)

%% remove data when out of water

minamp = 25;
ondeck = find(a1<minamp | a1<minamp | a3<minamp);
u(ondeck) = NaN;
v(ondeck) = NaN;
w(ondeck) = NaN;


%% look at quality metrics / ancillary data

% orientation
figure(1), clf
subplot(3,1,1), plot(t, pitch, '.', t, roll, '.'), datetick, legend('pitch','roll'), %set(gca,'YLim',[-10 10])
subplot(3,1,2), plot(t, heading, '.'), datetick, legend('heading'), set(gca,'YLim',[0 360])
subplot(3,1,3), plot(t, temp, '.'), datetick, legend('Temperature'), set(gca,'YLim',[6 14])
print -dpng SurfaceBuoy_AQD_Rosario_PRH.png

% amplitude
figure(2), clf
for i=1:3
    subplot(3,1,i), pcolor(t, d, eval(['a' num2str(i) ])'), shading flat, datetick, hold on, axis tight
    set(gca,'YDir','reverse'), ylabel('depth [m]')
    cb = colorbar; cb.Label.String = 'dB';
end
print -dpng SurfaceBuoy_AQD_Rosario_Amp.png


%% optional transformation to beam coordinates (see https://github.com/NortekSupport/coordinatetransforms/blob/main/Transform.m)
T = [1.5774 -0.7891 -0.7891; 0.0000 -1.3662 1.3662; 0.3677 0.3677 0.3677];  % see .hdr file
T_original = T; 
T(2,:) = -T(2,:);   T(3,:) = -T(3,:);  % downlooking T has negative of rows 2 and 3

for ti=1:length(t)
    hh = pi*(heading(ti)-90)/180; pp = pi*pitch(ti)/180; rr = pi*roll(ti)/180; % angles in radians
    H = [cos(hh) sin(hh) 0; -sin(hh) cos(hh) 0; 0 0 1] ; % heading matrix
    P = [cos(pp) -sin(pp)*sin(rr) -cos(rr)*sin(pp);  0  cos(rr)  -sin(rr);  sin(pp) sin(rr)*cos(pp)  cos(pp)*cos(rr)]; % tilt matrix
    HPT = H*P*T; % full transformation matrix
    for di=1:length(d)
        enu = [u(ti,di); v(ti,di); w(ti,di);];
        beam = inv(HPT)*enu;
        xyz = T_original*beam;
        u(ti,di) = beam(1); v(ti,di) = beam(2); w(ti,di) = beam(3);    coor = 'beam';
        %u(ti,di) = xyz(1); v(ti,di) = xyz(2); w(ti,di) = xyz(3);  coor = 'XYZ';
    end
end

%% scalar speed calculation from various coordinate sys options

beamangle = 25; % beam angle relative to vertical (25 deg, see manual)
linelength = 140; % moooring scope
maxmooringangle = atand(waterdepth/linelength)-10; % mooring (and sensor frame) angle relative to vertical
maxspeed = 3;

if coor(1:3) == 'enu' | coor(1:3) == 'xyz'
    speed = (u.^2 + v.^2 + w.^2).^.5;  % 3D pythagorean thm, only valid with orthogonal components (ENU or XYZ, not beam)
elseif coor(1:3) == 'bea'
    speed = abs( u + v ) ./ 2 ; % beams 1 & 2 should each point into the flow and measure negative along beam velocities.  Beam 3 may be contaminated.
    speed = speed ./ cosd(120/2); % projection from oblique cross-stream to streamwise (beams are 120 apart) 
    mooringangle = (speed(:,5)./maxspeed) .* maxmooringangle;  % ad hoc increase to full mooring angle based on speed
    speed = speed ./ sind(beamangle + mooringangle); % projection of beam 1 & 2 to horizontal 
    speedcheck = w./sind(beamangle - mooringangle); % projection of beam 3 to horizontal
end

%% component and speed plots (pcolors)

figure(5), clf
colormap(nawhimar)
subplot(3,1,1), pcolor(t, d, u'), shading flat, datetick, hold on, axis tight
set(gca,'YDir','reverse'), ylabel('depth [m]')
caxis([-2 2]), cb = colorbar; cb.Label.String = 'east [m/s]';
subplot(3,1,2), pcolor(t, d, v'), shading flat, datetick, hold on, axis tight
set(gca,'YDir','reverse'), ylabel('depth [m]')
caxis([-2 2]), cb = colorbar; cb.Label.String = 'north [m/s]';
subplot(3,1,3), pcolor(t, d, w'), shading flat, datetick, hold on, axis tight
set(gca,'YDir','reverse'), ylabel('depth [m]')
caxis([-2 2]), cb = colorbar; cb.Label.String = 'up [m/s]';
print -dpng SurfaceBuoy_AQD_Rosario_UVW.png

figure(6), clf
colormap(jet)
pcolor(t, d, speed'), shading flat, datetick, hold on, axis tight
set(gca,'YDir','reverse'), ylabel('depth [m]')
caxis([0 3.5]), cb = colorbar; cb.Label.String = 'horizontal speed [m/s]';
print -dpng SurfaceBuoy_AQD_Rosario_speed.png


%% tidal ellipses and u,v scatter

if corr(1:3)=='ENU'
    principalaxis = nanmean( atand( v ./ u ) );
    alld = ones(length(t),1)*d;
    
    figure(7), clf
    colormap copper
    scatter(u(:),v(:), 1, alld(:),'filled')
    cb = colorbar; cb.Label.String = 'depth [m]';
    xlabel('u, east [m/s]'), ylabel('v, north [m/s]')
    axis([-3 3 -3 3]), grid
    print -dpng SurfaceBuoy_AQD_Rosario_UVscatter.png
end


%% vertical profiles

lastbin = 17;

figure(8), clf
subplot(1,3,1), plot( smooth( nanmean(speed(:,1:lastbin)) ), d(1:lastbin) ), hold on, area([0 1.1],[3 3],0), 
plot([0 1.1],hubdepth*[1 1],'k--'), plot([0 1.1],(hubdepth+R)*[1 1],'k:'), plot([0 1.1],(hubdepth-R)*[1 1],'k:'), 
ylabel('depth [m]'), xlabel('avg speed [m/s]'), set(gca,'YDir','reverse'),
subplot(1,3,2), plot( smooth( 0.5 * 1030 * nanmean( speed(:,1:lastbin).^3 ) ), d(1:lastbin)), hold on, area([0 2000],[3 3],0), 
plot([0 2000],hubdepth*[1 1],'k--'), plot([0 2000],(hubdepth+R)*[1 1],'k:'), plot([0 2000],(hubdepth-R)*[1 1],'k:'), 
ylabel('depth [m]'), xlabel('avg avail. power density [W/m^2]'), set(gca,'YDir','reverse'),
% subplot(1,3,3), plot( 360+principalaxis, d), hold on, area([300 360],[3 3],0), 
% plot([0 1.1],hubdepth*[1 1],'k--'), plot([0 1.1],(hubdepth+R)*[1 1],'k:'), plot([0 1.1],(hubdepth-R)*[1 1],'k:'), 
% ylabel('depth [m]'), xlabel('principal axis [deg M]'), set(gca,'YDir','reverse'),set(gca,'XLim',[300 360])
print -dpng SurfaceBuoy_AQD_Rosario_profiles.png


%% AEP estimate for single O2 turbine, based on turbine parameters at top of script, 
% note that actual production would be double

% hubheight speed

[dz hubbin] = min( abs(d - hubdepth) ) ;
offset = 10; % guess at z adjustment for tilt.  Should not matter if profile is straight (minimal shear)
hubspeed = speed(:,hubbin-offset);
%hubspeed = filloutliers(hubspeed,'linear')

% annual definition
Tyear = 8760; % hours in a year
dt = median(diff(t)) * 24; % hourly spacing of data
T = range(t) * 24; % hours of valid data (103 days)

% 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 = ( Tyear / T ) * trapz( power ) * dt;
AEP_GWh = AEP./1e9;  % 1.9 GWh
AEP = 2*AEP % two rotors

%% Capacity Factor
ratedpower = mean( power(region3) );
CF = trapz(power) * dt ./ (T * ratedpower) % 0.17


%% speed and power histograms at hubheight

figure(9), clf
subplot(1,2,1), hist(hubspeed,[0:.1:5]), xlabel('Speed at hub depth [m/s]'), set(gca,'YLim',[0 1000]), ylabel('10 minute ensembles')
subplot(1,2,2), hist(power./1000,50), xlabel('turbine model output [kW]'),  set(gca,'YLim',[0 1000]), ylabel('10 minute ensembles')
print -dpng SurfaceBuoy_AQD_Rosario_hists.png


figure(10), clf
plot(t,hubspeed,'k-')
set(gca,'fontsize',16,'fontweight','demi'), datetick, ylabel('U_{hub} [m/s]')
print -dpng SurfaceBuoy_AQD_Rosario_hubspeed.png

%maxspeed = 


%% save clean output 
save SurfaceBuoy_AQD_Rosario.mat u v w t d hubspeed coor

[N bincenters] = hist(hubspeed,[0:.1:4]); output = [bincenters; N;]';
save -ASCII SurfaceBuoy_AQD_Rosario_HubSpeedHistogram.txt output

