COSMOS v7.655  COSMOSv7655
(AirShowerMC)
ksgamd.f
Go to the documentation of this file.
1 !c test ksgamd
2 ! include 'kgauss.f'
3 ! include 'kgamma.f'
4 ! include 'rnd.f'
5 !--------------------------------
6 ! implicit none
7 ! real*8 alfa, am, xm, x
8 ! integer i
9 !
10 ! alfa=15.
11 ! am=105.
12 ! xm=alfa/am/(am-alfa)
13 ! do i=1, 10000
14 !c activate one of the following calls.
15 !c call ksgmim(3, 300.d0, x)
16 !c call ksgmis(2, 10.d0, x)
17 !c call ksgmrm(3.5d0, 100.d0, x)
18 !c call ksgmrs(3.5d0, 20.d0, x)
19 !c call ksgmrs(0.2d0, 30.d0, x)
20 !c call ksgmrs(-0.45d0, 300.d0, x)
21 !c call ksgmrs(-0.95d0, 20.d0, x)
22 !c call ksgmrm(1.00d0, 2.d0, x)
23 ! write(*, *) sngl(x)
24 ! enddo
25 ! end
26 ! *****************************************************************
27 ! * sampling for the gamma distribution: (note power is not s-1)
28 ! *
29 ! * (x/a)**s exp(-x/a)/gamma(s+1)d(x/a).
30 ! *
31 ! * (s>-1.0; mean is a*(s+1)).
32 ! *
33 ! * ksgmim: sampling for an integer s with a given mean
34 ! * ksgmis: sampling for an integer s with a given scale parameter a
35 ! * ksgmrm: sampling for a real s with a given mean
36 ! * ksgmrs: sampling for a real s with a given scale parameter a
37 ! *
38 ! ****************************************************************
39 !
40 ! samples a random variable from gamma distribution, x**(s) exp(-x/a)
41 !
42 ! /usage/
43 ! call ksgmim(n, av, x)
44 ! n: input. integer>= 0. s=n.
45 ! av: input. real*8> 0. the average of the distribution
46 ! x: output. real*8. sampled random variable.
47 !
48 ! call ksgmis(n, a, x)
49 ! n: input. integer>=0. s=n
50 ! a: input. real*8. a in exp(-x/a)
51 ! the mean is given by a*(n+1)
52 ! x: output. real*8. sampled random variable.
53 !
54 ! call ksgmrm(s, av, x)
55 ! s: input. real*8 >-1.0
56 ! av: input. real*8. >0. the average of the distribution.
57 ! x: output. real*8. a sampled radonm varible.
58 ! if s is close to -1.0, x < 1.d-30 may not
59 ! be very accurate.
60 !
61 ! call ksgmrs(s, a, x)
62 ! s: input. real*8 > -1.
63 ! a: input. real*8 > 0. a in exp(-x/a)
64 ! the mean is given by a*(s+1)
65 ! x: output. real*8. sampled random variable.
66 ! if s is close to -1.0, x < 1.d-30 may not
67 ! be very accurate.
68 !
69 ! for s=3.5, 1653 msec is needed for calling 50000 times.
70 ! (by m780 35mips machine)
71 !
72  subroutine ksgmrs(s, a, x)
73 ! implicit none
74  real*8 s, a, x
75 !
76  real*8 alfa, u, r
77  integer n
78  logical more
79 !
80  if(s .le. -1.) then
81  write(*,*) ' input value to ksgmrs. s=',
82  * s,' invalid'
83  stop 9999
84  elseif(s .lt. 0.) then
85  call ksgamn(s, x)
86  else
87 ! sample a varible from x**s exp(-x)dx
88 ! decompose s=n+alfa, 1.0>alfa>=0.
89  n=s
90  alfa=s-n
91  if(alfa .eq. 0.) then
92 ! use integer case
93  call ksgmis(n, 1.d0,x)
94  else
95 ! importance sampling with weight function
96 ! rho= (1-alfa)*x**n exp(-x)/gamma(n+1) +
97 ! alfa*x**(n+1) exp(-x)/gamma(n+2)
98 ! =x**n exp(-x)/gamma(n+1) (1-alfa+alfa*x/(n+1))
99 ! then,
100 ! x**(n+alfa)exp(-x)/rho=
101 ! x**alfa/(1-alfa+alfa*x/(n+1))
102 ! which has a max value, (n+1)**alfa at x=n+1.
103 ! determine if x**n exp(-x) or x**(n+1) exp(-x)
104  more=.true.
105  do 100 while (more)
106 ! average trial # is 1.03 with max of 4 (for
107 ! s=3.5, 50000 cases).
108  call rndc(u)
109  if(u .le. alfa) then
110  call ksgmis(n+1,1.d0, x)
111  else
112  call ksgmis(n, 1.d0, x)
113  endif
114 ! rejection func.value
115  r= (x/(n+1))**alfa/(1.-alfa + alfa*x/(n+1))
116  call rndc(u)
117  more=u .gt. r
118  100 continue
119  endif
120  endif
121  x=a*x
122  end
123 !
124  subroutine ksgmrm(s, av, x)
125 ! implicit none
126  real*8 s, av, x
127 !
128  real*8 a
129 !
130  if(s .le. -1.) then
131  write(*,*) ' input value to ksgmrs. s=',
132  * s,' invalid'
133  stop 9999
134  endif
135  a=av/(s+1.)
136  call ksgmrs(s, a, x)
137  end
138 !
139  subroutine ksgmis(n, a, x)
140 ! implicit none
141  real*8 a, x
142  integer n
143 !
144  real*8 u, ui
145  integer i
146  logical more
147 !
148  if(n .le. -1) then
149  write(*,*) ' input value to ksgmis. n=',
150  * n,' invalid'
151  stop 9999
152  endif
153 !
154  more=.true.
155  do 100 while (more)
156  call rndc(u)
157  do 50 i=1, n
158  call rndc(ui)
159  u=u*ui
160  50 continue
161  more=u .le. 0.
162  100 continue
163  x=-a*log(u)
164  end
165 !
166  subroutine ksgmim(n, av, x)
167 ! implicit none
168  real*8 av, x
169  integer n
170 !
171  real*8 a
172 !
173  if(n .le. -1) then
174  write(*,*) ' input value to ksgmim. n=',n,' invalid'
175  stop 9999
176  endif
177 !
178  a=av/(n+1)
179  call ksgmis(n, a, x)
180  end
181 ! for -1.0<s<0.
182  subroutine ksgamn(s, x)
183 ! implicit none
184  real*8 s, x
185 !
186 !
187  real*8 pi/3.14159265/, rt2i/0.70710678/, eps/0.1d0/,
188  * err/1.d-3/, small/1.d-30/
189  real*8 u, r, us, br, kgamma, tmp, xold, acc
190  logical more, more2
191 !
192  if(s .ge. -0.5) then
193 ! importance sampling is completely possible by using weight
194 ! function of
195 ! rho=-2s* x**(-1/2)exp(-x)/gamma(1/2)+(2s+1)exp(-x).
196 ! then, x**s exp(-x)/rho=
197 ! x**s/( -2s x**(-1/2)/gamma(1/2) + (2s+1)).
198 ! this takes the
199 ! maximum, (1/gamma(1/2))**2s=pi**(-s) at x=(1/gamma(1/2))**2=1/pi
200 ! sampling for exp(-x)/root(x)dx is done by taking the square of
201 ! a gaussian random varible with mean 0 and variance 1/root(2).
202 !
203  more=.true.
204  do 100 while ( more )
205  call rndc(u)
206  if(u .le. (2*s+1.)) then
207  call rndc(u)
208  x=-log(u)
209  else
210  call kgauss(0.d0, rt2i, x)
211  x=x*x
212  endif
213 ! rejection function value
214  r=(pi*x)**(s+0.5)/(-2*s + (2*s+1.)*sqrt(pi*x))
215  call rndc(u)
216  more=u .gt. r
217  100 continue
218  else
219 ! complete rejection method not known for me(ptcl data
220 ! show some, but seems wrong).
221 ! introduce a cut at a small x(=eps).
222 ! divide the samling region into 2: x<eps and x> eps
223 ! at x>eps
224 ! use weight function exp(-x)/root(x)
225 ! x**s exp(-x)/rho= x**(s+1/2)
226 ! at x<eps
227 ! solve the usual equation approxmately. for small x,
228 ! int(0;x)(x**s exp(-x)/gamma(s+1))=
229 ! x**(s+1)/gamma(s+2)(1-(s+1)/(s+2)x+(s+1)/(s+3)/2x**2
230 ! then, the relative area of x< eps
231  br=eps**(s+1.)/kgamma(s+2.) * (1.- (s+1.)/(s+2.)*eps
232  * + (s+1.)/(s+3.)/2*eps*eps )
233  call rndc(u)
234  if(u .lt. br) then
235 ! take value x< eps
236  call rndc(u)
237  us=u**(1./(s+1.))*eps
238 ! initial guess at x<eps
239  x=us
240  tmp=
241  * (1.- (s+1.)/(s+2.)*eps+(s+1.)/(s+3)/2*eps*eps)
242  * **(1./(s+1.))
243  more=.true.
244  do 130 while( more )
245 ! average trial is < 2
246  xold=x
247  x=us*tmp/
248  * (1.-(s+1.)/(s+2.)*xold+(s+1.)/(s+3)/2*xold*xold)
249  * ** (1./(s+1.))
250  if(x .lt. small) then
251  acc=0.
252  else
253  acc=abs(x/xold-1.)
254  endif
255  more=acc .gt. err
256  130 continue
257  else
258  more=.true.
259  do 200 while ( more )
260  more2=.true.
261  do 150 while ( more2 )
262  call kgauss(0.d0, rt2i, x)
263  x=x*x
264  more2= x .lt. eps
265  150 continue
266 ! rejection function value
267  r=(x/eps)**(s+0.5)
268  call rndc(u)
269  more=u .gt. r
270  200 continue
271  endif
272  endif
273  end
274 
275 
276 
ksgmis
subroutine ksgmis(n, a, x)
Definition: ksgamd.f:140
ksgmrm
subroutine ksgmrm(s, av, x)
Definition: ksgamd.f:125
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
ksgamn
subroutine ksgamn(s, x)
Definition: ksgamd.f:183
rndc
subroutine rndc(u)
Definition: rnd.f:91
kgauss
subroutine kgauss(m, v, g1)
Definition: kgauss.f:10
d0
block data cblkEvhnp ! currently usable models data RegMdls ad *special data *Cekaon d0
Definition: cblkEvhnp.h:5
d
dE dx *! Nuc Int sampling table d
Definition: cblkMuInt.h:130
ksgmim
subroutine ksgmim(n, av, x)
Definition: ksgamd.f:167
ksgmrs
subroutine ksgmrs(s, a, x)
Definition: ksgamd.f:73