subroutine muon_pulse_shape(pulse, time_out)
use bmad
use parameters_bmad
implicit none

character*120 line
!character*18 filename(9)/'pulseshape_0.dat','pulseshape_1.dat','pulseshape_2.dat','pulseshape_3.dat','pulseshape_4.dat','pulseshape_5.dat','pulseshape_6.dat','pulseshape_7.dat', &
!                           'pulseshape_ave.dat'/    ! january 2018
character*45 filename(10)/'pulseshape_0_12543_mod.dat','pulseshape_1_12543_mod.dat','pulseshape_2_12543_mod.dat','pulseshape_3_12543_mod.dat','pulseshape_4_12543_mod.dat', &
                                          'pulseshape_5_12543_mod.dat','pulseshape_6_12543_mod.dat','pulseshape_7_12543_mod.dat',  'averageAllPulses_12543.dat', 'avgBeamPulsePDF.txt'/    ! january 2018
character*95 pulse_file
integer i,  pulse
integer, save :: i_tot, center, mina, minb
integer, save :: lun
integer, save :: save_pulse
integer ios
real(rp), save :: time(10000), a_pulse(10000), b_pulse(10000), int_apulse(10000), int_bpulse(10000)
real(rp), save :: int_pulse(10000)
real(rp)  time_out
real(rp), save ::   width
real(rp) di, frac_low, frac_high, i_interp
real(rp) rand
logical first/.true./
logical itexists

if(first .or. abs(pulse)/= save_pulse )then
  save_pulse = abs (pulse)
  lun =     lunget()
  pulse_file ='pulse_shapes/'//trim(filename(save_pulse))  
print *,save_pulse, pulse_file
 inquire (file=pulse_file, exist=itexists)
  if(.not. itexists)then
    print '(a,a,a)',' muon pulse shape file ', pulse_file,' is not found'
    stop
  endif
  open(unit=lun, file = pulse_file, status='old')
  print '(a,a)',' open file = ', pulse_file
 I=0
 a_pulse(:)=0
 b_pulse(:)=0
 time(:)=0
 do while(.true.)
  read(lun,'(a)', IOSTAT=ios)line
  if(ios  <0)exit
  if(index(line,'trace')/= 0 .or. index(line,'#')/=0 .or. line(1:10) == '          ' .or. index(line,'NaN')/=0 )cycle
   i=i+1
   read(line,*)time(i), a_pulse(i) !, b_pulse(i)
!   print *,line
!   print *, time(i), a_pulse(i) !, b_pulse(i)
   b_pulse(i) = a_pulse(i)   !this is temporary until b_pulse is made available in the pulse_shape file
   i_tot = i
   if(i_tot > 10000)then
      print *,' muon_pulse_shape: i_tot  10000 line 53'
      stop
   endif
   int_apulse(i) = 0
   int_bpulse(i) = 0
  end do
  close(unit=lun)
   mina = minval(a_pulse)
   minb = minval(b_pulse)
! integrate
  do i=1, i_tot
    int_apulse(i) = sum(a_pulse(1:i))-i*mina 
    int_bpulse(i) = sum(b_pulse(1:i))-i*minb
   end do
!normalize
   int_apulse(:) = int_apulse(:)/int_apulse(i_tot)
   int_bpulse(:) = int_bpulse(:)/int_bpulse(i_tot)
  if(pulse >0)then
     int_pulse(:) = int_apulse(:)
     center = maxloc((a_pulse(1:i_tot)),1)
  endif
  if(pulse <0)then 
     int_pulse(:) = int_bpulse(:)
     center = maxloc((b_pulse(1:i_tot)),1)
  endif
  lun = lunget()
  open(unit=lun, file = trim(directory)//'/muon_pulse_shape.dat')
  write(lun,'(a1,1x,a,a,1x,i12)')'#', pulse_file,' center = ',center
  do i=1,i_tot
  write(lun,'(3es12.4)')time(i), a_pulse(i), int_pulse(i)
  end do
  width = time(i_tot) - time(1)
  print *,'muon_pulse_shape: center =',center,' width = ',width, ' time(center) = ', time(center)

 first = .false.
endif !initialization
 
  ! Generate a random number in [0,1)
  call ran_uniform(rand)
!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
     IF (int_pulse(i)<=rand .and. rand<int_pulse(i+1)) THEN
      di = int_pulse(i+1) - int_pulse(i)        ! values in the int_apulse table (y-axis)
      frac_low  = ( int_pulse(i+1) - rand )/di       ! fraction of lower value to include (linear interpolation)
      frac_high = ( rand -  int_pulse(i)  )/di       ! fraction of upper value to include (linear interpolation)
      i_interp  = frac_low*i + frac_high*(i+1)            ! corresponding fractional index (x-axis)
      time_out  = (width/2. - center + (i_interp/i_tot - 0.5) * width)*1.e-9 ! map to "center" and "width"
      time_out  = (- center/i_tot*width + (i_interp/i_tot - 0.5) * width)*1.e-9 ! map to "center" and "width"
      EXIT
    ENDIF
  end do

return
end