%
%  This script develops a windfield that includes secondary eyewalls (if
%  present) and plots the evolution of the secondary wind speed maximimum.  
%
%  Copyright WindRiskTech, LLC, 2014
%  Last modified March, 2020
%-------------------------------------------------------------------------
%
months={'January','February','March','April','May','June','July','August','September','October','November','December'};
days=[31,28,31,30,31,30,31,31,30,31,30,31];
params
try
    load temp vsestore
catch
    disp('No secondary eyewall information for this event set, quitting')
    return
end    
load temp vnet vnetse hourstore daystore monthstore shape city_radius 
%if exist('vsestore','var') == 0
%    load('temp','vnet','vnetse','hourstore','daystore','monthstore','shape','city_radius')
%end    
load sorted
clf
%
%  Prompt, if necessary, for storm rank number
%
if ~exist('nt','var')
   nt=get_storm_num;
end   
%
vnetse1(1,:)=vnetse(nt,:);
[ni,nj]=size(vnet);
time=daystore(nt,1)+(0.5*hourstore(nt,1)+(0:nj-1))/12;
%
clear dir
FileInfo = dir('stats.txt');
TimeStamp = FileInfo.date;
%
if datenum(TimeStamp) < 7.3685e+05-1
    vold=1;
    m=0;
    vcon=zeros(nj,2);
    jcon=zeros(1,nj);
    for j=1:nj
        if vnetse(nt,j) == 0
            vnetse1(1,j)=NaN;
        end
        if j > 2 && vnetse(nt,j-1) > 0 && vnetse(nt,j-2) > 0 && ...
                abs(vnetse(nt,j)-vnetse(nt,j-1)) > 4*max(abs(vnetse(nt,j-1)- ...
                vnetse(nt,j-2)),4)
            vnetse1(1,j)=NaN;
        end    
        if j > 1 && vnetse(nt,j-1)-vnetse(nt,j) > 0.2*vnet(nt,j) ...
            && isnan(vold) == 0 && abs(vold-vnet(nt,j)) < 0.15*vnet(nt,j)
           vnetse1(1,j)=NaN; 
           m=m+1;
           jcon(m)=j;
           vcon(m,1)=vold;
           vcon(m,2)=vnet(nt,j);
        end
        vold=vnetse1(1,j);
    end
else
    vnetse1(vnetse1==0)=NaN;
end  
%
[~,jmin]=min(vnet(nt,:));
jmin=jmin-1;
%
j=plot(time(1:jmin),wunitsfac*vnet(nt,1:jmin),'b',time(1:jmin),wunitsfac*vnetse1(1,1:jmin),'r');
set(j,'linewidth',2)
%
if datenum(TimeStamp) < 7.3685e+05 -1
    hold on
    for i=1:m
        %j1=plot(time(jcon(i)-1:jcon(i)),vcon(i,:),'color',[0 0.5 0],'linewidth',2);
        j1=plot(time(jcon(i)-1:jcon(i)),wunitsfac*vcon(i,:),'r','linewidth',2);
    end    
    hold off
end
%        
u=get(gca,'xticklabel');
v=get(gca,'xtick');
un=str2double(u); 
a=1+mod(un-1,days(monthstore(nt,1)));
as=num2str(a);
set(gca,'xtick',v,'xticklabel',as)
set(gca,'xlim',[v(1),v(end)])
set(gca,'fontname',gfont,'fontsize',gfontsize,'fontweight',gfontweight)
if sum(a) == sum(un)
    xtit=months(monthstore(nt,1));
else
    monp=monthstore(nt,1)+1;
    if monp == 13
        monp=1;
    end    
    xtit=[char(months(monthstore(nt,1))),'-',char(months(monp))];
end   
if exist('yearstore','var') == 1 
    xtit=[char(xtit), ' ', char(num2str(yearstore(nt)))];
end  
xlabel(xtit)
ylabl=char(strcat('Maximum surface wind speed (',wunitslabel,')'));
ylabel(ylabl,'fontweight','bold','fontsize',12)
if sum(vnetse(nt,:)) > 0
    legend('Primary','Secondary','location','northwest')
end    
datea(1)=upper(datea(1));
ntstring=num2str(nt);
titstring2=['track number ' ntstring  ];
tit1=[datea ' ' titstring2];
title(tit1,'fontsize',12,'fontweight','bold','interpreter','none')
%
% Add button to access documentation of script
%
if strcmp(Docdisp,'y')
    tempname=mfilename;
    mname = strcat('graphics_docs/html/',tempname,'_doc.html');
    comstr=['open ',mname,';delete(hButton)'];
    labelStr = '<html><center><a href="">Doc';
    hButton = uicontrol('string',labelStr,'pos',[10,10,50,20],'callback',comstr); 
end
%
% Clean up
%
if strcmp(Clearvlbs,'y')
    clear vnetse1 jcon vcon vold time a as u un v
end    