subroutine compute_decayTime(muons,nmuons,unit)

 use bmad
 use muon_mod
 use random_mod
! use MuonDecayInterface

 implicit none
 integer :: seed = 0. 
 integer i,nmuons,unit
 real(rp) harvest
 real(rp) x,y,costheta,t
 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./
! real(rp) MuonToElectronEnergy,  MuonDecayAngDist
 

!generate distribution of electrons assuming muon at rest and polarization P in the +z direction

 !initialize randomizer's seed (harvest)                                                                                    
  ! 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(harvest)
 ! x=harvest
 ! x1=0.
 ! x2=1.
 ! tol = 0.0001             !precision

 !  y =   inverse(MuonToElectronEnergy,x, x1,x2,tol)

 ! call ran_uniform(harvest)

 ! x=harvest

 ! costheta = inverse(MuonDecayAngDist,x,-x2,x2,tol)

  call ran_uniform(harvest)
  
  x=harvest
 
  if (x /= n0) then
  t = -lifetime*log((n0-x)/n0)           !invert the distribution
  else 
  t= 0                                   !t = infinity treated as t = 0 
  end if

muons(i)%coord%path_len = t              !use the path_length attribute to store decay time

write(unit,'(es12.4)') t 

end do

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