%
%  This script calculates and graphs the distribution of maximum point
%  wind speed for a storm.
%  Note:  If discontinuities appear in the result, increase the parameter
%  "timelength" in params.m.
%       January, 2009
%       Copyright 2009 WindRiskTech L.L.C.
%  Last modified, March, 2020
%------------------------------------------------------------------------
%
params    %  Load parameters
clear yq xq va ua g1 lat long v scale
%
load temp shape
if strcmp(randfac,'y')
    if exist('tmr4','var')
        stop(tmr4)
        delete(tmr4)
        clear tmr4
    end   
    tmr4=timer;
    set(tmr4,'startdelay',2);
    set(tmr4,'timerfcn','try; delete(hms); catch err; end')
    hms=msgbox('Caution: Random draw for wind radii turned on; results will differ with each call');
    start(tmr4)
end    
clf('reset')
if exist('shearstore','var') == 0
    shape='circ';
    load('temp.mat')
elseif strcmp(shape,'poly')
    load temp polytype passe xint yint jint kint kfrac    
end
rmstore1=rmstore.*magfac;
rmsestore1=rmsestore.*magfac;
%  
% Randomize radii of maximum wind
%
if strcmp(randfac,'y')
    [nrm,mrm]=size(rmstore);
    rfac=outer(nrm,mrm);
    rmstore1=rmstore1.*rfac;
    rmsestore1=rmsestore1.*rfac;
end
clear utd vtd v rm vse rmse lat long
fs3='sorted';
load(fs3)
%load landtemp
pifac=acos(-1)./180;
%
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
ntstring=num2str(nt);
%
[n,m]=size(vstore);
%
[~,jmaxd]=min(vstore,[],2);
jmaxd(vstore(:,end)~=0)=m;
jmaxd=jmaxd-1;
jmaxd=transpose(jmaxd);
%
clear v lat long rm jmax utd vtd
[~,q]=min(abs(latstore(nt,:)));
q=q-1;
utd(1,:)=transpose(nonzeros(uinc(nt,:)));
vtd(1,:)=transpose(nonzeros(vinc(nt,:)));
lat(1,:)=transpose(nonzeros(latstore(nt,:)));
long(1,:)=transpose(nonzeros(longstore(nt,:)));  
v(1,:)=transpose(nonzeros(vstore(nt,:)));
qv=max(size(v));
rm(1,:)=transpose(nonzeros(rmstore1(nt,:)));
vse(1,:)=vsestore(nt,1:qv);
rmse(1,:)=rmsestore1(nt,1:qv);
if max(long) - min(long) > 330  && mean(nonzeros(longstore)) < 50 % These lines pertain to Mediterranean storms
    long=wrapTo180(long);
    utd(long<0.5)=[];
    vtd(long<0.5)=[];
    lat(long<0.5)=[];
    v(long<0.5)=[];
    rm(long<0.5)=[];
    vse(long<0.5)=[];
    rmse(long<0.5)=[];
    long(long<0.5)=[];
end  
jmax=max(size(long));
for i=1:q
    if long(1) > 200 && long (i) < 50
        long(i)=long(i)+360;
    end
end  
%
if strcmp(bound,'auto')
    if (strcmp(shape,'circ') && city_radius > 2000) || (strcmp(shape,'poly') && strcmp(polytype,'closed'))
        bxmin=min(nonzeros(long))-deltax;
        bxmax=max(nonzeros(long))+deltax;
        bymin=min(nonzeros(lat))-deltay;
        bymax=max(nonzeros(lat))+deltay;
    else    
        if strcmp(bound,'auto')
            if strcmp(shape,'circ')
                bxmin=clong-deltax;                
                bxmax=clong+deltax;
                bymin=clat-deltay;
                bymax=clat+deltay;
            else
                [~,xsize]=size(xint);
                vtemp=0;
                ixmax=1; 
                for ix=1:xsize
                    if kint(nt,ix) > 0 && vnet(nt,kint(nt,ix)) > vtemp
                        vtemp=vnet(nt,kint(nt,ix));
                        ixmax=ix;
                    end
                end    
                xm=xint(nt,ixmax);
                ym=yint(nt,ixmax);
                bxmin=xm-deltax;
                bxmax=xm+deltax;
                bymin=ym-deltay;
                bymax=ym+deltay;
            end
        end  
    end
end
%
bxmin=max(bxmin,0.5);
x=bxmin:dellatlongs:bxmax;
y=bymin:dellatlongs:bymax;
sx=max(size(x));
sy=max(size(y));
if sx == sy
    x(sx+1)=bxmax+dellatlongs;
    sx=sx+1;
end    
vmaxa=zeros(sx,sy);
[vs1, ~] = pointshortn(lat,long,v,rm,vse,rmse,utd,vtd,y,x,timeres);
vmaxa(:,:)=max(vs1(1,:,:,:),[],2);   
vmaxa=smoothb(vmaxa,sx-1,2,sy-1);
vmaxa=transpose(vmaxa);
%vmaxa=vmaxa.*1.15;  % Conversion to MPH
%
clf('reset')
m_proj(projection,'long',[bxmin bxmax], ...
   'lat',[bymin bymax]);
m_gridm('fontname',axisfont,'fontsize',axisfontsize,'fontweight',axisfontweight, ...
    'linestyle',gridline,'backgroundcolor',oceancolor)
ld=landcolor;
hold on
if cres == 1
      m_coast('patch',ld,'edgecolor',shorecolor);
elseif cres == 2  
      m_gshhs_l('patch',ld,'edgecolor',shorecolor);
elseif cres == 3  
      m_gshhs_i('patch',ld,'edgecolor',shorecolor);
elseif cres == 4   
      m_gshhs_h('patch',ld,'edgecolor',shorecolor);
end   
%
if strcmp(pstates,'y')
  m_states(bxmin,bxmax,bymin,bymax,'color',stcolor)
end
%
vvm=max(max(vmaxa));
vvmm=10*ceil(0.1*vvm);
cint=floor((vvm-20)./10);
if verLessThan('matlab', '8.5.0') && strcmp(wfill,'y') %  Apply this to pre 2015 (R2014b) matlab versions
    [C,h]=m_contourf(x,y,wunitsfac*vmaxa,wunitsfac*(20:cint:vvmm),'linewidth',2);
    alpha(0.99)
    hPatch = findobj(h,'Type','patch');
    set(hPatch,'facealpha',wtrans); 
    colormap(wmap)
else
    [C,h1]=m_contour(x,y,wunitsfac*vmaxa,wunitsfac*(20:cint:vvmm),'linewidth',1);
    colormap(wmap)
    if strcmp(wfill,'y')
        h=m_pcolor(x,y,wunitsfac*vmaxa);
        %
        a=colormap(wmap);
        for i=1:6
            a(i,:)=oceancolor;
        end
        colormap(a)
        %
        set(h,'EdgeAlpha',0,'FaceAlpha',wtrans);
        set(h,'FaceColor','interp');
    end
end    
if plottrack == 'y'
    m_plot(long(1:jmax),lat(1:jmax),'k','linewidth',3)
    ua=long(jmax)-long(jmax-1);
    va=lat(jmax)-lat(jmax-1);
    xq=long(jmax-1);
    yq=lat(jmax-1);
    scale=afactor/sqrt(ua^2+va^2);
    ua=ua*scale;
    va=va*scale;
    g1=m_quiver(xq,yq,ua,va,0);
    set(g1,'linewidth',3)
    set(g1,'color','k')     
end    
hold off
%
cb=colorbar;
if verLessThan('matlab', '8.5.0')  %  Apply this to pre 2015 (R2014b) matlab versions
    cb = findobj(gcf,'Type','axes','Tag','Colorbar');
    cbIm = findobj(cb,'Type','image');
    alpha(cbIm,wtrans)
end
set(cb,'fontname',axisfont,'fontweight','bold')
set(get(cb,'ylabel'),'String',wunitslabel,'fontweight','bold','rotation',0,'horizontalalignment','left');
%
if isempty(storm)
    datea(1)=upper(datea(1));
    tit1=datea;
    if exist('yearstore','var') == 1 
        %if exist('tit3','var')
        %    tit1=tit3;
        %else    
            tit1={tit1; ['Track number ' ntstring ', year ' char(num2str(yearstore(nt)))]};
        %end
    else
        tit1={tit1;['Track number ' ntstring]};
    end     
else
    tit1=[storm(1:2) ' ' storm(3:4) ' ' datea(1:4) ' ' datea(5:6) ' ' datea(7:8)...
    ' ' datea(9:10) ' GMT'];
    titstring2=['Track number ' ntstring  ];
    tit2=titstring2;
    tit1=[tit1 '  ' tit2];
end    
title(tit1,'fontsize',12,'fontweight','bold','interpreter','none')
z=vmaxa;
%
% 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')
    clearvars('-except',protected_variables{:})
end 