%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% - Legge i dati da un file testo di dati ondametrici che ha il formato
%   aaaa	MM	dd	hh	mm	SWH(m)  MWD(°N)	Tm(s)	Tp(s)	Temp(°C)
% - dati mancanti: NaN
% - identifica i dati con altezza superiore ad un valore
% - identifica inizio e fine di ogni blocco di dati con hs> soglia
% - identifica inizio e fine di mareggiate con determinate caratteristiche
%   di altezza e durata.% - Per ogni mareggiata si calcola la relativa energia e la classifica
% - per ogni mareggiata calcola i valori medi delle principali grandezze e
%   calcola i massimi (ed i corrispondenti valori delle altre variabili in
%   concomitanza del massimo trovato)
% 
% 20090408 Andrea Valentini avalentini@arpa.emr.it
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

fname='Nausicaa_070601_081231.txt';
dati=load(fname,'-ascii');
name='Nau_070601_081231_av';

aa=dati(:,1);       % anno
MM=dati(:,2);       % mese
dd=dati(:,3);       % giorno
hh=dati(:,4);       % ore
mm=dati(:,5);       % minuti
hs=dati(:,6);       % Alt. (sign.) m
dir=dati(:,7);      % Wave direction °N
tm=dati(:,8);       % Per. medio m
tp=dati(:,9);       % Per. picco m
temp=dati(:,10);    % Temp sup acqua °C

null=find(hs==0);
hs(null)=NaN;tm(null)=NaN;tp(null)=NaN;dir(null)=NaN;
null=find(tm==0);
tm(null)=NaN;
null=find(tp==0);
tp(null)=NaN;
null=find(dir==0);
dir(null)=NaN;


% Imposto lo step temporale di acquisizione dati
% MOLTO IMPORTANTE per il calcolo energetico
% 3    dati triorari
% 1    dati orari
% 0.5  dati semiorari
% 0.25 ogni 15 minuti
hstep=0.5;

jd=julian([aa MM dd hh mm mm*.00])-2440000;

nt=length(jd);
ind=(1:nt);
mask=zeros(nt,1); %preparo una maschera delle dimensioni dei dati

% seleziono gli eventi sopra una certa soglia
ind1=find(hs>=1.5);

% scrivo la data di tutti i dati che hanno altezza significativa sopra la
% soglia
datetot=gregorian(jd(ind1)+2440000);
filename=strcat(num2str(name),'_datetot.mat');
save (filename, 'datetot');
filename=strcat(num2str(name),'_datetot.txt');
save (filename, 'datetot', '-ascii');

mask(ind1)=1; % metto a 1 tutti i dati con hs>=1.5

% trova differenze superiori a 3h per catalogare mareggiate diverse
% in questo modo ammetto che ci sia un buco di 3h all'interno dello stesso
% blocco, buco dovuto o a dati mancanti o a valori inferiori alla soglia
diffe=find(diff(ind1)>6);

% gestisco il fatto che se l'ultimo blocco ha durata superiore a 6h
% devo fare in modo di inserire l'indice dell'estremo finale di ind1
% che è uguale a length(ind1) in diffe
if ind1(end)-ind1(diffe(end))>12
% av_20090407 fino ad oggi c'era 6 anziché 12, ma dalla documentazione
% sembra che debba essere il tempo di durata di una mareggiata (6h) e non il gap
% tra una mareggiata e l'altre (3H)
    diffe=cat(1,diffe,length(ind1));
end

date=ones(length(diffe)*2,1).*NaN;

% questo blocco isola le date di inizio e di fine di ogni gruppo di dati
% attigui che hanno hs>=1.5m
for i=1:length(diffe)-1;
    date(1)=jd(ind1(1)); %inizio prima
    date(2*i)=jd(ind1(diffe(i))); %fine prima e fine altre
    date(2*i+1)=jd(ind1(diffe(i)+1)); %inizio altre     
    date(2*length(diffe))=jd(ind1(end)); %fine ultima
end

% salvo queste date
dateg=gregorian(date+2440000);
filename=strcat(num2str(name),'_date.mat');
save (filename, 'dateg');
filename=strcat(num2str(name),'_date.txt');
save (filename, 'dateg', '-ascii')


% trova diff sup a 6h per catalogare mareggiate diverse che abbiano una
% durata di almeno 6h.
dur=find(diff(diffe)>12);
date2=ones(length(dur)*2,1).*NaN;
% qui gestisco il primo blocco in quanto diff fa la differenza tra indici
% consecutivi e quindi il primo risultato è la differenza tra secondo
% blocco e primo, non considerando la lunghezza del primo
if (ind1(diffe(1))-ind1(1)>=12); 
    for i=1:length(dur);
        date2(1)=jd(ind1(1)); %inizio prima
        date2(2)=jd(ind1(diffe(1))); %fine prima
        date2(2*i+1)=jd(ind1(diffe(dur(i))+1)); %inizio altre     
        date2(2*i+2)=jd(ind1(diffe(dur(i)+1))); %fine prima e fine altre
    end
else
    for i=1:length(dur);
        date2(2*i-1)=jd(ind1(diffe(dur(i))+1));%inizio
        date2(2*i)=jd(ind1(diffe(dur(i)+1)));%fine
    end
end
date2g=gregorian(date2+2440000);
filename=strcat(num2str(name),'_date_mareggiate.mat');
save (filename, 'date2g');
filename=strcat(num2str(name),'_date_mareggiate.txt');
save (filename, 'date2g', '-ascii');

% converto in data giuliana le date delle mareggiate
jd_mareg=julian(date2g)-2440000;
%potevo anche considerare direttamente date2

% Assignment of values for the Wave energy classification
% according to Mendoza & Jimenez, 2004
I=58.4;     %weak...energia in (m^2*hr)
II=127.9;   %moderate
III=389.7;  %significant
IV=706.9;   %severe
%V > 706.9   %extreme


%inizializzo gli array
id_evento=ones((length(date2g)/2),1).*NaN;
duration=ones((length(date2g)/2),1).*NaN;
energy=ones((length(date2g)/2),1).*NaN;
class=ones((length(date2g)/2),1).*NaN;

hs_mean=ones((length(date2g)/2),1).*NaN;
tm_mean=ones((length(date2g)/2),1).*NaN;
tp_mean=ones((length(date2g)/2),1).*NaN;
dir_mean=ones((length(date2g)/2),1).*NaN;

dateg_hs_Max=ones((length(date2g)/2),6).*NaN;
hs_Max=ones((length(date2g)/2),1).*NaN;
tm_hs_Max=ones((length(date2g)/2),1).*NaN;
tp_hs_Max=ones((length(date2g)/2),1).*NaN;
dir_hs_Max=ones((length(date2g)/2),1).*NaN;

dateg_tm_Max=ones((length(date2g)/2),6).*NaN;
tm_Max=ones((length(date2g)/2),1).*NaN;
hs_tm_Max=ones((length(date2g)/2),1).*NaN;
tp_tm_Max=ones((length(date2g)/2),1).*NaN;
dir_tm_Max=ones((length(date2g)/2),1).*NaN;

dateg_tp_Max=ones((length(date2g)/2),6).*NaN;
tp_Max=ones((length(date2g)/2),1).*NaN;
hs_tp_Max=ones((length(date2g)/2),1).*NaN;
tm_tp_Max=ones((length(date2g)/2),1).*NaN;
dir_tp_Max=ones((length(date2g)/2),1).*NaN;

tabella=ones((length(date2g)/2),51).*NaN;
tab_ciav=ones((length(date2g)/2),26).*NaN;


id=1;
for i=1:(length(date2g)/2);
    
%    id=strcat(num2str(i),'_',num2str(name));    %event index
%     id_evento(i)=i;

    if i==1
        id_evento(i)=id;
        id=id+1;
    else if date2g(2*i-1,1)==date2g(2*i-3,1)
        id_evento(i)=id;
        id=id+1;
    else
        id_evento(i)=1;
        id=2;
        end
    end
    
   
%    datag_start(i)=date2g(2*i-1,:); %event start date
    duration(i)=(etime(date2g(2*i,:),date2g(2*i-1,:)))/3600; %event duration in h
    
    ind_event=find(jd>=jd_mareg(2*i-1) & jd<=jd_mareg(2*i));
    
    hs1=hs(ind_event);
    energy(i)=nansum(hstep*((hs1.*hs1))); %energy of the event
    %NB!! if hourly measurements not 0.5 but 1 (maybe... control!)
    
    %classification of the storm according to Mendoza & Jimenez, 2004
    if I>energy(i)
        class(i)=1;
    else if II>energy(i)
        class(i)=2;
    else if III>energy(i)
        class(i)=3;
    else if IV>energy(i)
        class(i)=4;
    else
        class(i)=5;
    end
    end
    end
    end

      
    % mean values for the variables
    hs_mean(i)=nanmean(hs(ind_event)); %hs mean in the period
    tm_mean(i)=nanmean(tm(ind_event)); %mean of the mean period during the storm period
    tp_mean(i)=nanmean(tp(ind_event)); %mean of the level during the storm period
    dir_mean(i)=pmean(dir(ind_event)*pi/180)*180/pi; %mean of the wave direction

    % search of the maximum hs in the period and the respective variables
    ind_hs_Max=find(hs(ind_event)==max(hs(ind_event)));
    if (isempty(ind_hs_Max))
        hs_Max(i)=NaN;
        tm_hs_Max(i)=NaN;
        tp_hs_Max(i)=NaN;
        dir_hs_Max(i)=NaN;
        date_hs_Max(i)=NaN;
        dateg_hs_Max(i,:)=NaN;
    else
        hs_Max(i)=hs(ind_event(ind_hs_Max(1)));
        tm_hs_Max(i)=tm(ind_event(ind_hs_Max(1)));
        tp_hs_Max(i)=tp(ind_event(ind_hs_Max(1)));
        dir_hs_Max(i)=dir(ind_event(ind_hs_Max(1)));
        date_hs_Max(i)=jd(ind_event(ind_hs_Max(1)));
        dateg_hs_Max(i,:)=gregorian(date_hs_Max(i)+2440000);
    end
    
    % search of the maximum tm in the period
    ind_tm_Max=find(tm(ind_event)==max(tm(ind_event)));
    if (isempty(ind_tm_Max))
        tm_Max(i)=NaN;
        hs_tm_Max(i)=NaN;
        tp_tm_Max(i)=NaN;
        dir_tm_Max(i)=NaN;
        date_tm_Max(i)=NaN;
        dateg_tm_Max(i,:)=NaN;
    else
        tm_Max(i)=tm(ind_event(ind_tm_Max(1)));
        hs_tm_Max(i)=hs(ind_event(ind_tm_Max(1)));
        tp_tm_Max(i)=tp(ind_event(ind_tm_Max(1)));
        dir_tm_Max(i)=dir(ind_event(ind_tm_Max(1)));
        date_tm_Max(i)=jd(ind_event(ind_tm_Max(1)));
        dateg_tm_Max(i,:)=gregorian(date_tm_Max(i)+2440000);
    end

    % search of the maximum tp in the period
    ind_tp_Max=find(tp(ind_event)==max(tp(ind_event)));
    if (isempty(ind_tp_Max))
        tp_Max(i)=NaN;
        hs_tp_Max(i)=NaN;
        tm_tp_Max(i)=NaN;
        dir_tp_Max(i)=NaN;
        date_tp_Max(i)=NaN;
        dateg_tp_Max(i,:)=NaN;
    else
        tp_Max(i)=tp(ind_event(ind_tp_Max(1)));
        hs_tp_Max(i)=hs(ind_event(ind_tp_Max(1)));
        tm_tp_Max(i)=tm(ind_event(ind_tp_Max(1)));
        dir_tp_Max(i)=dir(ind_event(ind_tp_Max(1)));
        date_tp_Max(i)=jd(ind_event(ind_tp_Max(1)));
        dateg_tp_Max(i,:)=gregorian(date_tp_Max(i)+2440000);
    end

            
    tabella(i,:)=[id_evento(i) date2g(2*i-1,1) date2g(2*i-1,:) date2g(2*i,:) duration(i) energy(i) class(i) hs_mean(i) tm_mean(i) tp_mean(i) dir_mean(i) ...
        dateg_hs_Max(i,:) hs_Max(i) tm_hs_Max(i) tp_hs_Max(i) dir_hs_Max(i) ...
        dateg_tm_Max(i,:) tm_Max(i) hs_tm_Max(i) tp_tm_Max(i) dir_tm_Max(i) ...
        dateg_tp_Max(i,:) tp_Max(i) hs_tp_Max(i) tm_tp_Max(i) dir_tp_Max(i)];
    
    tab_ciav(i,:)=[id_evento(i) date2g(2*i-1,1) date2g(2*i-1,:) energy(i) class(i) duration(i) dir_hs_Max(i) ...
        hs_Max(i) hs_mean(i) dateg_hs_Max(i,:) ...
        tp_Max(i) tm_Max(i) tp_mean(i) tm_mean(i) tp_hs_Max(i) tm_hs_Max(i)];
    

end

filename=strcat(num2str(name),'_tabella.mat');
save (filename, 'tabella');
filename=strcat(num2str(name),'_tabella.txt');
save(filename, 'tabella', '-ascii');

filename=strcat(num2str(name),'_tab_ciav.mat');
save (filename, 'tab_ciav');
filename=strcat(num2str(name),'_tab_ciav.txt');
save(filename, 'tab_ciav', '-ascii');

