%
%  This script calculates and plots a histogram of peak wind speeds
%  experienced a a point given by cooridnates plat and plong owing to the
%  input ensemble of tracks.
%  K. Emanuel, June 2010
%  Copyright WindRiskTech, LLC, 2010
%  Modified for use by AER, October, 2010
%
%-------------------------------------------------------------------------
xhist=[20:10:200];
%-------------------------------------------------------------------------
meth=input('Discrete histogram (h) or Kernal-Smoothed (k)?','s');
%
shape='circ';
load('temp.mat');
load('sorted')
[n,m]=size(vstore);
[ut vt jmax]=utrans(latstore,longstore);
%
plat=1;
while plat ~= 0
posit=inputdlg({'Enter latitude','Enter longitude'},'Position Entry',2);
if isempty(posit)
    break
end    
plat=str2num(char(posit(1)));
plong=str2num(char(posit(2)));
if plat == 0
    break
end    
if plong < 0
    plong=plong+360;
end    
%
[vs, dir, dayk] = pointseries(latstore,longstore,vstore,rmstore, ...
    monthstore,daystore,hourstore,ut,vt,jmax,plat,plong);
vsmax=max(vs,[],2);
if meth == 'h'
    vhist=histc(vsmax,xhist);
    q=max(size(xhist));
    vhist=vhist./n;
    sumh=zeros(1,q);
    for i=1:q,
        sumh(i)=sumh(i)+sum(vhist(i:end));
    end    
    h=bar(xhist,sumh,1);
else    
    f = ksdensity(vsmax,xhist,'function','cdf');
    h=plot(xhist,1-f);
    set(h,'linewidth',3)
end    
set(gca,'fontweight','bold')
set(gca,'xlim',[xhist(1)-10 xhist(end)+10])
xlabel('Wind Speed (knots)','fontweight','bold')
ylabel('Exceedence Probability','fontweight','bold')
tit1=[storm(1:2) ' ' storm(3:4) ' ' datea(1:4) ' ' datea(5:6) ' ' datea(7:8)...
    ' ' datea(9:10) ' GMT'];
title(tit1,'fontsize',12,'fontweight','bold')
end

