%
%  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    
if exist('vsestore','var') == 0
    load('temp','vstore','vsestore','rmstore','rmsestore','hourstore','daystore', ...
        'monthstore','shape','city_radius')
end    
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') == 0
    if isempty(storm)
        prompt = {'Order by maximum at POI (1), maximum near POI (2), lifetime maximum (3), or by track number (4)?', ...
        'Track number'};
        dlg_title = 'Track Parameters';
        num_lines = 1;
        answer = inputdlg(prompt,dlg_title,num_lines);
        if isempty(answer) || isempty(answer{1})
            return
        end    
        nswitch=str2double(answer{1});
        if (strcmp(shape,'poly') || city_radius >= 2000) && nswitch == 1
            warndlg('Warning: Winds at POI not meaningful for this set; ordering by maximum near POI')
            nswitch=2;
        end                                
        nttemp=str2double(answer{2});
    else    
       nswitch=3;
       prompt={'Track Number'};
       dlg_title='Track Specification';
       num_lines = 1;
       answer = inputdlg(prompt,dlg_title,num_lines);
       if isempty(answer) || isempty(answer{1})
            return
       end
       nttemp=str2double(answer{1});
    end  
    if nswitch == 1
        nt=nsortc(nttemp);
    elseif nswitch == 2
       nt=nsortp(nttemp);
    elseif nswitch == 3   
       nt=nsort(nttemp);
    else
       nt=nttemp;
    end   
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