%
%  This script calculates and graphs the distribution of maximum point
%  wind speed, in knots, for a storm.
%  Note:  If discontinuities appear in the result, increase the parameter
%  "timelength" in the routine "pointshort.m".
%       K. Emanuel  1/1/2009
%       Copyright 2009 WindRiskTech L.L.C.
%  Last modified, October, 2010
%
%------------------------------------------------------------------------
%  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 buffers (degrees) if bound = 'auto'
    deltax=5;
    deltay=2.5;
% Specify map bounds if bounds='manu'
    xmin=20;
    xmax=380;
    ymin=-60;
    ymax=60;
% Spatial resolution
    delx=0.15;
%--------------------------------------------------------------------------
shape='circ';
load('temp.mat')
fs3=['sorted'];
load(fs3)
%
nttemp=input('Input track number');
if isempty(storm)
   nswitch=input('Order by maximum at POI (1), lifetime maximum (2), or random (3)?');
else    
   nswitch=2;
end   
%
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 v ut vt lat long rm jmax
[dum,q]=min(abs(latstore(nt,:)));
q=q-1;
ut(1,:)=utd(nt,:);
vt(1,:)=vtd(nt,:);
lat(1,:)=latstore(nt,:);
long(1,:)=longstore(nt,:);
v(1,:)=vstore(nt,:);
rm(1,:)=rmstore(nt,:);
jmax(1,1)=jmaxd(1,nt);
for i=1:q,
    if long(1) > 200 & long (i) < 50
        long(i)=long(i)+360;
    end
end  
%
if ~isempty(storm) | city_radius > 1000
    xmin=min(nonzeros(long))-deltax;
    xmax=max(nonzeros(long))+deltax;
    ymin=min(nonzeros(lat))-deltay;
    ymax=max(nonzeros(lat))+deltay;
else    
if char(bound) == char('auto')
    if char(shape) == char('circ')
        xmin=clong-deltax;
        xmax=clong+deltax;
        ymin=clat-deltay;
        ymax=clat+deltay;
    else
        file3='poly.in';
        a=load(file3);
        xa=a(:,1);
        ya=a(:,2);
        xb=a(:,3);
        yb=a(:,4);  
        [nseg,dum1]=size(a);
        for i=1:nseg,
            if xa(i) < 0
            xa(i)=xa(i)+360;
            end
            if xb(i) < 0
                xb(i)=xb(i)+360;
            end
        end    
        [pass,xint,yint,jint,kint] = boxm(long,lat,xa,ya,xb,yb);
        xmin=xint-deltax;
        xmax=xint+deltax;
        ymin=yint-deltay;
        ymax=yint+deltay;
    end
end  
end
%
x=[xmin:delx:xmax];
y=[ymin:delx:ymax];
sx=max(size(x));
sy=max(size(y));
vmax=zeros(sx,sy);
for i=1:sx,
    for j=1:sy,
        clear vs dir
        [vs, dir] = pointshort(lat,long,v,rm,ut,vt,jmax,y(j),x(i));
        vmax(i,j)=max(vs(1,:));
    end
end    
vmax=transpose(vmax);
%
clf('reset')
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(vmax));
vvmm=10.*ceil(0.1.*vvm);
m_contour(x,y,vmax,[40:10:vvmm],'linewidth',2)
m_gridm('fontname',axisfont,'fontsize',axisfontsize,'fontweight',axisfontweight,'linestyle',gridline)
if plottrack == 'y'
    m_plot(long(1:jmax),lat(1:jmax),'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;
    tit1={tit1;['Track number ' ntstring]};
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')
