%**********************  DIsplacement SEASonal analysis  ******************
%     Fits a given seasonal period to a displacement series d by one of   *
%     three methods,  then subtracts the seasonal fit from the original   *                            *
%     series.                                                             *
%         Longer periods should be identified and eliminated first!       *
%                                                                         *
%     Optionally, saves the filtered displacement                         *
%                                                                         *
%     Invoked as: diseas filename, where filename is the name of a .mat   *
%     file containing a displacements column vector named d.              *
%     The file may contain a column vector t with the times corresponding *
%     to the displacements; if no vector t is found, the program asks for *
%     the time corresponding to the first element of d.                   *
%                                                                         *
%                                                                         *
%                                                        28/JUL/2023 Nava *
%**************************************************************************
function diseas(filin)
%<<<<<<<<<<<<<<<<<<<<<<<<  Settable Parameters  <<<<<<<<<<<<<<<<<<<<<<<<<<<
swcut=1;    % Switch swcut=1 gives the option to cut the series' length
swpr=0;     % Switch swpr=1 prints intermediate results
swap=1;     % Switch swap=1 plots instantaneous amplitude and phase for variable phase
%<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
tpi=2*pi; 
rtd=180/pi;

fprintf(':: Displacement Seasonal Analysis ::           %s\n',datestr(now))

filin=upper(filin);
load(filin)
tit0=['\fontsize{14}',filin];
d=d';
N=length(d);  

dt=input('+ dt (yr)? [def=1/365 year]: ');
if isempty(dt), dt=1/365; end
dtd=dt*365;
fprintf('::  dt = %.5fyr = %.5fday ::\n',dt,dtd)

if exist('t')
    t=t'; 
else
    t0=input('+ t0 (yr)? [def=0.0 yr]: ');
    if isempty(t0), t0=0; end
    t=[0:(N-1)]*dt+t0;
end

tt=N*dt;
fprintf(':: Series with %i elements and dt=%fyear spans %f years (%f, %f) ::\n',...
                                                         N,dt,tt,t(1),t(N))
                                                     
if swcut   % Option to cut the series _____________________________________
    nk=input('   + Keep ? elements [def=all]: ');
    if ~isempty(nk)
        i1=input('    + starting at index ? [def=1]: ');
        if isempty(i1), i1=1; end
        i2=i1+nk-1;
        N=nk;
        t=t(i1:i2); d=d(i1:i2); tt=N*dt;
        fprintf(':: Series with %i elements and dt=%f spans %f years (%f, %f) ::\n',...
                                                         N,dt,tt,t(1),t(N))
    end
end
    
fignum=input('+ 1st. Figure number? [def=1]: ');
if isempty(fignum), fignum=1; end
fignum1=fignum;
figure(fignum1)
clf
subplot(2,1,1)
hold on
plot(t,d,'b')
tr=(t(N)-t(1))*.02; t1=t(1)-tr; t2=t(N)+tr;
dmn=min(d); dmx=max(d); dr=(dmx-dmn)*.05;
axs=[t1 t2 dmn-dr dmx+dr];
axis(axs)
set(gca,'FontSize',15)
box on
xlabel ('\fontsize{16}\itt  \fontsize{14}\rm(year)')
ylabel ('\fontsize{16}\itd_{ }\rm(\itt_{ }\rm), \its_{ }\rm(\itt_{ }\rm)  \fontsize{14}\rm(mm)')
title(tit0)
tab=axs(2)+.2*tr;
text(tab,axs(4),'\fontsize{15}(A)','horizontalalignment','left','verticalalignment','top')

Ts=input('+ Seasonal Period T ? [def=1yr]: ');  % Ts is nominal period
if isempty(Ts), Ts=1; end
NT=round(Ts/dt);
T=NT*dt;          % T is the actual period locatable with the given dt
fo=1/T;
tpif=tpi*fo;
fprintf('::  Period T=%fyr=%.3fday  fo=%f yr^-1  %i elements ::\n',T,T*365,fo,NT)
nqhm=0;
if tt>=T   %================ tt>T =========================================
    tpit=(tpif)*[0:NT-1]*dt;  
    ct=cos(tpit);
    st=sin(tpit);

    df=1/T; 
    iq=2;

    NTm=NT-1;
    m=N-NTm;
    ph=zeros(1,m);  a=zeros(1,m); 

    for j=1:m
        dj=d(j:j+NTm);
        
        Dre=sum(ct.*dj);  
        Dim=sum(st.*dj); 
        D1=Dre-i*Dim;
        phi=phase1(D1);
        a(j)=sqrt(Dre^2+Dim^2)*2/NT; 
        ph(j)=phi-tpif*(t(j)-t(1));
        while ph(j)>tpi, ph(j)=ph(j)-tpi; end
        while ph(j)<-tpi, ph(j)=ph(j)+tpi; end
        if swpr
            fprintf('i=%i i1s=%i D1=[%e,%ei] A=%f \n',j,i1s(j),real(D1),imag(D1),a(j))
            fprintf('     phi=%f (%fdeg) phase=%f (%fdeg)\n',phi,phi*r2d,ph(j),ph(j)*rtd)
            pause
        end
    end
    if m>1
        pha=ph(1);
      	for j=2:m
            dph=ph(j)-pha;
            if dph>pi
                ph(j)=ph(j)-tpi;
            elseif dph<-pi
                ph(j)=ph(j)+tpi;
            end
            pha=ph(j);
        end
    end

    am=mean(a);  asd=std(a); phm=mean(ph); phmd=phm*rtd; phsd=std(ph)*rtd;
    fprintf('*** T=%fyr: <A>=%fmm (StDA=%fmm), <Ph>=%fdeg (StDPh=%fdeg) ***\n',...
                    T,am,asd,phmd,phsd)   

    nqh=round(NT/2);  nqhm=nqh-1;
    ix=[1:m]+nqhm;
    ix1=nqh; ix2=N-nqhm;
    if swap
        figure(fignum+1)
        clf
        subplot(2,1,1)
        hold on
        ph1=min(ph)*rtd; ph2=max(ph)*rtd; 
        phr=(ph2-ph1); if phr<.00175, phr=.00175; end  % ~2deg
        phr=phr*.075;
        ph1=ph1-phr; ph2=ph2+phr;
        plot(t(ix1:ix2),ph*rtd,'b')
        plot([t(ix1) t(ix2)],[phmd phmd],'r')
        plot([t(ix1) t(ix2)],[phmd-phsd phmd-phsd],'--r')
        plot([t(ix1) t(ix2)],[phmd+phsd phmd+phsd],'--r')
        axis([axs(1) axs(2) ph1 ph2]) 
        box on
        grid on
        set(gca,'FontSize',15)
        xlabel ('\fontsize{16}\itt  \fontsize{14}\rm(year)')
        ylabel ('\fontsize{16}\it\phi   \fontsize{14}\rm(\circ)')
        text(tab,ph2,'\fontsize{15}(A)','horizontalalignment','left','verticalalignment','top')
        title(tit0)

        subplot(2,1,2)
        hold on
        a1=min(a); a2=max(a); ar=(a2-a1); if ar<.01, ar=.01; end
        ar=ar*.075;
        a1=a1-ar; a2=a2+ar;
        plot(t(ix1:ix2),a,'b')
        plot([t(ix1) t(ix2)],[am am],'r')
        plot([t(ix1) t(ix2)],[am-asd am-asd],'--r')
        plot([t(ix1) t(ix2)],[am+asd am+asd],'--r')
        axis([axs(1) axs(2) a1 a2]) 
        box on
        grid on
        set(gca,'FontSize',15)
        xlabel ('\fontsize{16}\itt  \fontsize{14}\rm(year)')
        ylabel ('\fontsize{16}\itA   \fontsize{14}\rm(mm-yr)')
        text(tab,a2,'\fontsize{15}(B)','horizontalalignment','left','verticalalignment','top')
    end
    
    ny=input(' Seasonal Variable/Average Phi [def=V/A]: ','s');
    if ~isempty(ny) && upper(ny)=='A'     %   Average seasonal-------------
         s=am*cos((t-t(1))*tpif+phm);
         swvar=0;
    else                                % Variable phi seasonal -----------
        if m>1
            phb=ones(1,nqhm)*ph(1);   
            phe=ones(1,nqhm)*ph(m);
            ph=[phb,ph,phe];  
            ab=ones(1,nqhm)*a(1);   
            ae=ones(1,nqhm)*a(m);
            a=[ab,a,ae]; 
        end    
        s=a.*cos((t-t(1))*tpif+ph);
        swvar=1;
    end                 %  Var phi seasonal -------------------------------
else   % ============== tt<T  =============================================
    fprintf('## Time=%f < Period=%f! -> seafit ##\n',tt,T) 
    [s,ph,am,m]=seafit(d,t(1),dt,T);
     phd=ph*rtd;
     fprintf('*** N=%i: A=%fmm, Ph=%fdeg  %i iterations ***\n',...
                                                 length(s),am,ph*rtd,m)   
     swvar=0; 
end  %=====================================================================

if exist('s0')  % If series was synthetic created by syntseas.m, s0 contains true seasonals
    [nts,~]=size(s0);
    for j=1:nts, Ts0(j)=s0{j,1}; end,   % Ts0 are nominal seasonal periods
    ix=find(Ts0==Ts); 
    if ~isempty(ix)
        s0T=s0{ix,2};
        ds=s-s0T(1:N);
        erms=sqrt(mean(ds.^2));
        if ~swvar
            fprintf('*** RMS err Identified-True = %fmm ***\n',erms)
        else
            ermst=sqrt(mean(ds(ix1:ix2).^2));
            fprintf('*** RMS err Identified-True = %fmm w/o tails = %fmm ***\n',erms,ermst)
        end
        figure(fignum), subplot(2,1,1)
        hold on,
        plot(t,s0T(1:N),':k','linewidth',2)
        s0(ix,:)=[];        %  Deletes synthetic component
        [nts,~]=size(s0);
    end
end

d=d-s;    %  Subtracts stational from displacement
figure(fignum), subplot(2,1,1)
hold on,
axs=axis;
smn=min(s); smx=max(s); sr=(smx-smn)*.05; smn=smn-sr; smx=smx+sr;
if smn<axs(3), axs(3)=smn; end
if smx>axs(4), axs(4)=smx; end
axis(axs)
plot(t,s,'r','linewidth',1)
if swvar
    plot([t(ix1) t(ix1)],[axs(3) axs(4)],':r')
    plot([t(ix2) t(ix2)],[axs(3) axs(4)],':r')
end

subplot(2,1,2)
plot(t,d,'k','linewidth',1)
dmn=min(d); dmx=max(d); 
if dmn<axs(3) || dmx>axs(4), dr=(dmx-dmn)*.05; axs=[t1 t2 dmn-dr dmx+dr]; end
axis(axs)
if swvar
    hold on
    plot([t(ix1) t(ix1)],[axs(3) axs(4)],':r')
    plot([t(ix2) t(ix2)],[axs(3) axs(4)],':r')
end

set(gca,'FontSize',15)
box on
xlabel ('\fontsize{16}\itt  \fontsize{14}\rm(year)')
ylabel ('\fontsize{16}\itd_{ }\rm(\itt_{ }\rm) - \its_{ }\rm(\itt_{ }\rm)  \fontsize{14}\rm(mm)')
    
title([tit0,'-S'])
text(tab,axs(4),'\fontsize{15}(B)','horizontalalignment','left','verticalalignment','top')

% Output of processed series______________________________________________
nwf=0;   % Switch nwf=1 means the name of the output series is changed
filout=[filin,'-S'];
fprintf('+ Save processed series as ? [def=%s / N]: ',filout)
cat2=input('','s'); 
if ~isempty(cat2), filout=cat2; nwf=1; end
if length(filout)~=1 || upper(filout(1))~='N' 
    t=t'; d=d';
    if exist('s0') && nts>0
        save(filout,'d','t','s0')
    else
        save(filout,'d','t')
    end
    fprintf('* Saved >>>>>>>>>>>>>> %s.mat *\n',filout)
    if nwf
        subplot(2,1,2)
        title(['\fontsize{14}',filout])
    end
end

fprintf('..................... Diseas ............................\n\n')

%****************************  SEAsonal FIT  ******************************
%	 Bootstrap scheme for fitting seasonals with period T from segments   *
%    shorter than T.                                                      *
%	                                                                      *
%                                                                         *
%                                                        22/Jul/2023 Nava *
%**************************************************************************
function [s,ph,aa,nit]=seafit(d,t1,dt,T)
%<<<<<<<<<<<<<<<<<<<<<<<<  Settable Parameters  <<<<<<<<<<<<<<<<<<<<<<<<<<<
tpi=2*pi; 
rtd=180/pi;
crit=.001;  % Convergence criterion: iteration stops when succesive results differ by <crit
nitx=5000;  % Max no. of iterations
swpl=1;     % Switch swpl=1 plots the fits onto the original series
swpr=0;     % Switch swpr=1 prints intermediate output
%<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
NT=round(T/dt);
t=[0:NT-1]*dt; 
tpit=(tpi/T)*t;  
ct=cos(tpit);
st=sin(tpit);

N=length(d);
if swpr, fprintf(':: seafit: N=%i, t1=%f ::\n',N,t1), end
if swpl, tt1=t+t1; end
n=N;
pha=10; aaa=-1; nit=0;
while 1
    Dre=sum(ct(1:n).*d); 
    Dim=sum(st(1:n).*d); 
    if swpr, fprintf('Dre=%f Dim=%f\n',Dre,Dim), end   % prova
    D1=Dre-i*Dim;
    ph=phase1(D1);
    aa=sqrt(Dre^2+Dim^2)*2/n; 
    nit=nit+1;
    if swpr, fprintf('*** nit=%i  N=%i: Measured Ph=%f, A=%f  ***\n',nit,n,ph*rtd,aa), end
    s=aa*cos(tpit(1:NT)+ph);  
    if swpl, plot(tt1(1:n),s(1:n),'y'), end  
    
    dph=abs(ph-pha)*rtd; daa=abs(aa-aaa); 
    if dph<crit && daa<crit, break, end
    if nit>=nitx, break, end
    n=NT;
    d=[d(1:N),s(N+1:n)];
    pha=ph; aaa=aa; 
end
s=s(1:N);
if swpl, 
    axs=axis; 
    if min(s)<axs(3), axs(3)=min(s); end
    if max(s)>axs(4), axs(4)=max(s); end
    axis(axs)
end

%*****************************  PHASE 1  **********************************
%              Finds phase [0,2*pi] for a complex number c.               *
%                                                                         *
%                                                        21/Dec/2021 Nava *
%**************************************************************************
function ph=phase1(c)
pih=pi/2; pi2=2*pi; 
rec=real(c); imc=imag(c);
if abs(rec)<2.5e-15, rec=0; end
if abs(imc)<2.5e-15, imc=0; end

if rec==0
    ph=pih; 
    if imc<0, ph=pih+pi; elseif imc==0, ph=0; end
else
    ph=atan(imc/rec);
    if rec<0, ph=ph+pi; end 
end
if ph<0, ph=ph+pi2; end


