subroutine fiber_energy_loss(thickness,co)
use bmad
use parameters_bmad
implicit none

type (coord_struct) co
character*120 line
character*18 filename/'eloss.dat'/

integer i
integer, save :: i_tot, center, mina, minb
integer lun
integer ios
integer, save :: prob(1000)
real(rp), save :: energy(1000), int_prob(1000)
real(rp) thickness, ave_thickness
real(rp)  eloss_out, eloss_ave
real(rp), save ::   width, save_loss
real(rp) di, frac_low, frac_high, i_interp
real(rp) rand
real(rp) fiber_radius/0.00025/
real(rp) :: p, betagamma, rapidity,E
real(rp) arbitrary_scaling_factor/1.2/ !so that change of average energy with time looks like harp measurement

logical first/.true./

if(first )then
  print *,'fiber_energy_loss'
  lun =     lunget()
  open(unit=lun, file = filename)
  print '(a,a)',' open file = ', filename
 I=0
 prob(:)=0
 energy(:)=0
 do while(.true.)
  read(lun,'(a)', IOSTAT=ios)line
  if(ios  <0)exit
  if(index(line,'eloss')/= 0)cycle
   i=i+1
   read(line,*)energy(i), prob(i)
   i_tot = i
   int_prob(i) = 0
  end do
   mina = minval(prob)
! integrate
  do i=1, i_tot
    int_prob(i) = sum(prob(1:i))-i*mina 
   end do
!normalize
   int_prob(:) = int_prob(:)/int_prob(i_tot)
!   int_bpulse(:) = int_bpulse(:)/int_bpulse(i_tot)
   center = maxloc((prob(1:i_tot)),1)
  do i=1,i_tot
  write(14,'(i10,1x,i12)')i, prob(i)
  write(13,'(i10,es12.4)')i,int_prob(i)
  end do
  width = energy(i_tot) - energy(1)
  ave_thickness = 8/twopi*fiber_radius
  write(14,'(a,1x,i12,es12.4,1x,i12,1x,i12)')'center, width, i_tot', center,width,i_tot,mina
 first = .false.
endif !initialization
  ! Energy loss depends on energy/momentum
  p         = ( 1+co%vec(6) )*pMagic  ! momentum magnitude (eV)
  betagamma = p/mmu                   ! Lorentz betagamma
  rapidity  = asinh(betagamma)        ! hyperbolic rotation angle
  E         = cosh(rapidity)*mmu      ! energy (eV)
 


  ! Generate a random number in [0,1)
  rand=-1
  do while(rand < int_prob(1) .or. rand > int_prob(i_tot))
  call ran_uniform(rand)
  end do
!print *,rand
  ! Find the random number in the pulse integral table, and convert to a number based on the user's center/width arguments
  do i = 1, i_tot-1
     IF (int_prob(i)<=rand .and. rand<int_prob(i+1)) THEN
      di = int_prob(i+1) - int_prob(i)        ! values in the int_apulse table (y-axis)
      frac_low  = ( int_prob(i+1) - rand )/di       ! fraction of lower value to include (linear interpolation)
      frac_high = ( rand -  int_prob(i)  )/di       ! fraction of upper value to include (linear interpolation)
      i_interp  = frac_low*i + frac_high*(i+1)            ! corresponding fractional index (x-axis)
      eloss_ave  =  i_interp/i_tot  * width ! map to "center" and "width"
      eloss_out = eloss_ave !* thickness/ave_thickness
!      write(112,'(3es12.4)')eloss_out, i_interp, i_interp/i_tot*width
      EXIT
    ENDIF
  end do

  ! Change the energy
  E = E  - eloss_out*1.e6 * arbitrary_scaling_factor

  ! Calculate the new momentum from the new energy
  betagamma = sinh(acosh(E/mmu))

  ! Convert to a momentum shift
  co%vec(6) = mmu*betagamma/pMagic - 1.
!write(112,'(3es12.4)')E,betagamma,co%vec(6)
return
end