clear
data = load('SIMAR_OR_EXT_1044066.txt');
WaveData = data(:,1:11);


% Remove lines with 999 values
WaveData(WaveData == 999) = NaN;
Hs = WaveData(:,8);
ValidLines = find(~isnan(Hs));
WaveData = WaveData(ValidLines, :);


Tempo = datenum(WaveData(:,1),WaveData(:,2), WaveData(:,3), WaveData(:,4), 0, 0);
    
Hs = WaveData(:,8);
Tp = WaveData(:,10); %Peak period
Azimute0 = WaveData(:, 11);

StormThreshold = 5; % Wave height threshold
DeltaTMin = 0.5;    % Maximum time between consecutive storms


nt = 0;
fimTemporal = 1;
i = 1;
while i < length(Hs)
    if(Hs(i) > StormThreshold)
        if (Tempo(i) - Tempo(fimTemporal)) > DeltaTMin || nt == 0
            nt = nt + 1;
            it(nt) = i;
        end
        
        while Hs(i) > StormThreshold
            i = i + 1;
        end
        ft(nt) = i;
        fimTemporal = i;
    end
    i = i + 1;
end


for itemporal = 1:nt
    HsTemporal = Hs(it(itemporal):ft(itemporal)-1); %correcção introduzida para obviar o excesso de alturas em relação à duração da tempestade (-1)
    Azimute0Temporal = Azimute0(it(itemporal):ft(itemporal));
    AzimuteMedioTemporal(itemporal) = 180/pi*atan2(sum(cos(Azimute0Temporal*pi/180)), sum(sin(Azimute0Temporal*pi/180))); 
    Hsmax(itemporal) = max(HsTemporal);
    EnergiaTemporal(itemporal) = sum(HsTemporal.^2*3); % x3 horas, cada intervalo de medição
%     y = Tempo(it(itemporal));
%     Historico(itemporal,1) = str2num(datestr(y,10));    %ano do temporal
%     Historico(itemporal,2) = (Tempo(ft(itemporal)) - Tempo(it(itemporal)))*24; %duração em horas
%     Historico(itemporal,3) = Hsmax(itemporal); 
%     Historico(itemporal,4) = Historico(itemporal,3).^2*Historico(itemporal,2); % energia "Dolan"
%     Historico(itemporal,5) = EnergiaTemporal(itemporal); %energia da tempestade
%     Historico(itemporal,6) = AzimuteMedioTemporal(itemporal);
    %Period and storm beginning
    %TpTemporal = Tp(it(itemporal):ft(itemporal));
    %Historico(itemporal,1) = Tempo (it(itemporal)); 
    %Historico(itemporal,2) = max (TpTemporal); 
    %Historico(itemporal,3) = EnergiaTemporal(itemporal);
    
end
    
%save 047startandpeak.txt Historico -ascii   

%res = [Hsmax' AzimuteMedioTemporal' EnergiaTemporal' Tempo(it)-693960 Tempo(ft)-693960];

%indexTemporal = find(Hs > 5);
%DataTemporais = Tempo(indexTemporal);
%HsTemporais = Hs(indexTemporal);
%plot(DataTemporais, HsTemporais);
%datetick('x',12)