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

% PURPOSE : Reading buoy directional moments, computing 2D spectrum & visualization
% HISTORY : - 2012 origination in matlab, F. Ardhuin
% - 01.2019 comments and rewriting, C. Peureux, cpeureux@ifremer.fr
% - 01.2019 updated for use of VideoWriter (avifile obsolete)
% REQUIREMENTS : CDIP_51201_2012_avg3h_freq.nc, readWWNC_FREQ.m, doppler_modified_land_2.ct
% MEM_calc.m, 
% DOC : Kuik et al. JPO 1988, Lygre & Krogstad JPO 1986

close all; % close all plotting windows
clear all; % clear all variables
clc; % clear console

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% 1. Load frequency spectrum and other spectral parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[lat,lon,freq,df,time,efi,th1mi,sth1mi,th2mi,sth2mi]=readWWNC_FREQ('CDIP_51201_2012_avg3h_freq.nc'); % from a buoy in Hawaii

%%% Compute spectral directional moments from Kuik et al.
ndates=521;
nf=size(efi,1); % Number of frequencies
efnow=efi(:,1,ndates); % Frequency spectrum [m^2/Hz]
th1now=th1mi(:,1,ndates); % Mean wave direction [deg] -> theta0 of Kuik et al.
sth1now=sth1mi(:,1,ndates); % Circular width parameter [deg]
th2now=th2mi(:,1,ndates); % Square root of mean wave squared direction [deg]
sth2now=sth2mi(:,1,ndates); % Square root of standard deviation of wave squared direction [deg]
d2r=pi/180; % degrres to radians conversion factor
m1=abs(1.-0.5*(sth1now(:)*d2r).^2);
m2=abs(1.-2*(sth2now(:)*d2r).^2);
a1=cos(th1now*d2r).*m1; % Directional distribution 2*pi periodic component (cosine component)
b1=sin(th1now*d2r).*m1; % Directional distribution 2*pi periodic component (sine component)
a2=cos(2.*th2now*d2r).*m2; % Directional distribution pi periodic component (cosine component)
b2=sin(2.*th2now*d2r).*m2; % Directional distribution pi periodic component (sine component)
Hs1 = 4*sqrt(sum(efnow.*df))

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% 2. Computes directional spectrum using MEM
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Maximum Entropy Method directional spectrum evaluation (Lygre and Krogstad
%%% 1986)
nth=72; % Number of directions for MEM calculation
[SD,Ef] = MEM_calc(a1,a2,b1,b2,efnow,nth); % input to MEM are line vectors. outputs are SD (directional distribution) and Ef (frequency spectrum for test))
Efth=SD.*repmat(efnow',nth,1)';
x1=0:(360/nth):360-360/nth;  
nth=size(Efth,2);
dir=x1'.*d2r;
%%% Reads the color palette for pretty pictures ... 
col=load('doppler_modified_land_2.ct');
col(:,254)=col(:,253);
col(:,255)=col(:,253);
col(:,256)=col(:,253);
%%% Changes dir from nautical (0 is North, clockwise) to trigonometric convention
dirmemo=dir;
dir=pi/2.-dir;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% 3. Plots
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% Plots the spectrum in polar coordinates
%%% Repeats the 1st direction at the end to wrap around nicely 
dirr=[dir' dir(1)]';
dir2=repmat(dirr',nf,1);
freq2=repmat(freq,1,nth+1);
Efth2=zeros(nf,nth+1);
Efth2(:,1:nth)=Efth;
Efth2(:,nth+1)=Efth(:,1);
x2=-cos(dir2).*freq2; % the - transforms "direction to" in "direction from"
y2=-sin(dir2).*freq2; % the - transforms "direction to" in "direction from"
%%% Plots 2D spectrum
figure(1);
set(gca,'FontSize',15);
colormap(col'./255);
pcolor(x2,y2,Efth2);axis equal;shading flat;
colorbar;
hold on;

for i=1:7 
    plot(0.1*i*cos(linspace(0,2*pi,25)),0.1*i*sin(linspace(0,2*pi,25)))
end
title('Spectrum estimated from MEM method');

%%% Plots 1D spectrum in frequency only: E(f)
figure(2);
clf
dth=2.*pi./real(nth);

Ef=sum(Efth,2)*dth;
plot(freq,Ef,'k-+','LineWidth',2);
set(gca,'FontSize',15);
xlabel('Frequency (Hz)');
ylabel('E(f) (m^2/Hz)');

%%% Computes and plots mean directions
dir=dirmemo;
dir2=repmat(dir',nf,1);
a1=sum(Efth.*cos(dir2),2).*dth./Ef;
b1=sum(Efth.*sin(dir2),2).*dth./Ef;

m1=sqrt(a1.^2+b1.^2);
sth1m=sqrt(2*(1-m1))*180/pi;

figure(5);
clf
th1m=mod(180-atan2(b1,a1)*180/pi,360);
plot(freq,th1m,'k-',freq,th1mi(:,1,ndates),'ro','LineWidth',2);
set(gca,'FontSize',15);
xlabel('Frequency (Hz)');
ylabel('Mean direction th1m (deg)');
legend('After MEM','before MEM')

figure(6);
clf
plot(freq,sth1m,'k-',freq,sth1mi(:,1,ndates),'ro','LineWidth',2);
set(gca,'FontSize',15);
xlabel('Frequency (Hz)');
ylabel('Spread sth1m (deg)');
legend('After MEM','before MEM')




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% 4. Computes and plots amps of surface elevation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% A surface elevation map is snthetized which frequency-direction
%%% spectrum is Efth

%%% Random draw of phases ... 
phases=rand(nf,nth)*2*pi;
nk=nf;
%%%  Amplitudes from frequency spectrum
df=freq.*(1.1-1/1.1)./2; % Frequency spacing (standar modelled frequency)
dth=2.*pi./real(nth); % Angular spacing
as=zeros(nk,nth); % Amplitudes [m]
for i=1:nf
    as(i,:)=sqrt(2.*Efth(i,:)*df(i)*dth); % By definition of the spectrum
end

%%% Define the area where the surface elevation is reconstructed
nx=201;
x=linspace(0,1000,nx);
y2=repmat(x',1,nx);
x2=repmat(x,nx,1);
y=x;

%%% Computes wave number
sig=2*pi.*freq;
g=9.81; % gravity
ks=sig.^2./g; % from dispersion relation

%%% Computes time vector for reconstruction
nt=20;
t=linspace(0,(nt-1)*0.5,nt);

%%% Changing color map to something nicer
col=load('doppler_modified_land_2.ct');
col(:,254)=col(:,253);
col(:,255)=col(:,253);
col(:,256)=col(:,253);
colormap(col'./255);

%%% Plot successive maps of surface elevation
figure(3);
clf
set(3,'Position',[1 1 3*nx+40 3*nx+40])
colormap(col'./255);
M=struct([]);
mov = VideoWriter('example_51001.avi');
open(mov) ;
for ii=1:nt
%
% initialise la surface a zero
%
zeta=zeros(nx,nx);
hs=0.;
E=0;
% indices 1->8, swell 1
% indices 9->14 swell 2
for i=1:nk
    for j=1:nth
        zeta(:,:)= zeta(:,:)+as(i,j)*cos(ks(i)*cos(pi/2-dir(j)).*x2 ...
            +ks(i)*sin(pi/2-dir(j)).*y2 + phases(i,j)-sig(i)*t(ii));
    E=E+as(i,j).^2./2.;
    end
end
pcolor(x2,y2,zeta);
caxis([-2.5,2.5]);
axis equal;shading interp;colorbar;
xlabel('x (m)');
ylabel('y (m)');
N(ii)=getframe;
F=getframe(gcf);
writeVideo(mov,F)
end

%%% Check that the value of Hs did not change (should be the same as Hs1)
Hs2=4*sqrt(E)

