%define basic variables
clear;
clc;
Rho=1030;
g=9.81;
Avolt=253;
Bvolt=-44000;
load rawdata_121002_2.txt;
rawdata=rawdata_121002_2;
%transformation from Volts to pressure using calibration constants Avolt
%and Bvolt
%Avolt has dimension Pa/V, Bvolt has dimension Pa

P=Avolt*rawdata+Bvolt;  
%interval is the sample frequency interval of the sensor
interval=.25;
n=numel(rawdata); %count number of samples
tottime=n*interval; %calculate total duration of observation in sec
ttime=(interval:interval:tottime); %create array with time
time=ttime.'; %change orientation of matrix
% calculate regression coefficients to compensate for change in waterlevel
% during the observations
%plot(time,P);
B = [ones(length(time),1) time] \ P;
Slope=B(2);
Intercept=B(1);
Pwave=P-Intercept'-time*Slope;

figure
plot(time,Pwave);
xlabel('time (sec)');
ylabel('pressure (Pa)');
title('Pressure time series Asparuhovo beach 2');

Pstatic=P-Pwave;
Depth=Pstatic/Rho/g;
AverageDepth=mean(Depth);

%figure
%plot(time,Depth);
%ylabel('depth (m)');
%xlabel('time (sec)');

counter=0;
Pmin=0;
Pmax=0;
period=0;
for k=1:n-1
    if(Pwave(k+1)*Pwave(k)<0 && Pwave(k+1)<0)  %new wave starts at k+1
        counter=counter+1;
        T(counter)=period;
        L0=g/2/pi*period*period;
        if (Depth(k)/L0 < 0.36) 
            L=sqrt(g*Depth(k))*(1-Depth(k)/L0)*period;
        else
            L=L0;
        end;
        H(counter)=(Pmax-Pmin)/Rho/g*cosh(2*pi/L/Depth(k));
        Pmax=0;
        Pmin=0;
        period=0;
    end;
    if(Pwave(k)>Pmax)  Pmax=Pwave(k);
    end;
    if (Pwave(k)<Pmin) Pmin=Pwave(k);
    end;
    period=period+interval;
end;
%counter = number of waves found
Hsorted=sort(H);
kk=0;
for k=1:counter
    if  Hsorted(k)>0.00
        kk=kk+1;
        Hsorted2(kk)=Hsorted(k);
        
    end;
end;
counter=kk;

Tsorted=sort(T);
Tmean = mean (Tsorted);
%loop to estimate water level variations
for k=1:n
        period=Tmean;
        L0=g/2/pi*period*period;
        if (Depth(k)/L0 < 0.36) 
            L=sqrt(g*Depth(k))*(1-Depth(k)/L0)*period;
        else
            L=L0;
        end;
        eta(k)=(Pwave(k))/Rho/g*cosh(2*pi/L/Depth(k));
end;

figure
plot(time,eta);
xlabel('time (sec)');
ylabel('waterlevel (m)');
title('Time series Asparuhovo beach 2');
%end waterlevel loop

%start loop to determine individual waves from waterlevel
etamin=0;
etamax=0;
counter=0;
for k=1:n-1
    if(eta(k+1)*eta(k)<0 && eta(k+1)<0)  %new wave starts at k+1
        counter=counter+1;
        Hzz(counter)=etamax-etamin;
        etamin=0;
        etamax=0;
    end;
    if(eta(k)>etamax)  etamax=eta(k);
    end;
    if (eta(k)<etamin) etamin=eta(k);
    end;
end;
%end loop to determine individual waves
Hsorted3=sort(Hzz);
kk=0;
for k=1:counter
    if  Hsorted3(k)>0.00
        
        kk=kk+1;
        Hsorted2(kk)=Hsorted3(k);
    end;
end;
counter=kk;

for k=1 : counter    
prob(k)= (counter-k+1)/counter;
end;

figure
probplot('rayleigh',Hsorted2(1:kk));
xlabel('wave height (m)');
title('Rayleigh Asparuhovo beach 2');

Hmean=mean(Hsorted2);
Hrms = sqrt(mean(Hsorted2.^2)) % rms wave height
sum = 0;
from=2*counter/3;
for k=from:counter
    sum=sum+Hsorted2(counter);
end
Hs=sum/(counter/3) % significant wave height
TanAlfa = 0.001;   %bedslope
Htr=(0.35+5.8*TanAlfa)*AverageDepth % Transistion wave height acc. to BG
HtrNorm=Htr/Hrms   %normalised transition height
H1Norm= 1.052;
H2Norm= 1.175;
H01Norm=1.796;
H001Norm=2.011;
H1=  H1Norm*Hrms
H2=  H2Norm*Hrms
H01= H01Norm*Hrms
H001=H001Norm*Hrms
PrHtr= 1-exp(-(Htr/H1)^2)
PrHtrx=1-exp(-(Htr/H2)^3.6)
PrH01= 1-exp(-(H01/H2)^3.6)
PrH001= 1-exp(-(H001/H2)^3.6)