!+
! Subroutine wall_hit_handler_custom (orb, ele, s, t)
!
! Dummy routine.
! This routine is called by the Runge-Kutta integrator odeint_bmad when a particle hits a wall.
! This routine can be replaced by a custom routine to do custom calculations.
!
! Modules needed:
!   use bmad
!
! Input:
!   orb   -- coord_struct: coordinates of particle.
!   ele   -- ele_struct: Element holding the aperture
!   s     -- real(rp): Longitudinal position from start of element.
!
! Output:
!   Any argument may be modified...
!-

subroutine wall_hit_handler_custom (orb, ele, s)

!use bmad_interface
use materials_mod  !contains quad_plate
use wall3d_mod ! , dummy => wall_hit_handler_custom

implicit none

type (coord_struct) :: orb, temp, orb_in
type (ele_struct) :: ele

real(rp) :: s
real(rp) height/0.0235/
real(rp) effective_thickness, thickness
real(rp) perp(3), d_radius
real(rp) vec(3), position(6)
real(rp) cost

character(*), parameter :: r_name = 'wall_hit_handler_custom'

!WRITE(*,'(a,a,3(a4,es10.3,4x))')'WALL_HIT_HANDLER_CUSTOM:',ele%name, 'x = ', orb%vec(1), 'y = ', orb%vec(3), 's = ', s
!print *,' WALL HIT HANDLER CUSTOM in: orb%state = ', coord_state_name(orb%state)

  if(.not. associated(ele%wall3d))return
  if(ele%aperture_at /= wall_transition$)return
  if (index(ele%name,'QUAD')==0)return
  
  orb%state=alive$
  if(.not. quad_plate)return
  if(orb%species == positron$) return
  if(abs(orb%vec(3)) > height .and. abs(orb%vec(1)) >height)return   ! misses the plates
  orb%state=lost$ 

    temp = orb
    position = [orb%vec(1:4),orb%s-ele%s_start,1.0_rp]
    thickness =ele%wall3d(1)%thickness
    vec(1:3) = (/orb%vec(2), orb%vec(4), 1._rp/)

    d_radius = wall3d_d_radius (position,ele, perp= perp)    
!    if((perp(1) < 0 .and. abs(perp(1)) < 0.1) .or. abs(perp(1)) < 0.1)return 
    cost = dot_product(vec,perp)/sqrt(dot_product(vec,vec)*dot_product(perp,perp))
    effective_thickness = abs(ele%wall3d(1)%thickness/cost)
! compute scattering angle only
       orb_in%vec=0
       call scatter(Al,effective_thickness,orb_in)
       orb%vec(2) = orb%vec(2) + orb_in%vec(2)
       orb%vec(4) = orb%vec(4) + orb_in%vec(4)
!       write(53,'(a,2es12.4,11es12.4)')'wall_hit_handler: thickness, effective_thickness, delta%vec, perp = ' &
!                                 ,thickness,effective_thickness,temp%vec(1:6)-orb%vec(1:6), orb%vec(1),orb%vec(3),perp

       orb%state = alive$
!print *,' WALL HIT HANDLER CUSTOM out: orb%state = ', coord_state_name(orb%state)
end subroutine wall_hit_handler_custom