%***********************  SYNThetic SEASonal series  **********************
%      Builds synthetic seasonal series to be analyzed by diseas.m.       *
%      Program asks for periods to be incorporated to the series, until   *
%      an empty period is entered. For each period, gives option to vary  *
%      the phase linearly. This constitutes the "seasonal" signal.        *
%                                                                         *
%      After all periods have been entered, gives the options to add      *
%      linear trend and random normal noise.                              *
%                                                                         *
%      Optionally, saves the synthetic series, together with the          *
%      seasonal(s) as a .mat file. to be analyzed with diseas.m.          *
%                                                                         *
%                            Invoked as: syntseas                         *
%                                                                         *
%                                                        28/JUL/2023 Nava *
%**************************************************************************
function syntseas
%<<<<<<<<<<<<<<<<<<<<<<<<  Settable Parameters  <<<<<<<<<<<<<<<<<<<<<<<<<<<
clr='brkg';  % Colors for succesive periods in plot
t1def=0;     % Deafult initial time for series (yr)
Ndef=1000;   % Default number of elements in series
Adef=5.;     % Default amplitude in series (mm)
phidef=33.;  % Default initiial phase in series (deg)
swpl=1;      % Switch swpl=1 plots the series and its components
%<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
tpi=2*pi; 
%rtd=180/pi;
dtr=pi/180;

fprintf(':: SYNThetic SEASonal series ::           %s\n',datestr(now))

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)

fprintf('+ t1 ? [def=%.5f]: ',t1def)
t1=input('');
if isempty(t1), t1=t1def; end
    
fprintf('+ Number of elements ? [def=%i]: ',Ndef)
N=input('');
if isempty(N), N=Ndef; end
t2=t1+(N-1)*dt;
t=t1:dt:t2;  
tt1=t-t1;

if swpl
    fignum=input('+ Figure number ? [def=1/0=No Fig.]: ');
    if isempty(fignum), fignum=1; end
    if fignum==0
        swpl=0;
    else
        figure(fignum)
        clf
        coln=0;
        hdp=0;
        mins=1000; maxs=-1000;
    end
end

d=zeros(1,N);   %  d is the total synthetic signal
nT=0;      %  nT is number of seasonals, each with a different period
while 1   %-------------- For each period ---------------------------------
    T=input('+ Period T (year)? [def=stop]: ');
    if isempty(T), break, end
    nT=nT+1;
    fprintf('   + Amplitude (mm)? [def=%.3f mm]: ',Adef)
    A=input('');
    if isempty(A), A=Adef; end
    fprintf('   + Initial Phase (deg)? [def=%.3f deg]: ',phidef)
    phid=input('');
    if isempty(phid), phid=phidef; end
    phi=phid*dtr;
    
    %  phase is varied as initial phase plus phimd*sine from 0 to max (or 
    %        min if negative) over time T4 (period/4).
    hdd0=0;
    phimd=input('     + Vary phase by 0 to +/- (deg)? [def=0]: ');
    if isempty(phimd), phimd=0; end
    if phimd~=0, 
        dd0=A*cos(tt1*(tpi/T)+phi);
        figure(fignum), subplot(2,1,1),  hdd0=plot(t,dd0,'y');
        %phi=phi+tt1*(phimd*dtr); 
        phimd=phimd*dtr;
        T4=input('        + Over period/4 (yr)? [def=1]: ');
        if isempty(T4), T4=1; end
        Tmd=T4*4; fmd=1/Tmd;
        phi=phi+phimd*sin(tpi*fmd*tt1);
        %figure(75), clf, plot(tt1,phi*rtd,'b'), pause, figure(fignum)  % prova
    end
    
    dd=A*cos(tt1*(tpi/T)+phi);  %  dd is the seasonal component for each period
       
    swkeep=1;
    if swpl % If components are being plotted, gives option to keep or reject each
        figure(fignum)
        subplot(2,1,1)
        hold on
        coln=coln+1; if coln>4, coln=1; end
        hdd=plot(t,dd,clr(coln));     % plots each component
        ny=input('       + Keep [def=Y/N]: ', 's');

        if ~isempty(ny) && upper(ny)=='N' 
            delete(hdd)
            if hdd0~=0, delete(hdd0), end
            swkeep=0; 
            nT=nT-1;
        end
    end
    
    if swkeep
        s0{nT,1}=T; s0{nT,2}=dd;  % Seasonal component is saved for comparison
        
        if min(dd)<mins, mins=min(dd); end
        if max(dd)>maxs, maxs=max(dd); end
    
        d=d+dd;    % Adds each component to total seasonal

        if swpl
            subplot(2,1,2)
            hold on
            if hdp~=0, delete(hdp), end
            hdp=plot(t,d,'b');   % Plots total seasonal
        end
    end
    
end   % -------------------------------------------------------------------

%  Seasonals have been defined. Gives options for trend and noise__________
m=input('+ Add trend with slope (mm/yr) ? [def=0]: ');
if isempty(m), m=0; end
if m~=0
    tr=m*tt1;
    if swpl
        subplot(2,1,1)
        hold on
        coln=coln+1; if coln>4, coln=1; end
        lt=['--',clr(coln)];
        plot(t,tr,lt)   %  plots trend
        d=d+tr;         % Adds trend to signal
        subplot(2,1,2)
        hold on
        if hdp~=0, delete(hdp), end
        hdp=plot(t,d,'b');  % plots total signal
    end
end    

nsd=input('+ Add Noise with St.Dev= (mm) ? [def=No]: ');
if ~isempty(nsd) && nsd~=0
    nnu=input('   + Normal or Uniform ? [def=N/U]: ','s');
    if isempty(nnu) || upper(nnu)~='U', nnu='N'; else nnu='U'; end
    fprintf('     '), randseed(1);
    if nnu=='N'
        rn=randn(1,N);
    else
        rn=(rand(1,N)-.5)*2;
    end
    rn=rn*nsd/std(rn);
    d=d+rn;     %  Adds noise to signal
    if swpl
        subplot(2,1,1)
        hold on
        coln=coln+1; if coln>4, coln=1; end
        plot(t,rn,clr(coln))   %  plots noise

        subplot(2,1,2)
        hold on
        if hdp~=0, delete(hdp), end
        hdp=plot(t,d,'b');   % plots total signal
    end
end

if swpl   % Formats and labels plots_______________________________________
    subplot(2,1,1)
    axs=axis;
    dr=(maxs-mins)*.05; axs(3)=mins-dr; axs(4)=maxs+dr; axis(axs)
    set(gca,'FontSize',15)
    box on
    xlabel ('\fontsize{16}\itt  \fontsize{14}\rm(year)')
    ylabel ('\fontsize{16}\itd_{ }\rm(\itt_{ }\rm)  \fontsize{14}\rm(mm)')
    tab=axs(2)+.02*(axs(2)-axs(1));
    text(tab,axs(4),'\fontsize{15}(A)','horizontalalignment','left','verticalalignment','top')

    subplot(2,1,2)
    axs=axis;
    mins=min(d); maxs=max(d);
    dr=(maxs-mins)*.05; axs(3)=mins-dr; axs(4)=maxs+dr; axis(axs)
    set(gca,'FontSize',15)
    box on
    xlabel ('\fontsize{16}\itt  \fontsize{14}\rm(year)')
    ylabel ('\fontsize{16}\itd_{ }\rm(\itt_{ }\rm)  \fontsize{14}\rm(mm)')
    text(tab,axs(4),'\fontsize{15}(B)','horizontalalignment','left','verticalalignment','top')
end

% Output of processed series______________________________________________
filout=input('+ Save processed series as ? [def=No]: ','s'); 
if ~isempty(filout)
    t=t'; d=d';
    save(filout,'d','t', 's0')
    fprintf('* Saved >>>>>>>>>>>>>> %s.mat *\n',filout)
    if swpl
        subplot(2,1,2)
        title(['\fontsize{14}',filout])
    end
end
fprintf('..................... SYNT SEAS ............................\n\n')

%***********************  RANDom number generator SEED  *******************
%     seed is a scalar INTEGER value from  0 to 2^32-1 = 4.2950e+009      *
%     Invoked as randseed;  when invoked as randseed(1), if no seed is    * 
%     entered manually,prints randomly generated one.                     *
% 21/02/2014                                             17/May/2018 Nava *
%**************************************************************************
function randseed(varargin)
mlver=version; mlver=str2num(mlver(1));
seed=input('+ Seed? [def=Clock]: ');
if isempty(seed)
 	cl=clock;
    seed=rand*abs(cl(5)-cl(6))*10^round(6*rand);
    while seed>4.295e9, seed=seed/3; end
    seed=round(seed);
    if nargin>0 && varargin{1}==1, fprintf('       Seed = %i \n',seed), end
end % if isempty(sdst)
if mlver>7, rng('default'), rng(seed), else rand('twister',seed), end %{{{
