function [vs, dir, dayk] = pointseries(latstore,longstore,vstore,rmstore, ...
    monthstore,daystore,hourstore,ut,vt,jmax,plat,plong)
%
%  This function take the nn X 400 matrices for the translation speeds 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. The routine utrans.m should be run before this one.
%      K. Emanuel  11/8/2008
%      Copyright 2009 WindRiskTech L.L.C.
%      Last modified January, 2011
%-------------------------------------------------------------------------
% Reduction and translation factors
%
reduc=1;
rawtrans=0;
trans=0.7;
vcrit=0;
%
%  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 2
%
ntime=timelength./timeres+1;
nsteps=2./timeres;
delj=floor(timelength./4);
%
% Outer radius (at which wind speed vanishes)
%
r0=1200;
%
%-------------------------------------------------------------------------
if plong < 0
    plong=plong + 360;
end    
%
[nn,m]=size(ut);
vs=zeros(nn,ntime);
dir=zeros(nn,ntime);
dayk=zeros(nn,ntime);
%
pifac=acos(-1)/180;
dfac=60.*1.852;
dtfac=1./(2.*1.852);
%
dx=dfac.*cos(pifac.*plat).*(plong-longstore);
dy=dfac.*(plat-latstore);
r=sqrt(dx.^2+dy.^2);
for n=1:nn,
    clear dx1 dx2 dy1 dy2 ep x1 y1 rmin rcity jmintot
    dx1=dx(n,1:jmax(n)-1);
    dx2=dx(n,2:jmax(n));
    dy1=dy(n,1:jmax(n)-1);
    dy2=dy(n,2:jmax(n));
    %
    ep(1:jmax(n)-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*(398-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).*vstore(n,j1)+frac.*vstore(n,j2);
       rmfrac=(1-frac).*rmstore(n,j1)+frac.*rmstore(n,j2);
       dxfrac=(1-frac).*dx(n,j1)+frac.*dx(n,j2);
       dyfrac=(1-frac).*dy(n,j1)+frac.*dy(n,j2);   
       rfrac=sqrt(dxfrac.^2+dyfrac.^2);
       ufrac=(1-frac).*ut(n,j1)+frac.*ut(n,j2);
       vfrac=(1-frac).*vt(n,j1)+frac.*vt(n,j2);
       v=windprofile(rfrac,rmfrac,r0,vmfrac);
       v=v.*latstore(n,1)./(abs(latstore(n,1))+1e-8);
       un=-v.*dyfrac./max(rfrac,1e-8)+trans.*ufrac;
       vn=v.*dxfrac./max(rfrac,1e-8)+trans.*vfrac;
   %
       rfac=rfrac/sqrt(rfrac.^2+4.*radcity.^2);
       vs(n,k)=sqrt(un.^2+vn.^2)-trans.*rfac.*sqrt(ufrac.^2+vfrac.^2);
       vs(n,k)=max(vs(n,k),0);
       tempd=360+atan2(-un,-vn)./pifac;
       dir(n,k)=mod(tempd,360);
       dayk(n,k)=datenum(2001,monthstore(n,jmintot),daystore(n,jmintot),hourstore(n,jmintot),0,0)+k.*timeres./24;
    end    
 %   
    end
end    