%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% WAVES 1 - LAB 1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% PURPOSE : plot Hs time series from two different types of measurements
% HISTORY : - origination in matlab
% - 01.2019 comments and rewriting, C. Peureux, cpeureux@ifremer.fr
% REQUIREMENTS : readWWNC_FREQ.m

clear all ;
close all ;

%%% PARAMETER
maxfreq = 0.1 ; % max frequency for Hs calculation [Hz]

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% 1. Computes Hs from buoy data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Reads buoy data
filebuoy='CDIP_51201_2012_avg3h_freq.nc' ;
[latb,lonb,freqb,dfb,timeb,efb,th1mb,sth1mb,th2mb,sth2mb] = readWWNC_FREQ(filebuoy) ;

%%% Computes Hs
ntb = length(timeb) ; % number of times
dfb2d = repmat(dfb,[1,ntb]) ;
efb = squeeze(efb) ;
Ib = freqb < maxfreq ;
Hsb = 4*sqrt(sum(efb(Ib,:).*dfb2d(Ib,:),1)) ;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% 2. Computes Hs from pressure gauge data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Reads buoy data
filepg='WAIMEA_2012_freq.nc' ;
[latg,long,freqg,dfg,timeg,efg] = readWWNC_FREQ(filepg) ;

%%% Computes Hs
ntg = length(timeg) ; % number of times
dfb2g = repmat(dfg,[1,ntg]) ;
efg = squeeze(efg) ;
Ig = freqg < maxfreq ;
Hsg = 4*sqrt(sum(efg(Ig,:).*dfb2g(Ig,:),1)) ;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% 3. Plots
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
t0d = datetime(2012,1,19) ;
t0 = datenum(t0d) ;
figure('position',[250 250 1000 200]) ;
plot(timeb-t0,Hsb,'-','LineWidth',2) ;
hold on ;
plot(timeg-t0,Hsg,'+') ;
legend('Buoy','Pressure gauge') ;
xlabel(['Days since ' num2str(year(t0d)) '-' num2str(month(t0d)) '-' num2str(day(t0d)) ' 00:00:00']) ;
xlim([0 60])
grid('on') ;
ylabel(['Hs (f<' num2str(maxfreq) 'Hz) [m]'])

%%% ANSWER TO LAB1
fsamp =  0.01:0.01:0.5  ;
nrmse = zeros(length(fsamp),1) ;
bias = zeros(length(fsamp),1) ;
for ifreq=1:length(fsamp)
    Ib = freqb<=fsamp(ifreq) ;
    hsb = 4*sqrt(trapz(freqb(Ib),dfb2d(Ib,:).*efb(Ib,:),1)) ;
    hsb = interp1(timeb,hsb,timeg) ;
    Ig = freqg<=fsamp(ifreq) ;
    hsg = 4*sqrt(trapz(freqg(Ig),dfb2g(Ig,:).*efg(Ig,:),1)) ;
    nrmse(ifreq) = sqrt(sum((hsb-hsg').^2)/sum(hsb.^2)) ;
    bias(ifreq) = sum(hsg'-hsb)/sum(hsb) ;
end

figure ;
plot(fsamp,bias*100,'o-')
hold on ;
plot(fsamp,nrmse*100,'+-')
grid on ;
xlabel('Frequency [Hz]')
ylabel('Errors [%]')
legend('Bias','NRMSE')
xlim([0.05,0.15])

