%
%  This script calculates and graphs the distribution of 
%  wind speed, in knots, for a given storm at a specified time.
%
%       K. Emanuel  1/20/2011
%       Copyright 2011 WindRiskTech L.L.C.
%
%------------------------------------------------------------------------
%  Options
%  Plot storm track?
    plottrack='y';
% Coastline resolution(0 for none, 1 for very low, 2 for low, 3 for medium,  
% 4 for high. Note that higher resolution takes more memory and CPU and 
% slows down the plot. You must separately download high resolution data bases
% for choices 3 and 4.)
	cres=3;
% Color of shoreline ('none' or MATLAB color options)
	shorecolor='none';
% Ocean color 
    oceancolor=[0.95 0.95 1];
% Color of land (only effective when no topography contoured)
    landcolor=[0.45 0.6 0.45];
% Axis label font
	axisfont='arial';
% Axis label font size
	axisfontsize=12;
% Axis label font weight
    axisfontweight='bold';
% Grid line style ('-',':','-.','--', or 'none')
    gridline='none';
% Map projection (Options include Mercator, Miller cylindrical, Equidistant
%  cylindrical, Transverse Mercator, Sinusoidal, Gall-Peters, 
%  Lambert conformal conic, Albers equal-area conic, Hammer-Aitoff,
%  and Mollweide. Other options are possible, but the call statement then
%  has to be changed; see M_Map documentation.)
    projection='gall-peters';
    % Plot ../chip/data/states (y or n). Need to download political boundaries from 
%  http://www.maproom.psu.edu/cgi-bin/ian/points/index.cgi
%  to use this
    pstates='y';
% Color of state boundaries (MATLAB color options)
    stcolor=[.4 .4 .4];
% Automatic (auto) or manual (manu) specification of map bounds
    bound='auto';
% Latitude and longitude distances of map boundaries from storm 
%   center (degrees) if bound = 'auto'
    deltax=4;
    deltay=3;
% Specify map bounds if bounds='manu'
    xmin=20;
    xmax=380;
    ymin=-60;
    ymax=60;
% Spatial resolution (degrees)
    delx=0.05;
%
%  Vortex radial profile
%
r0=1200;    % Outer radius (at which wind speed vanishes)
%
trans=0.7;  % Fraction of translation speed to add to storm winds
radcity=300;    % Distance (km) of storm from point of interest outside of
                % which storm ignored
magfac=1;   % Scale factor for radius of maximum winds
%
%--------------------------------------------------------------------------
%
r0=r0.*sqrt(magfac);
clf('reset')
shape='circ';
load('temp.mat')
fs3=['sorted'];
load(fs3)
%
if isempty(storm)
   nswitch=input('Order by maximum at POI (1), lifetime maximum (2), or random (3)?');
else    
   nswitch=2;
end  
nttemp=input('Input track number');
%
monthplot=input('Enter month');
dayplot=input('Enter day of the month');
hourplot=input('Enter GMT hour (may include decimal digits)');
%
if nswitch == 1
   nt=nsortp(nttemp);
elseif nswitch == 2   
   nt=nsort(nttemp);
else
   nt=nttemp;
end   
ntstring=num2str(nt);
%
[utd,vtd,jmaxd] = utrans(latstore,longstore);
%
[n,m]=size(utd);
%
clear timev V
dum=zeros(1,jmaxd(nt));
V=[dum+2000; monthstore(nt,1:jmaxd(nt)); daystore(nt,1:jmaxd(nt));...
    hourstore(nt,1:jmaxd(nt)); dum; dum];
V=transpose(V);
timev=datenum(V);
timeplot=datenum([2000 monthplot dayplot hourplot 0 0]);
[dum2,jplot]=min(max(timeplot-timev,0));
jplot=max(jplot,2);
weight=(timeplot-timev(jplot-1))./(timev(jplot)-timev(jplot-1));
%
latstorm=weight.*latstore(nt,jplot)+(1-weight).*latstore(nt,jplot-1);
longstorm=weight.*longstore(nt,jplot)+(1-weight).*longstore(nt,jplot-1);
vstorm=weight.*vstore(nt,jplot)+(1-weight).*vstore(nt,jplot-1);
rmstorm=magfac.*(weight.*rmstore(nt,jplot)+(1-weight).*rmstore(nt,jplot-1));
utstorm=weight.*utd(nt,jplot)+(1-weight).*utd(nt,jplot-1);
vtstorm=weight.*vtd(nt,jplot)+(1-weight).*vtd(nt,jplot-1);
%
pifac=acos(-1)/180;
dfac=60.*1.852;
%
if char(bound) == char('auto')
    xmin=longstorm-deltax;
    xmax=longstorm+deltax;
    ymin=latstorm-deltay;
    ymax=latstorm+deltay;
end
%
x=[xmin:delx:xmax];
y=[ymin:delx:ymax];
sx=max(size(x));
sy=max(size(y));
vspeed=zeros(sx,sy);
for i=1:sx,
    dx=dfac.*cos(pifac.*latstorm).*(x(i)-longstorm);
    for j=1:sy,
        dy=dfac.*(y(j)-latstorm);
        r=sqrt(dx.^2+dy.^2);
        vcirc=windprofile(r,rmstorm,r0,vstorm);
        vcirc=vcirc.*latstorm./(abs(latstorm)+1e-8);
        un=-vcirc.*dy./max(r,1e-8)+trans.*utstorm;
        vn=vcirc.*dx./max(r,1e-8)+trans.*vtstorm;
   %
        rfac=r/sqrt(r.^2+4.*radcity.^2);
        vspeed(i,j)=sqrt(un.^2+vn.^2)-trans.*rfac.*sqrt(utstorm.^2+vtstorm.^2);
        vspeed(i,j)=max(vspeed(i,j),0);
    end
end    
vspeed=transpose(vspeed);
%
clf
m_proj(projection,'long',[xmin xmax], ...
   'lat',[ymin ymax]);
set(gca,'color',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 pstates == 'y'
   m_states(xmin,xmax,ymin,ymax,'color',stcolor)
end
%
vvm=max(max(vspeed));
vvmm=10.*ceil(0.1.*vvm);
vinc=10;
if vvmm <= 100
    vinc=5;
end    
m_contour(x,y,vspeed,[20:vinc:vvmm],'linewidth',2)
m_gridm('fontname',axisfont,'fontsize',axisfontsize,'fontweight',axisfontweight,'linestyle',gridline)
if plottrack == 'y'
    longplot=zeros(1,jplot);
    longplot(1,1:jplot-1)=longstore(nt,1:jplot-1);
    longplot(1,jplot)=longstorm;
    latplot=zeros(1,jplot);
    latplot(1,1:jplot-1)=latstore(nt,1:jplot-1);
    latplot(1,jplot)=latstorm;
    m_plot(longplot,latplot,'k','linewidth',3)
end    
hold off
cb=colorbar;
set(cb,'fontweight','bold')
set(get(cb,'Ylabel'),'String','Knots','fontweight','bold')
if isempty(storm)
    tit1=datea;
    adate=datestr(timeplot,'mmmm dd, HH:MM');
    tit1={tit1;['Track number ' ntstring ', ' adate ' GMT']};
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')
