program MuonDecay

 use bmad
 use random_mod
! use MuonDecayInterface

 implicit none
 integer :: seed = 0. 
 integer i
 real(rp) harvest
 real(rp) x,y,costheta
 real(rp) energy
 real(rp) x1, x2, tol

! 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)

  write(11,'(2es12.4)') y,costheta

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

end program