%
%  Make exceedence 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 December, 2020
%  Copyright 2009 WindRiskTech L.L.C.
%-------------------------------------------------------------------------
params    %  Load parameters
load temp kmax2 vstore
%  
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;
kmax3=max(kmax3,1);
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: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); % Replaced histc with histcounts, May, 2017
    vhistp=freq*Nf./n;
    sumtotal=sum(Nf);
end
[Nfsyn,~]=histcounts(vtemp,xb);
y=zeros(1,histsize);
if strcmp(histf,'pdf')
    y=vhistp;
else    
    for i=1:histsize
        y(i)=sum(vhistp(1,i:end));
    end   
end
h=bar(wunitsfac*xplot,y);
set(gca,'fontname',gfont,'fontsize',gfontsize,'fontweight',gfontweight)
if strcmp(hbins,'ss') 
    set(gca,'xticklabel',{'TS' '1' '2' '3' '4' '5'})
    xlabel('Saffir-Simpson Category','fontweight','bold')
else    
    clear u us us2
    set(gca,'xlim',[min(wunitsfac*x)-10 max(wunitsfac*x)+10])
    xlabl=char(strcat('Maximum Wind Speed (',wunitslabel,')'));
    xlabel(xlabl,'fontweight','bold')
end    
%
if strcmp(histf,'cum')
    ylabel('Annual Exceedence Frequency','fontweight','bold')
else
    ylabel('Annual Frequency Density','fontweight','bold')
end   
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
%
% Clean up
%
if strcmp(Clearvlbs,'y')
    clearvars('-except',protected_variables{:})
end 