!+
! Subroutine track1_custom (start_orb, ele, param, end_orb, track, err_flag, finished)
!
! For scattering in injection line "markers". And checking special apertures like inflector "D"
!  And for recording electrons that intercept TRACKER planes (defined as markers) 
!
! Note: This routine is not to be confused with track1_custom2.
! See the Bmad manual for more details.
!
! General rule: Your code may NOT modify any argument that is not listed as
! an output agument below."
!
! Modules Needed:
!   use bmad
!
! Input:
!   start_orb  -- Coord_struct: Starting position.
!   ele    -- Ele_struct: Element.
!   param  -- lat_param_struct: Lattice parameters.
!
! Output:
!   end_orb   -- Coord_struct: End position.
!   track     -- track_struct, optional: Structure holding the track information if the 
!                 tracking method does tracking step-by-step.
!   err_flg  -- Logical: Set true if there is an error. False otherwise.
!-

subroutine track1_custom (start_orb, ele, param, end_orb, err_flag,  finished, track)

use bmad_struct, only:ele_struct, coord_struct, lat_param_struct, track_struct
use bmad_interface, dummy => track1_custom
use materials_mod
use parameters_bmad   !includes eloss, us_scatter, ds_scatter, tilt

implicit none

type (coord_struct) :: start_orb
type (coord_struct) :: end_orb
type (coord_struct) temp_orb, delta,temp_orb2
type (ele_struct) :: ele
type (lat_param_struct) :: param
type (track_struct), optional :: track
logical hit/.false./,hit_y/.false./
logical err_flag
logical energy_loss
logical withinInflectorAperture
logical first/.true./
!logical first1/.true./,first2/.true./,first3/.true./
!logical fbfirstx1/.true./,fbfirstx2/.true./,first_plane_detector/.true./
logical :: first_calo = .true., first_harp = .true.
logical finished
!logical first_muonx1/.true./,first_muony1/.true./,first_muonx2/.true./,first_muony2/.true./

real(rp) limit_in/6.905/,limit_out/7.010/   !radial limits for the trackers
real(rp), save :: t_first,t_first2,t_first3,t_first4 !first time a muon hit a fiber monitor
real(rp) radius/7.112/                      !radius of the ring
real(rp) apRadius/0.045/                    !radius of apeture
real(rp) fbdia/0.0005/,spacing/0.013/       !diameter and length of intervals of the fibers
real(rp),dimension(1:7):: xfiberPosition   !positions of fibers measuring x
real(rp)::tcycle = 2*Pi*rMagic/(0.9994*speedoflight) !time for muon to travel for one turn
real(rp) thickness

character(32) :: r_name = 'track1_custom'
character(26) :: filename
character(3)  :: turn
integer lun,k,j
integer nturn/128/
integer, save:: lun51,lun52

integer caloNum ! calorimeter number (1 to 24)
integer harpNum ! harp monitor number (1 to 4)
integer muonNum
integer fiber
real(rp) :: x_hit, y_hit ! coords where particle hits calorimeter
real(rp), parameter :: calo_edge_x1 = -0.31481 ! coords to define edges of calo front face
real(rp), parameter :: calo_edge_x2 = -0.08981
real(rp), parameter :: calo_edge_y1 = -0.075
real(rp), parameter :: calo_edge_y2 =  0.075

!

!call out_io (s_fatal$, r_name, 'THIS DUMMY ROUTINE SHOULD NOT HAVE BEEN CALLED IN THE FIRST PLACE.')
err_flag = .false.

! Remember to also set end_orb%t

withinInflectorAperture = .true.
energy_loss = eloss 
end_orb = start_orb
temp_orb = start_orb

if(first)then
  lun51=lunget()
  open(unit=lun51,file = trim(directory)//'/'//'UpstreamInflectorScattering.dat') 
  lun52=lunget()
  open(unit=lun52,file = trim(directory)//'/'//'DownstreamInflectorScattering.dat') 
  write(lun51,'(a)')'Change in coordinate vector due to upstream inflector scattering'
  write(lun51,'(6a12)')'x','xp','y','yp','z','zp'
  write(lun52,'(a)')'Change in coordinate vector due to downstream inflector scattering'
  write(lun52,'(6a12)')'x','xp','y','yp','z','zp'

 lun=lunget()
 open (unit=lun, file = trim(directory)//'/'//"log.dat",access='append')
 write(lun, '(a)')'TRACK1_CUSTOM: Scattering in inflector ends'
 write(lun, '(1x,a15,l)')' energy_loss = ', energy_loss
 write(lun,'(1x,a16)')ele%name
 first=.false.
 close(unit=lun)
endif
!print '(a,1x,2(a,2es12.4),2es12.4)',ele%name,'start_orb = ',start_orb%vec(1:2), ' end_orb = ', end_orb%vec(1:2), start_orb%s, ele%s

if(ele%name == 'MARK_CRYO_US' .and. us_scatter)then
      call offset_particle(ele, param, set$, temp_orb)
      call scatter( Al,                0.001_rp, temp_orb)
      call offset_particle(ele, param, unset$, temp_orb)
      end_orb = temp_orb
      end_orb%s = ele%s 
   return
endif
if(ele%name == 'MARK_INFLECTOR_US')then
     withinInflectorAperture = .false.
!     print '(a2,1x,6es12.4)','bs',temp_orb%vec(1:6)
     if(us_scatter)then
      call offset_particle(ele, param, set$, temp_orb)
      call inflector_scatter1(temp_orb, .true., .false., energy_loss)
      call offset_particle(ele, param, unset$, temp_orb)
     endif
     delta%vec = temp_orb%vec - start_orb%vec
     write(lun51,'(6es12.4)')delta%vec(1:6)
!     print '(a2,1x,12es12.4,a,es12.4)','us',temp_orb%vec(1:6),delta%vec(1:6), ' inflector_width= ', inflector_width
!    if(inflector_width  <= 0.00901)then
!      call E821InflectorAperture(temp_orb%vec(1),temp_orb%vec(3),withinInflectorAperture,.true.,.false.)
!     else
!      call E989InflectorAperture(temp_orb%vec(1),temp_orb%vec(3),withinInflectorAperture,.true.,.false.)
!    endif
!    if(.not. withinInflectorAperture)temp_orb%state=lost$
    end_orb = temp_orb
    end_orb%s = ele%s 
 return
endif
if(ele%name == 'MARK_INFLECTOR_DS')then
    withinInflectorAperture = .false.
    if(ds_scatter)then
      call offset_particle(ele, param, set$, temp_orb)
      call inflector_scatter1(temp_orb, .false., .true., energy_loss)
      call offset_particle(ele, param, set$, temp_orb)
     endif
     delta%vec = temp_orb%vec - start_orb%vec

     write(lun52,'(6es12.4)')delta%vec(1:6)
!     print '(a2,1x,6es12.4)','ds',delta%vec(1:6)
!    if(inflector_width <= 0.00901)then
!      call E821InflectorAperture(temp_orb%vec(1),temp_orb%vec(3),withinInflectorAperture, .false.,.true.)
!     else
!      call E989InflectorAperture(temp_orb%vec(1),temp_orb%vec(3),withinInflectorAperture, .false.,.true.)
!    endif
!    if(.not. withinInflectorAperture)temp_orb%state = lost$
    end_orb = temp_orb
    end_orb%s = ele%s 
 return
endif



!==========================================
! fiber monitor section
! added by dlr, September 2016
! revised for multiple scattering January 2017 
!==========================================

if (index(ele%name, 'FIBER')/=0) then

  ! check coords to see if particle hit calorimeter front face
  x_hit = start_orb%vec(1)
  y_hit = start_orb%vec(3)

    if (first_harp) then ! create file and write header info
      lun = lunget()
      open(unit=lun, file=trim(directory)//'/'//"harp_plane_hits.dat", status="new", action="write")
      write(unit=lun, fmt='(a6,8a12)') " harp", &
                                           "t         ","s         ", &
                                           "x         ","px        ", &
                                           "y         ","py        ", &
                                           "z         ","pz        "
      close(unit=lun)
      open(unit=lun, file=trim(directory)//'/'//"fiber_hits.dat", status="new", action="write")
      write(unit=lun, fmt='(a6,1x,a6,9a12)') " harp","fiber", &
                                           "t         ","s         ", &
                                           "x         ","px        ", &
                                           "y         ","py        ", &
                                           "z         ","pz        ",'thickness'
      close(unit=lun)
      first_harp = .false.
    end if

    ! get harp num (integer) from element name
    read (ele%name(6:6),'(i1)') harpNum
    call which_fiber(harpNum, start_orb, fiber, thickness)

    ! write information about this hit to the file
    lun = lunget()
    open(unit=lun, file=trim(directory)//'/'//"harp_plane_hits.dat", status="old", access="append", action="write")
    write(unit=lun, fmt='(i5,8es18.8)')     harpNum, & 
                                            start_orb%t, start_orb%s, &
                                            start_orb%vec(1), start_orb%vec(2), &
                                            start_orb%vec(3), start_orb%vec(4), &
                                            start_orb%vec(5), start_orb%vec(6)
    close(unit=lun)

    if(fiber /= 0)then
       open(unit=lun, file=trim(directory)//'/'//"fiber_hits.dat", status="old", access="append", action="write")
       write(unit=lun, fmt='(i5,1x,i5,9es18.8)')     harpNum, fiber, & 
                                            start_orb%t, start_orb%s, &
                                            start_orb%vec(1), start_orb%vec(2), &
                                            start_orb%vec(3), start_orb%vec(4), &
                                            start_orb%vec(5), start_orb%vec(6), thickness
       close(unit=lun)
       end_orb = start_orb
       if(FiberScatter)then
          call scatter(SciFiber, thickness, end_orb)
          call EnergyLoss(SciFiber,thickness,end_orb)
       endif
       write(111,'(7es12.4)')thickness,start_orb%vec - end_orb%vec
    endif ! fiber hit
 return
end if ! harp plane


!========================
! end fiber monitor section
!========================
!==========================================
! calorimeter section
! added by Robin Bjorkquist, September 2015 
!==========================================

if (start_orb%species == positron$) then
if (index(ele%name, 'CALO')/=0) then

  ! check coords to see if particle hit calorimeter front face
  x_hit = start_orb%vec(1)
  y_hit = start_orb%vec(3)
  if (calo_edge_x1 < x_hit .and. x_hit < calo_edge_x2 .and. calo_edge_y1 < y_hit .and. y_hit < calo_edge_y2) then

    if (first_calo) then ! create file and write header info
      lun = lunget()
      open(unit=lun, file=trim(directory)//'/'//"calo_hits.dat", status="new", action="write")
      write(unit=lun, fmt='(a10,a5,8a12)') "  muon_num"," calo", &
                                           "t         ","s         ", &
                                           "x         ","px        ", &
                                           "y         ","py        ", &
                                           "z         ","pz        "
      close(unit=lun)
      first_calo = .false.
    end if

    ! get calo num (integer) from element name
    read (ele%name(5:6),'(i2)') caloNum

    ! get muon num, which is stored in the path_len field of the positron's coord_struct
    muonNum = int(start_orb%path_len)

    ! write information about this hit to the file
    lun = lunget()
    open(unit=lun, file=trim(directory)//'/'//"calo_hits.dat", status="old", access="append", action="write")
    write(unit=lun, fmt='(i10,i5,8es12.4)') muonNum, &
                                            caloNum, & 
                                            start_orb%t, start_orb%s, &
                                            start_orb%vec(1), start_orb%vec(2), &
                                            start_orb%vec(3), start_orb%vec(4), &
                                            start_orb%vec(5), start_orb%vec(6)
    close(unit=lun)

    ! kill positron track
    if(start_orb%species == positron$)temp_orb%state = lost$

  end if ! particle hit front face
end if ! CALO
end_orb = temp_orb
end_orb%s = ele%s
return
end if ! positron

!========================
! end calorimeter section
!========================

!if(end_orb%state == lost$)print '(a12,1x,a16,7es12.4)','Lost in',ele%name,end_orb%s,end_orb%vec(1:6)

end subroutine