!+
!subroutine compute_decayTime(muons,nmuons,unit)
!
!compute and assign to each muon its decay time according to the decay time distribution.
!
!Note:the time is stored in the muons(i)%coord%r attribute.
!
!Modules Needed:
!
! use bmad
! use muon_mod
! use random_mod
!
!Input:
!      muons -- muon_struct:the muon array
!      nmuons -- integer:number of muons in the array
!      unit   -- integer:unit to record the times
!
!Output:
!      muons -- muon_struct:the muon array with decay time attached
!- 

subroutine compute_decayTime(muons,nmuons,unit)

 use bmad
 use muon_mod
 use random_mod
use nr
! use MuonDecayInterface

 implicit none
 integer :: seed = 0. 
 integer i,nmuons,unit
 real(rp) gamma/29.3/
 real(rp) temp, gam
! real(rp), allocatable :: harvest(:)
 real(rp) harvest, dum
 real(rp) x,y,costheta,t,z
 real(rp) energy
 real(rp) x1,x2,tol
 type (muon_struct), dimension(nmuons):: muons
 real(rp)::lifetime = 2.2e-6,n0=1
 logical  first/.true./

 ! ran_seed = 0 <--> seed randomizer w/ current CPU time 

write(unit,'(a20)')    'Muon Decay Time'
                                                                   
call ran_seed_put(seed)

do i=1,nmuons

call ran_uniform(dum)

  t = -lifetime*log(dum)

temp = 1.-muons(i)%coord%beta **2
if(temp > 0.)then
  gam = 1./sqrt(temp)
 else
  print '(a,es12.4,a,a)',' compute_decayTime: 1-beta**2 = ', temp,' which is less than zero ',' set to default'
  gam = gamma
endif

muons(i)%coord%r = t*gam

!write(199,'(i8,2es12.4)') i,muons(i)%coord%r,dum

end do
if(nmuons == 1)muons(1)%coord%r = 6.8e-7   !4 turns
close(unit = unit)
contains

!+
! Function MuonToElectronEnergy(x)
!
! Routine to compute the probability that the muon decay electron will have fraction energy x=E/Emax
!
! Defining relation:
!     MuonToElectronEnergy(x) = x**3*(2-x)
!

function MuonToElectronEnergy(x)

  use precision_def

  implicit none

  real(rp) MuonToElectronEnergy, x

  MuonToElectronEnergy = x**3 * (2-x)
  
  energy = x

end function

!+
! Function MuonDecayAngDist(x,y)
!
! Routine to compute the probability that the decay electron momenum is at and angle acos(x) with respect to muon spin
!
! Defining relation 
!    MuonDecayAngDist(x) = 0.5*(1-x)*(1-a-a*x)
!    where a = 0.5*(2*y-1)/(3-2*y)
!

function MuonDecayAngDist(x)

   use precision_def

   implicit none

   real(rp) MuonDecayAngDist, x,y
   real(rp) a

   a=0.5*(2*energy-1)/(3-2*energy)
   MuonDecayAngDist = 0.5*(1-x)*(1-a-a*x)

end function

!+
! Function MuonDecayTDist(t)
!
! Routine to compute the portion of decayed muons at time t

function MuonDecayTDist(t)

   use precision_def

   implicit none

   real(rp) MuonDecayTDist,t
   real(rp)::lifetime = 2.2e-6,N0=1

  MuonDecayTDist=N0-N0*exp(-t/lifetime)

end function

end subroutine