function [vs, dir, dayk] = bestseries(lat,long,vmax,rm, ...
    month,day,hour,plat,plong)
%
%  This function takes a particular best track storm and,
%  using the latitude (plat) and longitude (plong) of the point of interest,
%  calculates the time series of wind speed and direction at the point of
%  interest. 
%      K. Emanuel  7/5/2010
%      Copyright 2010 WindRiskTech L.L.C.
%
%-------------------------------------------------------------------------
transfac=0.7;   % Fraction of translation speed to add to circular winds
%
%  Parameters for making detailed wind records at fixed points
%
radcity=300;    % Distance (km) of storm from point of interest outside of which storm ignored
timelength=24;  % Total length of time series (hours)
timeres=0.25;   % Time resolution of time series (hours) should divide evenly into 6
%
ntime=timelength./timeres+1;
nsteps=6./timeres;
delj=floor(timelength./12);
nj=max(size(lat));
%
%  Vortex radial profile
%
% Outer radius (at which wind speed vanishes)
%
r0=1200;
%
% Shape parameters
%
b=0.25;
nb=0.9;
mb=1.6;
%-------------------------------------------------------------------------
if plong < 0
    plong=plong + 360;
end    
%
mb2=2.*mb;
fac1=(1-b).*(mb+nb);
fac2=b.*(1+2.*mb);
fac3=2.*(mb+nb);
fac4=2.*mb+1;
%
[dum,jmax]=min(vmax);
jmax=jmax-1;
vs=zeros(1,ntime);
dir=zeros(1,ntime);
dayk=zeros(1,ntime);
%
vmax=vmax.*3600./1852;
[ut,vt,jmax]=utransbest(lat,long,month,day,hour);
vmax=vmax-transfac.*sqrt(ut.^2+vt.^2);
%
pifac=acos(-1)/180;
dfac=60.*1.852;
%
dx=dfac.*cos(pifac.*plat).*(plong-long);
dy=dfac.*(plat-lat);
r=sqrt(dx.^2+dy.^2);
for n=1:1,
    clear dx1 dx2 dy1 dy2 ep x1 y1 rmin rcity jmintot
    dx1=dx(1,1:jmax-1);
    dx2=dx(1,2:jmax);
    dy1=dy(1,1:jmax-1);
    dy2=dy(1,2:jmax);
    %
    ep(1:jmax-1)=((dx1.*(dx1-dx2)+dy1.*(dy1-dy2)))./(((dx1-dx2).^2+(dy1-dy2).^2)+1.0e-8);
    ep=max(ep,0);
    ep=min(ep,1);
    x1=dx1+ep.*(dx2-dx1);
    y1=dy1+ep.*(dy2-dy1);
    rcity=sqrt(x1.^2+y1.^2);
    [rmin,jmin]=min(rcity);
    jmintot=jmin-delj;
    jmintot=max(jmintot,1);
    %
    ntime=min(ntime,(1+nsteps*(nj-2-jmintot)));
%
    if rmin <= radcity
%        
    for k=1:ntime,
       j1=jmintot+floor((k-1)./nsteps);
       j2=j1+1;
       frac=(k-1)./nsteps-j1+jmintot;
       vmfrac=(1-frac).*vmax(j1)+frac.*vmax(j2);
       rmfrac=(1-frac).*rm(j1)+frac.*rm(j2);
       dxfrac=(1-frac).*dx(j1)+frac.*dx(j2);
       dyfrac=(1-frac).*dy(j1)+frac.*dy(j2);
       ufrac=(1-frac).*ut(j1)+frac.*ut(j2);
       vfrac=(1-frac).*vt(j1)+frac.*vt(j2);       
       rfrac=sqrt(dxfrac.^2+dyfrac.^2);
       r1=min(rfrac,r0);
       rat=r1./max(rmfrac,1);
       v=vmfrac.*((r0-r1)./(r0-rmfrac)).*sqrt(rat.^mb2.*(fac1./(nb+mb.*rat.^fac3)+fac2./(1+mb2.*rat.^fac4)));
       v=v.*lat(1,1)./(abs(lat(1,1))+1e-8);
       un=-v.*dyfrac./max(rfrac,1e-8)+transfac.*ufrac;
       vn=v.*dxfrac./max(rfrac,1e-8)+transfac.*vfrac;
   %
       rfac=rfrac/sqrt(rfrac.^2+4.*radcity.^2);
       vs(k)=sqrt(un.^2+vn.^2)-transfac.*rfac.*sqrt(ufrac.^2+vfrac.^2);
       vs(k)=max(vs(k),0);
       tempd=360+atan2(-un,-vn)./pifac;
       dir(k)=mod(tempd,360);
       dayk(k)=datenum(double(2001),double(month(jmintot)),double(day(jmintot)),...
           double(hour(jmintot)),double(0),double(0))+k.*timeres./24;
    end    
 %   
    end
end    
end
