integer i,j,n_wds,ntotal_ocls,num_wds,num_ocls,n_match_p1,k,l
      integer n_match_p2,n_match_p3,n_match_p4,n_match_p5,n_match_p6
      integer flag
      
      parameter (ntotal_ocls=2000)
      parameter (ntotal_wds=400000)
      
      double precision alpha_ocl(ntotal_ocls),delta_ocl(ntotal_ocls),
     &                 rad_50(ntotal_ocls),pi_ocl(ntotal_ocls),
     &                 s_pi_ocl(ntotal_ocls),pma_ocl(ntotal_ocls),
     &                 s_pma_ocl(ntotal_ocls),pmd_ocl(ntotal_ocls),
     &                 s_pmd_ocl(ntotal_ocls)
      double precision alpha_wd,delta_wd,pi_wd,s_pi_wd,pma_wd,s_pma_wd,
     &                 pmd_wd,s_pmd_wd
      double precision pi,sigma_diam,sigma_pms,sigma_pi
      double precision gamma,diam_distance,limit_diam,limit_pi
      double precision pma_distance,pmd_distance,pi_distance
      double precision limit_pma,limit_pmd
      double precision x1_test,y1_test,test_distance
      double precision epsilon,angle,test_area
      
      parameter (pi=3.14159265359)
      parameter (sigma_diam=4.5,sigma_pms=3.,sigma_pi=3.)
      parameter (epsilon=0.05)
      
      character*23 name_wd,name_ocl(ntotal_ocls)
      
      open(unit=1,file="list_ocls.txt",status="old")
      
3     format(a23,1x,a23,1xi5,1x,i7,1x,8(f7.3,1x),i2)
      
      num_ocls=0
      num_wds=0
      n_match_p1=0
      n_match_p2=0
      n_match_p3=0
      n_match_p4=0
      n_match_p5=0
      n_match_p6=0
      
      do 10 i=1,ntotal_ocls
      
      read(1,*,end=11) name_ocl(i),alpha_ocl(i),delta_ocl(i),rad_50(i),
     &                 pi_ocl(i),s_pi_ocl(i),pma_ocl(i),s_pma_ocl(i),
     &                 pmd_ocl(i),s_pmd_ocl(i)
     
      alpha_ocl(i)=alpha_ocl(i)*(pi/180.)
      delta_ocl(i)=delta_ocl(i)*(pi/180.)

      num_ocls=num_ocls+1
      
10    continue

11    close(1)

      write(*,*) num_ocls," OCLs read"
      
      open(unit=1,file="list_wds.txt",status="old")
      
      open(unit=2,file="list_match.txt",status="replace")
      
      do 100 j=1,ntotal_wds
      
      read(1,*,end=110) name_wd,alpha_wd,delta_wd,pi_wd,s_pi_wd,pma_wd,
     &                  s_pma_wd,pmd_wd,s_pmd_wd
     
      alpha_wd=alpha_wd*(pi/180.)
      delta_wd=delta_wd*(pi/180.)
      
      num_wds=num_wds+1
      
      do 50 i=1,num_ocls
      
      flag=0

      gamma=(dsin(delta_wd)*dsin(delta_ocl(i)))+(dcos(delta_wd)*
     &       dcos(delta_ocl(i))*dcos((alpha_wd-alpha_ocl(i))))
      
      diam_distance=dacos(gamma)*(180/pi)
      limit_diam=sigma_diam*rad_50(i)
      
      if(diam_distance .ge. limit_diam)then
      goto 50
      endif
      
      n_match_p1=n_match_p1+1
      
      pi_distance=abs(pi_wd-pi_ocl(i))
      limit_pi=(sigma_pi*s_i_wd)+(sigma_pi*s_pi_ocl(i))

      if(pi_distance .ge. limit_pi)then
      goto 50
      endif
      
      n_match_p2=n_match_p2+1
      
      pma_distance=abs(pma_wd-pma_ocl(i))
      limit_pma=(sigma_pms*s_pma_wd)+(sigma_pms*s_pma_ocl(i))
      
      if(pma_distance .ge. limit_pma)then
      goto 50
      endif

      n_match_p3=n_match_p3+1
      
      pmd_distance=abs(pmd_wd-pmd_ocl(i))
      limit_pmd=(sigma_pms*s_pmd_wd)+(sigma_pms*s_pmd_ocl(i))

      if(pmd_distance .ge. limit_pmd)then
      goto 50
      endif

      n_match_p4=n_match_p4+1
      
      do 40 k=1,7200
      
      angle=(dble(k)*0.05)*(pi/180.)

      x1_test=pma_ocl(i)+(sigma_pms*s_pma_ocl(i)*dcos(angle))
      y1_test=pmd_ocl(i)+(sigma_pms*s_pmd_ocl(i)*dsin(angle))
      
      test_area=(((x1_test-pma_wd)/(sigma_pms*s_pma_wd))**2)+
     &          (((y1_test-pmd_wd)/(sigma_pms*s_pmd_wd))**2)-1

      if(test_area .le. 0)then

      n_match_p5=n_match_p5+1

      test_area=(((pma_wd-pma_ocl(i))/(sigma_pms*s_pma_ocl(i)))**2)+
     &          (((pmd_wd-pmd_ocl(i))/(sigma_pms*s_pmd_ocl(i)))**2)-1

      if(test_area .le. 0)then
      
      flag=1
      
      n_match_p6=n_match_p6+1
      
      endif

      write(2,3) name_wd,name_ocl(i),i,j,diam_distance,limit_diam,
     &           pma_distance,limit_pma,pmd_distance,limit_pmd,
     &           pi_distance,limit_pi,flag
     
      goto 50
      
      endif
      
40    continue

50    continue
      
100   write(*,*) j," done ",n_match_p1,n_match_p2,n_match_p3,n_match_p4,
     &           n_match_p5,n_match_p6," matches"

110   close(1)

      close(2)

      write(*,*) num_wds," WDs read"

      stop
      end