function [storms recovery]=selectStorms(varargin)
% addpath fullfile(fileparts(which(mfilename)),'..','InputData','Buoy')

OPT = struct(...
    'datapath', fullfile(fileparts(which(mfilename)),'..','InputData','Buoy'),...
    'netcdffilename', fullfile(fileparts(which(mfilename)),'..','InputData','Buoy', 'data.nc'),...
    'visualize', 0,...
    'modeltrain_number',5,...
    'threshold_factor',2,...
    'durationThreshold',36,...
    'location',4,...
    'method',1);

OPT = setproperty(OPT, varargin{:});

outputdir = fullfile(fullfile(fileparts(which(mfilename)),'..','InputData'),['Modeltrain_number_',num2str(OPT.modeltrain_number,'%02i')]);

addpath(genpath(fullfile(cd,OPT.datapath)));

%% from netCDF
if OPT.method==1
    waveheight = nc_varget(OPT.netcdffilename,'WVHT');
    waveheight(waveheight==nc_attget(OPT.netcdffilename, 'WVHT', 'FillValue'))=NaN;
    Time = nc_varget(OPT.netcdffilename,'time')+datenum(1970,1,1);
end
%% from .mat files
if OPT.method==2
    load('42007.mat');
    waveheight = D.stdmet.wave_height;
    Time = D.stdmet.time;
end


%% Threshold value,

% threshold
mean_WVHT = nanmean(waveheight);
threshold=mean_WVHT*OPT.threshold_factor;

% remove the values before a certain date 
time_interval = datenum(2002,1,1,0,0,0);
waveheight(Time<time_interval)=[]; 
Time(Time<time_interval)=[];
record_changed_WVHT=zeros(1,length(waveheight));

% replace the NaNs with the previous value, this makes it possible to look
% at a duration of the storm above a threshold
for i=1:length(waveheight)
    % replace the NaN
    if (abs(waveheight(i))>=0)<1
        try
        waveheight(i)= waveheight(i-1);
        catch
        waveheight(i)=-2;
        end
        record_changed_WVHT(i)=1;
    else
        record_changed_WVHT(i)=0;
    end
end

% crossings from under threshold/above threshold and v.v.
WVHT_above_threshold = waveheight>=threshold;
crossingthreshold=[0; diff(WVHT_above_threshold)];

pospoints=find(crossingthreshold==1);
negpoints=find(crossingthreshold==-1);
ii = min(length(pospoints),length(negpoints));
for i=1:ii
   duration(i)=(Time(negpoints(i))-Time(pospoints(i)));
end

duration_threshold = datenum(0,0,0,OPT.durationThreshold,0,0);
duration_above_threshold = duration>=duration_threshold;

Storms = [pospoints(find(duration_above_threshold)),negpoints(find(duration_above_threshold))];
for i=1:size(Storms,1)
    storms{i}.waveheight = waveheight(Storms(i,1):Storms(i,2));
    storms{i}.time = Time(Storms(i,1):Storms(i,2));
end

% get hurricane katrina
[katrina steepness Hm0]=getKatrina();

% Locate the maximum value of hurricane Katrina
[x y]=max(Hm0);
% time in seconds of the maximum afther the beginning
yt=y*katrina.input.dt;
yt=double(yt);
yt2=(length(Hm0)-y)*katrina.input.dt;
yt2=double(yt2);
yt3=katrina.input.dt;
yt3=double(yt3);

dt_start = datenum(0,0,0,0,0,yt);
dt_end = datenum(0,0,0,0,0,yt2);
dt_Katrina = datenum(0,0,0,0,0,yt3);


% find the maximums of the storms from the time serie
for i=1:size(storms,2)
    [xx(i) yy(i)]=max(storms{1,i}.waveheight);
    time_max_WVHT(i)= storms{1,i}.time(yy(i),1);
end

% all the start - points of the hurricanes Katrina
start_Katrinas = time_max_WVHT - dt_start;
end_Katrinas = time_max_WVHT + dt_end;

for i=1:size(storms,2)
    storms{i}.Katrinas.time = [start_Katrinas(i):dt_Katrina:(end_Katrinas(i)-dt_Katrina)];
    storms{i}.Katrinas.waveheight = Hm0;
end

[WL fieldname] = readTides('visualize',OPT.visualize,'location',OPT.location);

for istorm = 1:length(storms)
    [v1 v2] = find(abs(WL.(fieldname{OPT.location,1}{1,1}).time - storms{1,istorm}.Katrinas.time(1))<=0.05,1);
    [v3 v4] = find(abs(WL.(fieldname{OPT.location,1}{1,1}).time - storms{1,istorm}.Katrinas.time(end))<=0.05);
    if v1>0 & v3>0;
    storms{1,istorm}.surge.time = WL.(fieldname{OPT.location,1}{1,1}).time(v1:v3(end));
    storms{1,istorm}.surge.wl = WL.(fieldname{OPT.location,1}{1,1}).measurements(v1:v3(end));
    end
    [w1 w2] = max(storms{1,istorm}.surge.wl);
    [k1 k2] = max(katrina.surge.surge1);
    storms{istorm}.Katrinas.surge.one = (w1/k1)*katrina.surge.surge1(1:end);
    [w1 w2] = max(storms{1,istorm}.surge.wl);
    [k1 k2] = max(katrina.surge.surge2);
    storms{istorm}.Katrinas.surge.two = (w1/k1)*katrina.surge.surge2(1:end);
    [w1 w2] = max(storms{1,istorm}.surge.wl);
    [k1 k2] = max(katrina.surge.surge3);
    storms{istorm}.Katrinas.surge.three = (w1/k1)*katrina.surge.surge3(1:end);
    [w1 w2] = max(storms{1,istorm}.surge.wl);
    [k1 k2] = max(katrina.surge.surge4);
    storms{istorm}.Katrinas.surge.four = (w1/k1)*katrina.surge.surge4(1:end);
end


for i=1:size(storms,2)-1
    recovery.time(i) = storms{i+1}.Katrinas.time(1,1) - storms{i}.Katrinas.time(1,end);
end
recovery_small=((recovery.time)/30)<1;
for i=1:size(storms,2)
    max_wvht_Storms(i)=max(storms{1,i}.waveheight);
end
for i=1:length(recovery_small)
    if recovery_small(i)>0
        if max_wvht_Storms(i)>=max_wvht_Storms(i+1)
            max_wvht_Storms(i+1)=0;
        else
            max_wvht_Storms(i)=0;
        end
    end
end
storms = storms(max_wvht_Storms>0);
clear recovery;
for i=1:size(storms,2)-1
    recovery.time(i) = storms{i+1}.Katrinas.time(1,1) - storms{i}.Katrinas.time(1,end);
    t2 = storms{i+1}.Katrinas.time(1,1);
end
i=size(storms,2);
recovery.time(i) = datenum(now) - t2;

for i=1:length(storms)
    scaling(i)=(max(Hm0)/max(storms{1,i}.waveheight));
    [regime(i)]=scaleKatrina(katrina,i,scaling(i),steepness,Hm0,OPT.modeltrain_number);
end

threshold_waveheight =['threshold waveheight = ', [num2str(mean_WVHT),'(mean WVHT)'],'*',num2str(OPT.threshold_factor),'  [m]'];
threshold_duration =['threshold duration = ',num2str(duration_threshold), '        days'];
number_of_storms = ['number of storms = ',num2str(size(storms,2))];

fid = fopen(fullfile(outputdir,'readme.txt'),'wt');
fprintf(fid,'%s\n','This readme-file is automatically generated');
fprintf(fid,'%s\n','The details of the modeltrain are listed below');
fprintf(fid,'%s\n','');
fprintf(fid,'%s\n',threshold_waveheight);
fprintf(fid,'%s\n',threshold_duration);
fprintf(fid,'%s\n',number_of_storms);
fclose(fid);

for i=1:length(regime)
    WVHTS_storms(:,i) = regime(i).waves';
end

for i = 1:length(storms)
    if abs(storms{i}.Katrinas.surge.one(1,1))>=0
        k(i)=0;
        l(i)=i;
    else
        k(i)=1;
        l(i)=i;
    end
end
for i = 1:length(storms)
    max_storms(i) = max(storms{i}.waveheight);
end

Storms_wo_surge = l(k>0);
for i = 1:length(Storms_wo_surge)
    % find the storm with simular waveheight
    stormHeight = max_storms(Storms_wo_surge(i));
    max_storms(i) = Storms_wo_surge(i)+10;
    [wa1 wa2] = min(abs(max_storms(Storms_wo_surge(end)+1:end)-stormHeight));
    storm_nummertje(i)=Storms_wo_surge(i);
    wordt_storm(i)=(wa2+Storms_wo_surge(end)-1);
    storms{Storms_wo_surge(i)}.Katrinas.surge = storms{(wa2+Storms_wo_surge(end)-1)}.Katrinas.surge;
end
for i = 1:length(storms)
    max_storms(i) = max(storms{i}.waveheight);
end

for i = 1:length(storms)
    fid = fopen(fullfile(outputdir,['Storm_',num2str(i,'%02i')],['Scaled_Surge_Storm_',num2str(i,'%02i'),'.dat']),'w');
    for j=1:length(katrina.input.files.list)+1
        fprintf(fid,'%6.2f %6.2f %6.2f %6.2f %6.2f\n',katrina.surge.time(j),storms{i}.Katrinas.surge.one(j),storms{i}.Katrinas.surge.two(j),storms{i}.Katrinas.surge.three(j),storms{i}.Katrinas.surge.four(j));
    end
end

dlmwrite(fullfile(outputdir,'WVHTS_storms.txt'),WVHTS_storms, 'delimiter', '\t');

save('storms','storms')
save('recovery','recovery')
movefile('storms.mat',outputdir)
movefile('recovery.mat',outputdir)

%% visualize 

if OPT.visualize==1
    figure('Name','Selected Storms');
    plot(Time,waveheight);
    hold on;
    plot(Time(record_changed_WVHT>0),waveheight(record_changed_WVHT>0),'.g','MarkerSize',1)
    for i=1:size(storms,2)
        plot(storms{i}.time,storms{i}.waveheight, 'r','LineWidth',2),hold on
        datetick('x')
        plot(storms{i}.Katrinas.time,regime(i).waves(1,:),'--k', 'LineWidth',2)
        plot(storms{i}.surge.time,storms{i}.surge.wl,'m')
        plot(storms{i}.Katrinas.time,storms{i}.Katrinas.surge.one(2:end),'g')
        text(storms{i}.Katrinas.time(1),(max(storms{i}.waveheight)+0.2),num2str(i),'FontSize',14,'FontWeight','bold'),hold on
    end
end