program source_point

   use bmad

      type (lat_struct) lat
      type (ele_struct) ele 

      real(rp) x,y,x0,y0,dx,dy,deltas,s, theta_image
      real(rp) theta_0,theta
      real(rp) y_image_0/2.219067/,x_image_0/-9.013192/
      real(rp) rho, r11, r12, r21, r22

      real(rp) tmin, x_tan, y_tan
      real(rp) s_tan
      real(rp) xoff/0./, yoff/0./
      real(rp) x_image, y_image
      real(rp) s_off
      real(rp) orbit_x, orbit_xp, dx0, dy0 
      real(rp) length(2)
      integer ix
      integer nargs,iargc
      integer i
      integer j,k

      character*140 lat_file
      character*120 line, last_line, changes(1:3)/3*' '/
      character*16 element_name(2)/'B05E#1','B05E#2'/

      logical start_flag/.true./

      nargs = cesr_iargc()
      if(nargs == 1)then
         call cesr_getarg(1,lat_file)
         print *, 'Using ', trim(lat_file)
       else
      lat_file = 'bmad.'
      type '(a,$)',' Lattice file name ? (default= bmad.) '
      read(5,'(a)') line
       call string_trim(line, line, ix)
       lat_file = line
       if(ix == 0) lat_file = 'bmad.'
       type *, ' lat_file = ', lat_file
      endif

       call bmad_and_xsif_parser (lat_file, lat)

       print '(a,$)',' Lattice element ?'
       read *,element_name
       print '(a,$)',' x,y offset from nominal detector location ? '
       read *,xoff,yoff
       print '(a,$)',' x,xp at start of element ? '
       read *,orbit_x, orbit_xp

!      do j =-5,5
!       do k=-5,5

       x_image = x_image_0 + xoff
       y_image = y_image_0 + yoff

      tmin = twopi
     do i=1,2
       
       call element_locator(element_name(i),lat,ix)
         print '(i,1x,a,1x,i)', i,element_name(i), ix
       length(i) = lat%ele(ix)%value(l$)
    
       print '(a6,i,a16,f12.4)',' ix = ', ix, lat%ele(ix)%name, lat%ele(ix)%s
       orbit_x = (i-1)*orbit_xp * length(1) + orbit_x
       print '(2(a12,f12.4))',' orbit_x = ', orbit_x,' orbit_xp = ', orbit_xp 
       theta_0 = lat%ele(ix-1)%floor%theta + orbit_xp !note that if xp>0, the magnitude of theta_0 (which is always <0) decreases

      rho=lat%ele(ix)%value(rho$) 
      deltas = 0.001
      s=0
      r11 = cos(theta_0)
      r12 = -sin(theta_0)
      r21 = -r12
      r22 = r11

       dx0 = r12*orbit_x !if theta_0=-270, r12=-1
       dy0 = -r22*orbit_x !if theta_0=-180, r22=-1
 
       y0 = lat%ele(ix-1)%floor%x
       x0 = lat%ele(ix-1)%floor%z


!       print '(a6,f12.4,a6,f12.4,a10,f12.4)',' x0 = ',x0,' y0 = ', y0,' theta0 = ', theta_0
!       print '(5a12)','s','x','y','theta','theta_image' 
       do while(s < lat%ele(ix)%value(l$))
         s= s+deltas 
         theta = -s/rho
         dx = -(rho) * sin(theta)
         dy = -(rho) *(1- cos(theta))
         x = x0 + r11*dx + r12*dy + dx0
         y = y0 + r21*dx + r22*dy + dy0
         theta_image = atan2((y_image-y) ,(x_image-x))
!         if(theta + theta_0 > twopi/2) theta = theta-twopi
!        print '(5f12.4)', s,x,y,theta+theta_0, theta_image         
        if(abs(theta_image-(theta+theta_0+twopi)) < tmin)then
         tmin = abs(theta_image-(theta+theta_0+twopi))
         x_tan = x
         y_tan = y
         s_off = s
         s_tan = s+lat%ele(ix-1)%s
        endif
        write(11,'(6f12.4)') s_tan,x,y,theta+theta_0, theta_image, theta+theta_0+twopi
       end do
!       print '(3f12.4)', lat%ele(ix)%floor%z, lat%ele(ix)%floor%x,lat%ele(ix)%floor%theta
!       print '(3f12.4)', x,y,theta
   end do
       write(11,'(/,a,/)')"#"
       write(11,'(3f12.4)')x_image, y_image, s_off
       write(11,'(3f12.4)')x_tan,y_tan, s_tan
       write(6,'(/,a,/)')"#"
       write(6,'(3f12.4)')x_image, y_image, s_off
       write(6,'(3f12.4)')x_tan,y_tan, s_tan
        write(11,'(/,a,/)')"#"
      if(start_flag)then
      write(12, '(a)')'source location and image offset'
      write(12, '(5a12)')'s','x','y','x_image','y_image'
      start_flag = .false.
      endif
      write(12, '(5f12.4)')s_tan, x_tan, y_tan, j*xoff, k*yoff
      write(12,'(a15,2f12.4)')' inital x, xp = ',orbit_x,orbit_xp
!    end do
!   end do
end