%
%  Make return period plot for lifetime maximum wind speed or maximum wind 
%  speed within city_radius of point  of interest, if shape ='circ', else 
%  for maximum winds crossing line segements, if shape = 'poly'. In the first
%  case, quantity plotted is set by the parameter 'htype' in the params.m
%  file. Intensity limits and resolution can also be adjusted in params. 
%  February, 2009; last modified September, 2021
%  Copyright 2009 WindRiskTech L.L.C.
%-------------------------------------------------------------------------
params    %  Load parameters
%  
clf('reset')
if exist('shearstore','var') == 0
    clear vmax
    shape='circ';
    load temp
end
load('sorted','storm','datea')
[n,m]=size(vstore);
%
nswitch=0;
if strcmp(shape,'circ') && city_radius < 2000
    load('temp.mat','vsmax')
end    
if (strcmp(shape,'poly') || city_radius >= 2000) && htype == 3
    htype=2;
end  
kmax3=kmax2;
if htype == 1
    vtemp=vnetmax;
    [~,kmax3]=max(vnet,[],2);
elseif htype == 2
    vtemp=vmax;
else
    vtemp=vsmax;
end    
%
if strcmp(wreduc,'y') && htype < 3
    logfac=log(wheight/500)/0.35; % Factor to Reduce gradient wind to 10 m wind
    %
    % Load neutral drag coefficients
    %
    load C_Drag500 cd
    load bathymetry bathy
    mincd=min(min(cd));
    cd(bathy<0)=mincd; % This statement added to set ocean point drag coefficients to 
                       % their proper value; this treats a few bogus points
                       % in the drag coefficient data set resulting from
                       % ECMWF's interpolation from a 1 x 1 grid
    clear bathy
    rat=1./(1+sqrt(mincd)*logfac); % Factor to ensure correct wind speeds over water    
    %
    % Interpolate drag coefficient to storm position
    %
    for i=1:n
        ib=1+floor(4.*longstore(i,kmax3(i)));
        ibp=ib+1;
        if ibp > 1440
            ibp=1;
        end    
        jb=1+floor(4.*(latstore(i,kmax3(i))+90));
        b1=cd(ib,jb);
        b2=cd(ib,jb+1);
        b3=cd(ibp,jb);
        b4=cd(ibp,jb+1);
        %
        %  This block interpolates drag coefficient to latitude and longitude of max wind
        %
        dely=1+4.*(latstore(i,kmax3(i))+90)-jb;
        delx=1+4.*longstore(i,kmax3(i))-ib;
        d1=(1.-delx).*(1.-dely);
        d2=dely.*(1.-delx);
        d3=delx.*(1.-dely);
        d4=delx.*dely;
        cdrag=1./(d1./b1+d2./b2+d3./b3+d4./b4); % Conservative estimate; weights small roughness
        cdfac=max(1+sqrt(cdrag).*logfac,0);
        vtemp(i)=vtemp(i).*rat.*cdfac;
        %
    end
end  
%
if strcmp(hbins,'ss')
    x=[39,64,83,96,113,137]; % Changed 40 to 39, July 2020
    xb=x;
    xb(end+1)=max(vend,250); %  Arbitrary upper bound for Cat 5 needed for plotting pdfs
elseif strcmpi(hbins,'auto')
    load temp vcrit
    x=str2double(vcrit)-1:vres:10*ceil(0.1*max(vtemp));
    xb=x;
    xb(end+1)=xb(end)+vres;
else    
    x=vbeg-1:vres:vend;      % Changed vbeg to vbeg-1 July 2020
    xb=x;
    xb(end+1)=vend+vres;
end    
xplot=x;
if strcmp(histf,'pdf')
    for i=1:max(size(x))
        xplot(i)=0.5*(xb(i)+xb(i+1));
    end
end    
if strcmp(hbins,'ss')
    xplot=1:6;
end  
%
histsize=max(size(x));
%
%  Next lines added February, 2014, to account for varying frequency of
%  events
%
ny=1;
ib=1;
if exist('yearstore','var')
    atemp=ismember((min(yearstore):max(yearstore)),yearstore);
    atemp=cast(atemp,'like',yearstore);
    yearset=nonzeros(atemp.*(min(yearstore):max(yearstore)));
    yeartemp=yearstore(1);
    ny=max(size(freqyear));
    vhisttemp=zeros(ny,max(size(x)));
    nss=zeros(1,ny);
    mss=ceil(n/ny);
    vnew=zeros(ny,mss);
    for i=1:n 
        [~,nyy]=ismember(yearstore(i),yearset);
        nss(nyy)=nss(nyy)+1;
        vnew(nyy,nss(nyy))=vtemp(i);
    end  
    sumtotal=0;
    for i=1:ny
        [Nf,~]=histcounts(vnew(i,:),xb);  % Replaced histc with histcounts, May, 2017
        vhisttemp(i,:)=freqyear(i)*Nf;
        sumtotal=sumtotal+sum(Nf);
    end    
    vhistp=sum(vhisttemp,1)/n;
else
    [Nf,~]=histcounts(vtemp,xb);nn% Replaced histc with histcounts, May, 2017
    vhistp=freq*Nf./n;
    sumtotal=sum(Nf);
end
y=zeros(1,histsize);
if strcmp(histf,'pdf')
    y=vhistp;
else    
    for i=1:histsize
        y(i)=sum(vhistp(1,i:end));
    end   
end
if strcmp(hbins,'ss') 
    h=bar(xplot,(1./(1-exp(-y))),0.5);
    set(gca,'yscale','log')
    ylimits=get(gca,'ylim');
    ylimits(1)=1;
    ylimits(2)=10^(0.1*ceil(10.*log10(ylimits(2))));
    set(gca,'ylim',ylimits)
    set(gca,'fontname',gfont,'fontsize',gfontsize,'fontweight',gfontweight)
    set(gca,'xticklabel',{'TS','1' '2' '3' '4' '5'})
    xlabel('Saffir-Simpson Category','fontweight','bold')
else    
    h=semilogy(wunitsfac*xplot,1./(1-exp(-y)),'.b','markersize',24);
    set(gca,'fontname',gfont,'fontsize',gfontsize,'fontweight',gfontweight)
    set(gca,'xlim',wunitsfac*[min(x)-10 max(x)+10])
    xlabl=char(strcat('Maximum Wind Speed (',wunitslabel,')'));
    xlabel(xlabl,'fontweight','bold')
    if strcmp(rgridx,'y')
        set(gca,'xgrid','on')
    end
    if strcmp(rgridy,'y')
         set(gca,'ygrid','on')
    end      
end    
%
ylabel('Return Period (years)','fontweight','bold')
if isempty(storm)
    datea(1)=upper(datea(1)); 
    tit1=datea;
else
    tit1=[storm(1:2) ' ' storm(3:4) ' ' datea(1:4) ' ' datea(5:6) ' ' datea(7:8)...
    ' ' datea(9:10) ' GMT'];
end 
if strcmp(shape,'circ') 
    if htype == 2
        tit2=['All events within ' num2str(city_radius) ' km of POI, '];
        fsize=10;
    elseif htype == 3
        tit2='All events at POI, ';
        fsize=11;
    else
        tit2='Storm Lifetime Maximum Wind, All Events in ';
        fsize=11;
    end    
else
    if strcmp(polytype,'closed')
        tit2='All events within closed polygon, event set ';
        fsize=11;
    else    
        tit2='All events at points of intersection with line segments in poly.in, ';
        fsize=10; 
    end
end
titnet=[tit2 tit1];
title(titnet,'fontsize',fsize,'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
clear cd
xadj=wunitsfac*xplot;
yadj=1./(1-exp(-y));
pfield='wind';
%
% Clean up
%
if strcmp(Clearvlbs,'y')
    clearvars('-except',protected_variables{:})
end 