COSMOS v7.655  COSMOSv7655
(AirShowerMC)
Functions/Subroutines
kbetar.f File Reference

Go to the source code of this file.

Functions/Subroutines

subroutine kbetar (a, b, x)
 

Function/Subroutine Documentation

◆ kbetar()

subroutine kbetar ( real*8  a,
real*8  b,
real*8  x 
)

Definition at line 19 of file kbetar.f.

References rndc(), and true.

Referenced by cprtsampe().

19 ! **************************
20  implicit none
21  real*8 a ! input. see below
22  real*8 b ! input. see below
23  real*8 x ! output. sampled random variable
24 
25 ! generate a random variable x following beta function kernel
26 ! i.e. density function is x**(a-1)* (1-x)**(b-1)dx
27 !
28 ! both of a and b cannot be < 1.0
29 ! if a and b unbalance largely, rejection rate will
30 ! become large: (a=2, b=5--->1/6.6 is accepted
31 ! a=2.1,b=1.2-->30 % rejection)
32 ! To sample 100000 variables, 4.226 sec is needed.
33 ! (wiht a=1.5, b=2.5 with facom m780)
34 ! this is sometimes faster than ksbeta for a class of a,b
35 ! but has limitations.
36 
37  real*8 ai, bm, u, bi, am
38 
39  if(a .lt. 1. and. b .lt. 1.) then
40  write(*,*) ' kbetar cannot accept a<1,b<1:'
41  write(*,*) ' use, ksbeta in such case'
42  write(*,*) ' a=',a, ' b=',b
43  stop
44  endif
45  if( (a .ge. b .and. b .ge. 1.) .or. a .lt. 1.) then
46  ai=1./a
47 
48  bm=b-1.
49 ! *** until loop***
50  do while (.true.)
51  call rndc(u)
52 
53  x=u**ai
54  call rndc(u)
55  if ( u .le. (1.-x)**bm)
56  * goto 100
57  enddo
58  100 continue
59  else
60  bi=1./b
61 
62  am = a-1.
63 ! *** until loop***
64  do while (.true.)
65  call rndc(u)
66  x=1.-u**bi
67  call rndc(u)
68  if ( u .le. x**am)
69  * goto 200
70  enddo
71  200 continue
72  endif
true
block data cblkElemag data *AnihiE ! Eposi< 1 TeV, anihilation considered *X0/365.667/, ! radiation length of air in kg/m2 *Ecrit/81.e-3/, ! critical energy of air in GeV *MaxComptonE/1./, ! compton is considered below 1 GeV *MaxPhotoE/1.e-3/, ! above this, PhotoElectric effect neg. *MinPhotoProdE/153.e-3/, ! below 153 MeV, no gp --> hadrons ! scattering const not MeV *Knockon true
Definition: cblkElemag.h:7
rndc
subroutine rndc(u)
Definition: rnd.f:91
nrl_atmos::a
real(4), save a
Definition: cNRLAtmos.f:20
nrl_atmos::b
real(4), save b
Definition: cNRLAtmos.f:21
x
! structure defining a particle at production ! Basic idea of what is to be contained in ! the particle structue is that dynamical ones should be included those derivable from the particle code ! is not included ******************************************************type fmom momentum sequence union map real e endmap map real * x
Definition: Zptcl.h:21
Here is the call graph for this function:
Here is the caller graph for this function: