function [vs, dir] = pointshort(latstore,longstore,vstore,rmstore, ...
   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.
%  Choice of three wind profiles:  Be sure to comment/uncomment appropriate
%  code in both parameter list and body of code. 
%      K. Emanuel  11/8/2008
%      Copyright 2009 WindRiskTech L.L.C.
%      Last modified February, 2011
%-------------------------------------------------------------------------
% Reduction and translation factors
%
reduc=1;
rawtrans=0;
trans=0.7;
vcrit=0;
%
%  Parameters for making detailed wind records at fixed points
%
radcity=400;    % Distance (km) of storm from point of interest outside of which storm ignored
timelength=60;  % Total length of time series (hours)
timeres=0.5;    % 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;
%
%  Coriolis
%
f=5.0e-5;
%
%bs=1.8;
%rn=300./20;     % Holland parameters
%xn=0.8;
%
%  Emanuel (2006) parameters
%b=0.25;
%nb=0.9;
%mb=1.6;  
%
%mb2=2.*mb;
%fac1=(1-b).*(mb+nb);
%fac2=b.*(1+2.*mb);
%fac3=2.*(mb+nb);
%fac4=2.*mb+1; 
%
%-------------------------------------------------------------------------
if plong < 0
    plong=plong + 360;
end  
%
[nn,m]=size(ut);
vs=zeros(nn,ntime);
dir=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)+1e-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);
       %
       % Holland (2010) wind profile
       %
       %q=max(size(r));
       %
       %x=zeros(1,q);
       %rh=rfrac./rmfrac;
       %for i=1:q,
       %     if rh(i) <= 1
       %         x(i)=0.5;
       %     else
       %         x(i)=0.5+(rh(i)-1).*(xn-0.5)./(rn-1);
       %     end
       %end
       %
       %v=vmfrac.*(rh.^-bs.*exp(1-rh.^-bs)).^x;
       %
       % Emanuel (2006) wind profile
       %
       %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)));
       %
       % ER 2011 wind profile
       %
       Mm=rmfrac.*vmfrac+0.5.*f.*rmfrac.^2;
       rn=rfrac./rmfrac;
       M=Mm.*(2.*rn.^2./(1+rn.^2));
       v=(M-0.5.*f.*rfrac.^2)./(rfrac+1e-8);
       %
       % End of wind profile options
       %
       v=max(v,0);
       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+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);
    end    
 %   
    end
end    