%SCRIPS PARA ANÁLISIS ARMÓNICO DE MAREA
m_posorja=xlsread('mareas_apg.xlsx','Posorja');
m_Gquil=xlsread('mareas_apg.xlsx','Guayaquil 1');
m_Enfilada=xlsread('mareas_apg.xlsx','Enfilada C2');

t_pos=datenum(m_posorja(:,1:6));
t_gye=datenum(m_Gquil(:,1:6));
t_C2=datenum(m_Enfilada(:,1:6));

scrsz = get(0,'ScreenSize');

% Define inference parameters.
%       infername=['P1';'K2'];
%       inferfrom=['K1';'S2'];
%      infamp=[.33093;.27215];
%      infphase=[-7.07;-22.40];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%ANÁLISIS ARMÓNICO DE POSORJA%

% The call (see t_demo code for details).
[tidestruc_pos,pout_pos]=t_tide(m_posorja(:,end),...
    'interval',1, ...                     % hourly data
    'start',t_pos(1),...               % start time is datestr(tuk_time(1))
    'latitude',-2.7061,...               % Latitude of obs
    'synthesis',1);                       % Use SNR=1 for synthesis.

figure('Position',[1 1 (3.*scrsz(3))/5 (5.*scrsz(4))/6])
%clf;orient tall;
subplot(211);
plot(t_pos,[m_posorja(:,end) pout_pos]);
line(t_pos,m_posorja(:,end)-pout_pos,'linewi',2,'color','r');
set(gca,'Fontsize',12)
%xlabel('Days in 1975');
ylabel('Altura (m)','FontSize',12,'FontWeight','bold');
gtext('Serie de tiempo Original','color','b');
gtext('Predicción de marea','color',[0 .5 0]);
gtext('Residual','color','r');
title('Análisis Armónico - Posorja','FontSize',18);
%ajustar ejes y aplicar datetick
datetick('x',1)

subplot(212);
fsig=tidestruc_pos.tidecon(:,1)>tidestruc_pos.tidecon(:,2); % Significant peaks
semilogy([tidestruc_pos.freq(~fsig),tidestruc_pos.freq(~fsig)]',[.0005*ones(sum(~fsig),1),tidestruc_pos.tidecon(~fsig,1)]','.-r');
line([tidestruc_pos.freq(fsig),tidestruc_pos.freq(fsig)]',[.0005*ones(sum(fsig),1),tidestruc_pos.tidecon(fsig,1)]','marker','.','color','b');
line(tidestruc_pos.freq,tidestruc_pos.tidecon(:,2),'linestyle',':','color',[0 .5 0]);
set(gca,'ylim',[.0005 1],'xlim',[0 .5],'Fontsize',12);
xlabel('frequencia (cph)','FontSize',12,'FontWeight','bold');
text(tidestruc_pos.freq,tidestruc_pos.tidecon(:,1),tidestruc_pos.name,'rotation',45,'vertical','base');
ylabel('Amplitud (m)','FontSize',12,'FontWeight','bold');
title('Líneas analizadas con el 95% de confianza','FontSize',12,'FontWeight','bold'); %titulo del subplot
gtext('Componentes Significantes','color','b');
gtext('Componentes no significantes','color','r');
gtext('Nivel del 95% Confianza','color',[0 .5 0]);

print('-dpng','-r400','armonic_pos1.png')

%%%%%%%%%%%%%%%%%%%%%
figure('Position',[1 1 (3.*scrsz(3))/5 (5.*scrsz(4))/6])

subplot(211);
errorbar(tidestruc_pos.freq(~fsig),tidestruc_pos.tidecon(~fsig,3),tidestruc_pos.tidecon(~fsig,4),'.r');
hold on;
errorbar(tidestruc_pos.freq(fsig),tidestruc_pos.tidecon(fsig,3),tidestruc_pos.tidecon(fsig,4),'o');
hold off;
set(gca,'ylim',[-45 360+45],'xlim',[0 .5],'ytick',[0:90:360],'FontSize',12);
xlabel('frequencia (cph)','FontSize',12,'FontWeight','bold');
ylabel('Fase de Greenwich (grados)','FontSize',12,'FontWeight','bold');
title('Ángulos de fase analizados con 95% IC','FontSize',12,'FontWeight','bold');
gtext('Componentes Significantes','color','b');
gtext('Componentes no Significantes','color','r');

subplot(212);
ysig=m_posorja(:,end);
yerr=m_posorja(:,end)-pout_pos;
%nfft=389;
nfft=2^(nextpow2(length(m_posorja(:,end))));

bd=isnan(ysig);
gd=find(~bd);
bd([1:(min(gd)-1) (max(gd)+1):end])=0;
ysig(bd)=interp1(gd,ysig(gd),find(bd)); 
[Pxs,F]=psd(ysig(finite(ysig)),nfft,1,[],ceil(nfft/2));
%[Pxs,F]=pmtm(ysig(finite(ysig)),4,4096,1);
yerr(bd)=interp1(gd,yerr(gd),find(bd)); 
[Pxe,F]=psd(yerr(finite(ysig)),nfft,1,[],ceil(nfft/2));
%[Pxe,F]=pmtm(yerr(finite(ysig)),4,4096,1);

semilogy(F,Pxs);
set(gca,'Fontsize',12)
line(F,Pxe,'color','r');
xlabel('frequencIA (cph)','FontSize',12,'FontWeight','bold');
ylabel('m^2/cph','FontSize',12,'FontWeight','bold');
title('Espectro Estimado antes y después de remover la energía espectral','FontSize',12,'FontWeight','bold');
gtext('Serie Original (interpolada)','color','b');
gtext('Energía No-mareal analizada','color','r');

print('-dpng','-r400','armonic_pos2.png')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%ANÁLISIS ARMÓNICO DE Puerto de Guayaquil%

% The call (see t_demo code for details).
[tidestruc_gye,pout_gye]=t_tide(m_Gquil(:,end),...
    'interval',1, ...                     % hourly data
    'start',t_gye(1),...               % start time is datestr(tuk_time(1))
    'latitude',-2.255946,...               % Latitude of obs
    'synthesis',1);                       % Use SNR=1 for synthesis.

figure('Position',[1 1 (3.*scrsz(3))/5 (5.*scrsz(4))/6])
%clf;orient tall;
subplot(211);
plot(t_gye,[m_Gquil(:,end) pout_gye]);
line(t_gye,m_Gquil(:,end)-pout_gye,'linewi',2,'color','r');
set(gca,'Fontsize',12)
%xlabel('Days in 1975');
ylabel('Altura (m)','FontSize',12,'FontWeight','bold');
gtext('Serie de tiempo Original','color','b');
gtext('Predicción de marea','color',[0 .5 0]);
gtext('Residual','color','r');
title('Análisis Armónico - Puerto de Guayaquil','FontSize',18);
%ajustar ejes y aplicar datetick
datetick('x',1)

subplot(212);
fsig=tidestruc_gye.tidecon(:,1)>tidestruc_gye.tidecon(:,2); % Significant peaks
semilogy([tidestruc_gye.freq(~fsig),tidestruc_gye.freq(~fsig)]',[.0005*ones(sum(~fsig),1),tidestruc_gye.tidecon(~fsig,1)]','.-r');
line([tidestruc_gye.freq(fsig),tidestruc_gye.freq(fsig)]',[.0005*ones(sum(fsig),1),tidestruc_gye.tidecon(fsig,1)]','marker','.','color','b');
line(tidestruc_gye.freq,tidestruc_gye.tidecon(:,2),'linestyle',':','color',[0 .5 0]);
set(gca,'ylim',[.0005 1],'xlim',[0 .5],'Fontsize',12);
xlabel('frequencia (cph)','FontSize',12,'FontWeight','bold');
text(tidestruc_gye.freq,tidestruc_gye.tidecon(:,1),tidestruc_gye.name,'rotation',45,'vertical','base');
ylabel('Amplitud (m)','FontSize',12,'FontWeight','bold');
title('Líneas analizadas con el 95% de confianza','FontSize',12,'FontWeight','bold'); %titulo del subplot
gtext('Componentes Significantes','color','b');
gtext('Componentes no significantes','color','r');
gtext('Nivel del 95% Confianza','color',[0 .5 0]);

print('-dpng','-r400','armonic_gye1.png')

%%%%%%%%%%%%%%%%%%%%%
figure('Position',[1 1 (3.*scrsz(3))/5 (5.*scrsz(4))/6])

subplot(211);
errorbar(tidestruc_gye.freq(~fsig),tidestruc_gye.tidecon(~fsig,3),tidestruc_gye.tidecon(~fsig,4),'.r');
hold on;
errorbar(tidestruc_gye.freq(fsig),tidestruc_gye.tidecon(fsig,3),tidestruc_gye.tidecon(fsig,4),'o');
hold off;
set(gca,'ylim',[-45 360+45],'xlim',[0 .5],'ytick',[0:90:360],'FontSize',12);
xlabel('frequencia (cph)','FontSize',12,'FontWeight','bold');
ylabel('Fase de Greenwich (grados)','FontSize',12,'FontWeight','bold');
title('Ángulos de fase analizados con 95% IC','FontSize',12,'FontWeight','bold');
gtext('Componentes Significantes','color','b');
gtext('Componentes no Significantes','color','r');

subplot(212);
ysig=m_Gquil(:,end);
yerr=m_Gquil(:,end)-pout_gye;
%nfft=389;
nfft=2^(nextpow2(length(m_Gquil(:,end))));

bd=isnan(ysig);
gd=find(~bd);
bd([1:(min(gd)-1) (max(gd)+1):end])=0;
ysig(bd)=interp1(gd,ysig(gd),find(bd)); 
[Pxs,F]=psd(ysig(finite(ysig)),nfft,1,[],ceil(nfft/2));
%[Pxs,F]=pmtm(ysig(finite(ysig)),4,4096,1);
yerr(bd)=interp1(gd,yerr(gd),find(bd)); 
[Pxe,F]=psd(yerr(finite(ysig)),nfft,1,[],ceil(nfft/2));
%[Pxe,F]=pmtm(yerr(finite(ysig)),4,4096,1);

semilogy(F,Pxs);
set(gca,'Fontsize',12)
line(F,Pxe,'color','r');
xlabel('frequencIA (cph)','FontSize',12,'FontWeight','bold');
ylabel('m^2/cph','FontSize',12,'FontWeight','bold');
title('Espectro Estimado antes y después de remover la energía espectral','FontSize',12,'FontWeight','bold');
gtext('Serie Original (interpolada)','color','b');
gtext('Energía No-mareal analizada','color','r');

print('-dpng','-r400','armonic_gye2.png')


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%ANÁLISIS ARMÓNICO DE ENFILADA  C2

% The call (see t_demo code for details).
[tidestruc_C2,pout_C2]=t_tide(m_Enfilada(:,end),...
    'interval',1, ...                     % hourly data
    'start',t_C2(1),...               % start time is datestr(tuk_time(1))
    'latitude',-2.586314,...               % Latitude of obs
    'synthesis',1);                       % Use SNR=1 for synthesis.

figure('Position',[1 1 (3.*scrsz(3))/5 (5.*scrsz(4))/6])
%clf;orient tall;
subplot(211);
plot(t_C2,[m_Enfilada(:,end) pout_C2]);
line(t_C2,m_Enfilada(:,end)-pout_C2,'linewi',2,'color','r');
set(gca,'Fontsize',12)
%xlabel('Days in 1975');
ylabel('Altura (m)','FontSize',12,'FontWeight','bold');
gtext('Serie de tiempo Original','color','b');
gtext('Predicción de marea','color',[0 .5 0]);
gtext('Residual','color','r');
title('Análisis Armónico - Enfilada C2','FontSize',18);
%ajustar ejes y aplicar datetick
datetick('x',1)

subplot(212);
fsig=tidestruc_C2.tidecon(:,1)>tidestruc_C2.tidecon(:,2); % Significant peaks
semilogy([tidestruc_C2.freq(~fsig),tidestruc_C2.freq(~fsig)]',[.0005*ones(sum(~fsig),1),tidestruc_C2.tidecon(~fsig,1)]','.-r');
line([tidestruc_C2.freq(fsig),tidestruc_C2.freq(fsig)]',[.0005*ones(sum(fsig),1),tidestruc_C2.tidecon(fsig,1)]','marker','.','color','b');
line(tidestruc_C2.freq,tidestruc_C2.tidecon(:,2),'linestyle',':','color',[0 .5 0]);
set(gca,'ylim',[.0005 1],'xlim',[0 .5],'Fontsize',12);
xlabel('frequencia (cph)','FontSize',12,'FontWeight','bold');
text(tidestruc_C2.freq,tidestruc_C2.tidecon(:,1),tidestruc_C2.name,'rotation',45,'vertical','base');
ylabel('Amplitud (m)','FontSize',12,'FontWeight','bold');
title('Líneas analizadas con el 95% de confianza','FontSize',12,'FontWeight','bold'); %titulo del subplot
gtext('Componentes Significantes','color','b');
gtext('Componentes no significantes','color','r');
gtext('Nivel del 95% Confianza','color',[0 .5 0]);

print('-dpng','-r400','armonic_enfilad1.png')

%%%%%%%%%%%%%%%%%%%%%
figure('Position',[1 1 (3.*scrsz(3))/5 (5.*scrsz(4))/6])

subplot(211);
errorbar(tidestruc_C2.freq(~fsig),tidestruc_C2.tidecon(~fsig,3),tidestruc_C2.tidecon(~fsig,4),'.r');
hold on;
errorbar(tidestruc_C2.freq(fsig),tidestruc_C2.tidecon(fsig,3),tidestruc_C2.tidecon(fsig,4),'o');
hold off;
set(gca,'ylim',[-45 360+45],'xlim',[0 .5],'ytick',[0:90:360],'FontSize',12);
xlabel('frequencia (cph)','FontSize',12,'FontWeight','bold');
ylabel('Fase de Greenwich (grados)','FontSize',12,'FontWeight','bold');
title('Ángulos de fase analizados con 95% IC','FontSize',12,'FontWeight','bold');
gtext('Componentes Significantes','color','b');
gtext('Componentes no Significantes','color','r');

subplot(212);
ysig=m_Enfilada(:,end);
yerr=m_Enfilada(:,end)-pout_C2;
%nfft=389;
nfft=2^(nextpow2(length(m_Enfilada(:,end))));

bd=isnan(ysig);
gd=find(~bd);
bd([1:(min(gd)-1) (max(gd)+1):end])=0;
ysig(bd)=interp1(gd,ysig(gd),find(bd)); 
[Pxs,F]=psd(ysig(finite(ysig)),nfft,1,[],ceil(nfft/2));
%[Pxs,F]=pmtm(ysig(finite(ysig)),4,4096,1);
yerr(bd)=interp1(gd,yerr(gd),find(bd)); 
[Pxe,F]=psd(yerr(finite(ysig)),nfft,1,[],ceil(nfft/2));
%[Pxe,F]=pmtm(yerr(finite(ysig)),4,4096,1);

semilogy(F,Pxs);
set(gca,'Fontsize',12)
line(F,Pxe,'color','r');
xlabel('frequencIA (cph)','FontSize',12,'FontWeight','bold');
ylabel('m^2/cph','FontSize',12,'FontWeight','bold');
title('Espectro Estimado antes y después de remover la energía espectral','FontSize',12,'FontWeight','bold');
gtext('Serie Original (interpolada)','color','b');
gtext('Energía No-mareal analizada','color','r');

print('-dpng','-r400','armonic_enfilad2.png')