clear; clc; close all; 

%% Carga de datos

x=[0 x12 x13 x14 x15]; %Ingrese las distancias de los wavegauges al primer sensor. Ej.[0 1.3 3 4 5.5]
h=1; %Profundidad
rate2sample=freq; %frecuencia de muestreo

SL=; %Cargar matriz con las medidas de nivel de todos los sensores. Primer sensor en la primer columna, etc. 
t=; %Cargar vector de tiempo asociado a la serie de oleaje

%% Descompone entre incidente y reflejada

ARR=[1 2 3]; % El metodo usa 3 sensores. Elegir un arreglo adecuado de acuerdo a la longitud de onda (Mansard & Funke, 1980) . Ej. 1-2-3; 1-3-4 y  2-3-5
Zi_tot=[];Zr_tot=[];

X1=x(ARR(1)); X2=x(ARR(2)); X3=x(ARR(3)); 

eta=SL(:,ARR)/100; x2=[0 X1-X2 X1-X3];
    
if mod(length(eta), 2) == 1 %Corrección por si el numero de datos es impar
    etax=[]
    for i=ARR
     aux=eta(:,i); aux2=aux(length(eta))
    aux(length(eta)+1)=aux2;
    etax=[etax aux];
    end
    eta=etax;
end
    

dt=1/rate2sample; eps=1e-10;f1=0.05; f2=4;
[f,Zi,Zr]= clsquare_s(eta,x2,h,dt,eps,f1,f2);
[etai,etar,eta2]  = regen_s(f,Zi,Zr,x2,h,25); 


figure %ejemplo
t2=0:1/rate2sample:(numel(SL(:,1))-1)*1/rate2sample;
plot(t2,eta(:,1),'b',t2,eta2(:,1),'r')
legend('Medida','Reconstruida');ylabel('SL (m)');xlabel('t (s)');title('Sensor1')


%% Calcula potencia incidente y reflejada

rho=1000; g=9.81; b=1.5;
for i=1:numel(f)
    ki(i)= wavnum(1/f(i),h);
    kr(i)= wavnum(1/f(i),h);
    ci=sqrt(g/ki(i)*tanh(ki(i)*h));
    cr=sqrt(g/kr(i)*tanh(kr(i)*h));
    Cgi=0.5*(1+(2*ki(i)*h)/sinh(2*ki(i)*h))*ci;
    Cgr=0.5*(1+(2*kr(i)*h)/sinh(2*kr(i)*h))*cr;
    Pi(i)=1/2*(2*(abs(Zi(i))))^2*rho*g*Cgi*b;
    Pr(i)=1/2*(2*(abs(Zr(i))))^2*rho*g*Cgr*b;
end
%
Pi_tot=nansum(Pi);
Pr_tot=nansum(Pr);



