clear 
close all
warning('off','all')
%
g = 9.8; rho = 1.025e3;
%
compi = complex(0,1);
compi = 1i;
%
set(0,'DefaultAxesFontName','Arial')
set(0,'DefaultAxesFontSize',16)
set(0,'DefaultAxesColorOrder',[1 0 0;0 0 1; 0 0 0],'DefaultAxesLineStyleOrder','-|--|-.|:|.')
%
N = 128*2*2;
Nm = 4*2;
%
Nt = N*2*2*2;
Nf = N;
%
Mf = Nf;
%
oml = 0.18;omh = 1.57;
%
omeg = zeros(Nf,1);
t = zeros(Nt,1);
%
delom = (omh-oml)/double(Nf);
%
tend = 600.0;
tstart = 0.0;
%
delt = (tend-tstart)/double(Nt);
%
Hs = 1.0;
Te = 9.0;
%
% Sept. 28, 2020; ~ 3 pm; MBARI Paragon measurements.
Hs = 2.7;
Te = 9.5;
%
omeg = zeros(Mf,1);
%
for m = 1:Mf
    omeg(m,1) = oml + double(m-1)*delom;
end
%
t = zeros(Nt,1);
%
for n = 1:Nt
    t(n,1) = double(n-1)*delt;
end
%
%
% Pieson_Moskowitz spectrum
%
multcof = 2.0*131.5*Hs^2/Te^4;
multarg = 1054.0/(Te^4);
%
%etacf = zeros(M,1);
%
spec = zeros(Mf,1);
etaAf = zeros(Mf,1);
%
for m = 1:Mf
    scof = multcof/(omeg(m,1)^5);
    st = exp(-multarg/(omeg(m,1)^4));
    spec(m,1) = scof*st;
end %m
%
% phase = 2*pi*rand(Mf,1);
% %
% amp = zeros(Mf,1);
% %
% for m = 1:Mf;
%     amp(m,1) = sqrt(2*spec(m,1)*delom); % use if simply summing varous frequency components to form irregular signal, not inv FT as here.    
%  %   amp(m,1) = 1.0;
%     %   amp(m,1) = sqrt(spec(m,1)/2); % Inverse Fourier transformation with division by delt accounts for dom here.
% end %m
%
etaAt = zeros(Nt,1);
% %
% for n = 1:Nt
%     sum = 0.0;
xA = 1000.0;
xB = 0.0;
%
Nx = Nt;
delx = (xA-xB)/double(Nx);
%
for n = 1:Nx
    x(n,1) = double(n-1)*delx;
end
%
for m = 1:Mf
    k(m,1) = omeg(m,1)^2/g;
%     etaf1 = exp(compi*phase(m,1));
%     etaAf(m,1) = amp(m,1)*exp(-compi*(k(m,1)*xA))*etaf1;
end
%
% for n = 1:Nt
%     sum = 0.0;
%     for m = 1:Mf
%         sum = sum + real(etaAf(m,1)*exp(compi*omeg(m,1)*t(n,1)));
%     end %m
%     etaAt(n,1) = sum;
% end
% %
%
% figure(1)
% plot(t,etaAt,'-','LineWidth',3)
% xlabel('Time t(s)','FontSize',16)
% ylabel('\eta_A (m)','FontSize',16)
% grid on
%
% %
% load crests_t.csv;
% fname1 = crests_t;
% %
% etaMt = zeros(Mf,1);
% %
% %tm = 100.0;
% tm1 = 100.0;
% %
% for m = 1:Mf
%     etaf = exp(compi*phase(m,1));
%     etaMt(m,1) = amp(m,1)*exp(compi*omeg(m,1)*tm1)*etaf;
% end %m
% %
% etaMx = zeros(Nx,1);
% %
% for n = 1:Nx
%     sum = 0.0;
%     for m = 1:Mf
%         sum = sum + real(etaMt(m,1)*exp(-compi*k(m,1)*x(n,1)));
%     end
%     etaMx(n,1) = sum;
% end
% %
% fname1 = 'crest_meas1';
% figure(2)
% plot(x,etaMx,'-r','LineWidth',3)
% xlabel('Distance x (m)','FontSize',16)
% ylabel('\eta_M (m)','FontSize',16)
% grid on
% %
% saveas(gcf,fname1,'epsc')
% saveas(gcf,fname1,'jpeg')
%
Mt = 8;
tm = zeros(Mt,1);
%
delt = 2.5 + 0.6; %time between scans + 1/2 the time needed for 1 scan.
delt = 1.26;
for mt = 1:Mt
    tm(mt,1) = double(mt-1)*delt;
end %mt
%
load crests_x_mbari_multi.csv;
fname1 = crests_x_mbari_multi;
%
%Mc = length(fname1);
Mc = 3;
xc = zeros(Mc,Mt,1);
yc = zeros(Mc,Mt,1);
%
for m = 1:Mc
    for mt = 1:Mt
    xc(m,mt,1) = -fname1(m,mt);
    yc(m,mt,1) = 0.5;
    end % mt
end %m
%
Ac = zeros(Mf,1);
Bc = zeros(Mf,1); % set to unity
alpha = zeros(Mf,1); % phases
C = zeros(Mf,1); % amplitudes
%

for n = 1:Mf
    sumnum = 0.0;
    sumden = 0.0;
    for m = 1:Mc
        sumnt = 0.0;
        %sumdt = 0.0;
        for mt = 1:Mt
        cterm = cos(k(n,1)*xc(m,mt,1) - omeg(n,1)*tm(mt,1));
        sumnt = sumnt + cterm;
        end %mt
        sumnum = sumnum + sumnt/double(Mt);
    end %m
    %
    for m = 1:Mc
        sumdt = 0.0;
        for mt = 1:Mt
        sterm = sin(k(n,1)*xc(m,mt,1) - omeg(n,1)*tm(mt,1));
        sumdt = sumdt + sterm;
        end %mt
        sumden = sumden + sumdt/double(Mt);
    end
    %
    Ac(n,1) = sumnum/sumden;
    alpha(n,1) = atan2(1.0,Ac(n,1));
end %n
%
% spectrum used above is borrowed below.
for n = 1:Mf
    C(n,1) = sqrt(2*spec(n,1)*delom);
 %   C(n,1) = amp(n,1);
end
%
%
fname1 = 'phases_frm_crests.csv';
csvwrite(fname1,alpha);
% %
% etam = zeros(Nx,1);
% for n = 1:Nx
%     sumx = 0.0;
%     for m = 1:Mf
%  %       for mt = 1:Mt
%         wprof = exp(-compi*(k(m,1)*x(n,1) - omeg(m,1)*tm(1,1) + alpha(m,1)));
%         wprof1 = cos(k(m,1)*x(n,1) - omeg(m,1)*tm(1,1) - alpha(m,1));
%         sumx = sumx + C(m,1)*real(wprof);
%         sumx = sumx + C(m,1)*wprof1;
%     end %m
%     etam(n,1) = sumx;
% end %n
% %
% %
% fname1 = 'wave_full_from_crests';
% figure(3)
% plot(x,etaMx,'-r',x,etam/4.0,'--b','LineWidth',3)
% xlabel('x (m)','FontSize',16)
% ylabel('\eta_x, \eta_m (m)','FontSize',16)
% grid on
% legend('crest locations from','full wave','location','best')
% saveas(gcf,fname1,'epsc')
% saveas(gcf,fname1,'jpeg')
% %
% %
% fname1 = 'wave_crest_loc_wave_elev';
% figure(4)
% plot(xc,yc,'*b',x,etam/4.0,'-r','LineWidth',3)
% xlabel('{\it x} (m)','FontSize',16)
% ylabel('Crests, \eta_m (m)','FontSize',16)
% grid on
% legend('crest locations','fitted full wave','location','best')
% saveas(gcf,fname1,'epsc')
% saveas(gcf,fname1,'jpeg')
% %
% 
%     



    