subroutine which_fiber(harpNum, orb, fiber, thickness)

use bmad
implicit none
type (coord_struct) orb
 
 integer harpNum, fiber
 integer i

 real(rp) xfiber(-3:3)/-0.039,-0.026,-0.013,0.,0.013,0.026,0.039/
 real(rp) yfiber(-3:3)/-0.039,-0.026,-0.013,0.,0.013,0.026,0.039/
 real(rp)  radius/0.00025/
 real(rp) delta
 real(rp) thickness

 ! initialize harp positions

    fiber = 0
    do i=-3,3,1
     if(harpNum == 1 .or. harpNum == 3)delta = orb%vec(1) - xfiber(i)
     if(harpNum == 2 .or. harpNum == 4)delta = orb%vec(3) - yfiber(i)
       if(abs(delta) < radius)then
         fiber = i + 4  ! 1-7
         thickness = 2*sqrt(radius**2 - delta**2)
         exit
       endif
    end do
  return
 end subroutine