%
%  This script develops a windfield that includes secondary eyewalls (if
%  present) and plots the evolution of the radii of maximum winds and an 
%  outer radius. 
%
%  Copyright WindRiskTech, LLC, 2014
%  Last modified June, 2017
%-------------------------------------------------------------------------
%
vouterthresh=25;   % Circular wind component (knots) defining outer storm radius
%
%-------------------------------------------------------------------------
%
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 vstore vsestore rmstore rmsestore hourstore daystore monthstore ...
    shape city_radius
%
rmstore1=magfac*rmstore;
rmsestore1=magfac*rmsestore;
clear vse rmse
%  
% Randomize radii of maximum wind
%
if strcmp(randfac,'y')
    [nrm,mrm]=size(rmstore);
    rfac=outer(nrm,mrm);
    rmstore1=rmstore1.*rfac;
    rmsestore1=rmsestore1.*rfac;
end
load sorted
clf
%
%  Prompt, if necessary, for storm rank number
%
if ~exist('nt','var')
   nt=get_storm_num;
end   
%
vse(1,:)=vsestore(nt,:);
rmse(1,:)=rmsestore1(nt,:);
[ni,nj]=size(vstore);
r=0:2:1000;
qr=max(size(r));
V=zeros(qr,nj);
router=zeros(1,nj);
%
for j=1:nj
    %V(:,j)=windprofilem(vstore(nt,j),rmstore1(nt,j),vsestore(nt,j), ...
    %    rmsestore1(nt,j),r,3);
    V(:,j)=windprofilem(vstore(nt,j),rmstore1(nt,j),0,0,r,3);
    [Vmax,kVmax]=max(V(:,j));
    [~,kout]=min(abs(V(kVmax:end,j)-0.95*vouterthresh));
    if Vmax > 1.2*vouterthresh
        router(j)=r(kout);
        router(j)=max(router(j),rmstore(nt,j));
    end    
end
router(router==0)=NaN;
%
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
    rcon=zeros(nj,2);
    vcon=zeros(nj,2);
    jcon=zeros(1,nj);
    rmold=1;
    vold=1;
    m=0;
    for j=1:nj
        if rmsestore1(nt,j) == 0
            rmse(1,j)=NaN;
            vse(1,j)=NaN;
        end
        if j > 2 && vsestore(nt,j-1) > 0 && vsestore(nt,j-2) > 0 && ...
            abs(vsestore(nt,j)-vsestore(nt,j-1)) > 4*max(abs(vsestore(nt,j-1)- ...
            vsestore(nt,j-2)),4)
          vse(1,j)=NaN;
          rmse(1,j)=NaN;
        end    
        if j > 1 && vsestore(nt,j-1)-vsestore(nt,j) > 0.2*vstore(nt,j) ...
            && isnan(vold) == 0 && abs(vold-vstore(nt,j)) < 0.15*vstore(nt,j)
           rmse(1,j)=NaN;
           vse(1,j)=NaN;
           m=m+1;
           jcon(m)=j;
           rcon(m,1)=rmold;
           rcon(m,2)=rmstore1(nt,j);
           vcon(m,1)=vold;
           vcon(m,2)=vstore(nt,j);
        end
        rmold=rmse(1,j);
        vold=vse(1,j);
    end
else
    rmse=rmsestore1(nt,:);
    rmse(rmse==0)=NaN;
end    
%
[~,jmin]=min(vstore(nt,:));
jmin=jmin-1;
%
if datenum(TimeStamp) < 7.3685e+05 -1
    h=plot(time(1:jmin),rmstore1(nt,1:jmin),'b',time(1:jmin), ...
    rmse(1,1:jmin),'r'); 
    set(h,'linewidth',2)
else    
    h=plot(time(1:jmin),rmstore1(nt,1:jmin),'.b',time(1:jmin), ...
       rmse(1,1:jmin),'.r','markersize',9); 
end   
%
hold on
if datenum(TimeStamp) < 7.3685e+05 -1
    for i=1:m
        h1=plot(time(jcon(i)-1:jcon(i)),rcon(i,:),'r','linewidth',2);
    end    
end
jj=plot(time(1:jmin),router(1:jmin),'c','linewidth',1.5);
hold off
%        
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)
ylabel('Radii (km)')
legtit=['Radius of ',num2str(vouterthresh),' kt wind'];
if sum(vsestore(nt,:)) > 0
    legend([h;jj],'Primary','Secondary',legtit,'location','best')
else
    legend([h(1,1);jj],'Primary',legtit,'location','best')
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
%
clear rmsestore1 rmstore1 rmse vmse V router r jcon rcon vcon rmold vold u v