program MuonDecay

 use bmad
 use random_mod
! use MuonDecayInterface

 implicit none
 integer :: seed = 0. 
 integer i
 real(rp) harvest
 real(rp) x,y,costheta,t
 real(rp) energy
 real(rp) x1, x2,t2,tol
 real(rp):: y2 = 0.4
 real(rp):: err=0.02 ! the specific value of y when looking at angular distribution, and the allowed error
 real(rp)::lifetime = 2.2e-6,n0=1
! 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                                                                    
  call ran_seed_put(seed)

do i=1,100000

  call ran_uniform(harvest)
  x=harvest
  x1=0.
  x2=1.
  tol = 0.0001

  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
  t2=5*lifetime

  if (x /= n0) then
  t = -lifetime*log((n0-x)/n0)
  else 
  t=-1
  end if

  !if ( (y2-err <= y) .and. (y2+err >= y)) then  
  write(11,'(3es12.4)') y,costheta,t
  ! endif 

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 program