function [pass, vmax, jmax] = best_filter(vsm,latm,longm,city_radius,clat,clong)
%
%  Thus function finds those tracks that pass within city_radius of
%  [clat, clong] and also finds the maximum wind speed (vmax) within that
%  circle, occuring at jmax(i). Pass= 1 if track meets the criterion, 0 if  
%  it does not. Use best_filter_poly if line segment filtering is desired.
%  Copyright 2009 WindRiskTech L.L.C.
%  Modified August 2009 to include jmax.
%  Last modified March, 2012
%-------------------------------------------------------------------------
fac=acos(-1)/180;
fac2=1.852*60;
cosfac=cos(fac*clat);
if clong < 0
    clong=clong+360;
end    
[iq,~]=size(vsm);
vmax=zeros(1,iq);
pass=zeros(1,iq);
jmax=zeros(1,iq)+1;
[~,kmax]=min(abs(latm),[],2);
kmax=kmax-1;
kmax=max(kmax,1);
%    
for i=1:iq
      pass2=0;
      vmax2=0;
      jmax2=1;
      for k=1:kmax(i)
            dx1=cosfac*(clong-longm(i,k));
            dy1=clat-latm(i,k);
            kp=k+1;
            kp=min(kp,kmax(i));
            dx2=cosfac*(clong-longm(i,kp));
            dy2=clat-latm(i,kp);
            ep=((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);
            v1=vsm(i,k)+ep*(vsm(i,kp)-vsm(i,k));
            %v1=max(vsm(i,k),vsm(i,kp));
            dis=fac2*sqrt(x1^2+y1^2);
            dis1=fac2*sqrt(dx1^2+dy1^2);
            if dis <= city_radius && v1 > vmax2
                pass2=1;
                vmax2=v1;
                jmax2=k;
            end
            if dis1 <= city_radius && vsm(i,k) > vmax(i)
                pass(i)=1;
                vmax(i)=vsm(i,k);
                jmax(i)=k;
            end        
      end 
      if vmax2 > vmax(i)
          vmax(i)=vmax2;
          pass(i)=pass2;
          jmax(i)=jmax2;
      end    
end
% 