COSMOS v7.655  COSMOSv7655
(AirShowerMC)
ksampPw.f File Reference

Go to the source code of this file.

Functions/Subroutines

subroutine ksamppw (ini, coef, power, node, n, x, fx)
 
subroutine ksamppw2 (coef, power, xp, cum, intgv, np, x, fx)
 
subroutine ksamppwx (c, p, x0, n, xp, c2)
 
subroutine ksamppwmrg (c1, p1, x1, n1, c2, p2, x2, n2, c, p, x, n)
 

Function/Subroutine Documentation

◆ ksamppw()

subroutine ksamppw ( integer  ini,
real*8, dimension(n coef,
real*8, dimension(n power,
real*8, dimension(n+1)  node,
integer  n,
real*8  x,
real*8  fx 
)

Definition at line 95 of file ksampPw.f.

References d0, parameter(), and rndc().

Referenced by cmbreme1(), cmbreme2(), and main().

95  implicit none
96 !
97  integer ini ! input. if coef, power, node, n are
98  ! different from the previous call, give 1.
99  ! If all of them are the same as the
100  ! previous call, give non 1
101  integer n ! input. number of segments
102  real*8 coef(n) ! input. c of cX**(-p)dX in each segment.
103  real*8 power(n) ! input. p of cX**(-p)dX in each segment.
104  real*8 node(n+1) ! input. node(1) and node(2) are the left and
105  ! right boundaries of segment 1.,
106  ! node(2) and node(3) are for segment 2,etc
107  real*8 x ! output. A sampled x. If the input condition
108  ! is ng, x = 0 will result.
109  real*8 fx ! output. function value at x, i.e., cx**-p
110 !
111  integer nmx ! max number of nodes
112  parameter( nmx = 20 )
113 
114  real*8 inte(nmx) ! integral value of each segment region
115  real*8 ci(nmx) ! normalized cummulative integral values of inte.
116  integer nsave ! for saving n
117 
118  save inte, ci, nsave
119  real*8 u, sum, temp
120  integer i, j
121 
122  real*8 cum(*), intgv(*)
123  integer np
124 
125 !
126  if(ini .eq. 1) then
127  if(nmx .lt. n) then
128  x = 0.
129  else
130 ! get integral value of each region and sum of them
131  sum = 0.
132 ! c/(1-p) x**(1-p) from x1 to x2
133  do i = 1, n
134  if(power(i) .ne. 1.d0) then
135  inte(i) = coef(i)/(1.0-power(i))*
136  * (node(i+1)**(1.-power(i)) - node(i)**(1.-power(i)))
137  else
138  inte(i) = coef(i) * log(node(i+1)/node(i))
139  endif
140  sum = sum + inte(i)
141  enddo
142 ! normalize and make cummulative data
143  do i = 1, n
144  ci(i) = inte(i)/sum
145  enddo
146  do i = 2, n
147  ci(i) = ci(i) + ci(i-1)
148  enddo
149  ci(n) = 1.0 ! for safety
150  endif
151  nsave = n
152  endif
153 ! --------------------------------------------------
154 ! assume inte and node are ready
155 ! choose which segment, j
156  call rndc(u)
157  do i = 1, n
158  if(u .le. ci(i)) then
159  j = i
160  goto 100
161  endif
162  enddo
163 ! never come here
164 !
165  100 continue
166 
167 
168  call rndc(u)
169  if(power(j) .ne. 1.d0) then
170  temp = node(j+1)**(1.-power(j)) -
171  * inte(j)*u/( coef(j)/(1. - power(j)) ) ! = x **(1-power(j))
172 
173  x = temp**(1./(1.-power(j)))
174  else
175  x = (node(j+1)/node(j))**u * node(j)
176  endif
177  fx = coef(j) * x **(-power(j))
178  1000 continue
179  return
180 !
181 !
182 ! ********************************
183  entry ksamppwq(cum, intgv, np)
184 ! *******************************
185 ! inquire the current inner variables for the
186 ! present spectrum.
187 !
188 
189  np = nsave
190 
191  do i = 1, np
192  cum(i) = ci(i)
193  intgv(i) = inte(i)
194  enddo
integer npitbl real *nx parameter(n=101, npitbl=46, nx=n-1) real *8 uconst
nodes i
subroutine rndc(u)
Definition: rnd.f:91
********************block data cblkHeavy ********************integer j data *HeavyG2symbol *data *HeavyG2code kiron data j
Definition: cblkHeavy.h:36
real(4), dimension(:), allocatable, save temp
Definition: cNRLAtmos.f:29
block data cblkEvhnp ! currently usable models data RegMdls ad *special data *Cekaon d0
Definition: cblkEvhnp.h:5
integer n
Definition: Zcinippxc.h:1
! 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:

◆ ksamppw2()

subroutine ksamppw2 ( real*8, dimension(np)  coef,
real*8, dimension(np)  power,
real*8, dimension(np+1)  xp,
real*8, dimension(np)  cum,
real*8, dimension(np)  intgv,
integer  np,
real*8  x,
real*8  fx 
)

Definition at line 199 of file ksampPw.f.

References rndc().

199 !
200 ! This interface is used to skip the initial calculations
201 ! by using output from ksampPwQ
202 !
203  implicit none
204 !
205 
206  integer np ! input. number of segments
207  real*8 coef(np) ! input. c of X**(-p)dX in each segment.
208  real*8 power(np) ! input. p of X**(-p)dX in each segment.
209  real*8 xp(np+1) ! input. output from ksampPwQ
210  real*8 cum(np) ! input. //
211  real*8 intgv(np) ! //
212 
213  real*8 x ! output. sampled x.
214  real*8 fx ! output. funcion value at x.
215 
216  real*8 u, temp
217  integer i, j
218 ! choose which segment, j
219  call rndc(u)
220  do i = 1, np
221  if(u .le. cum(i)) then
222  j = i
223  goto 100
224  endif
225  enddo
226 !
227  100 continue
228 
229  call rndc(u)
230  temp = xp(j+1)**(1.-power(j)) -
231  * intgv(j)*u/( coef(j)/(1. - power(j)) ) ! = x **(1-power(j))
232 
233  x = temp**(1./(1.-power(j)))
234  fx = coef(j) * x**(-power(j))
nodes i
subroutine rndc(u)
Definition: rnd.f:91
********************block data cblkHeavy ********************integer j data *HeavyG2symbol *data *HeavyG2code kiron data j
Definition: cblkHeavy.h:36
real(4), dimension(:), allocatable, save temp
Definition: cNRLAtmos.f:29
! 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:

◆ ksamppwmrg()

subroutine ksamppwmrg ( real*8, dimension(n1)  c1,
real*8, dimension(n1)  p1,
real*8, dimension(n1+1)  x1,
integer  n1,
real*8, dimension(n2)  c2,
real*8, dimension(n2)  p2,
real*8, dimension(n2+1)  x2,
integer  n2,
real*8, dimension(*)  c,
real*8, dimension(*)  p,
real*8, dimension(*)  x,
integer  n 
)

Definition at line 273 of file ksampPw.f.

Referenced by cmbreme1(), and cmbreme2().

273 !
274  implicit none
275  integer n1 ! input. number of segments of 1 st function
276  ! c1* x**(-p1) dx
277  real*8 c1(n1) ! input.
278  real*8 p1(n1) ! input.
279  real*8 x1(n1+1) ! input ! nodal point of segments of 1 st func.
280  integer n2 ! input. ! the same as the 2nd func.
281  real*8 c2(n2) ! input.
282  real*8 p2(n2) ! input.
283  real*8 x2(n2+1) ! input
284  integer n ! output ! those for f1 x f2
285  real*8 c(*) ! output.
286  real*8 p(*) ! output
287  real*8 x(*) ! output (dim >= nx+1 )
288 
289 !
290 
291  integer i, j
292 !
293  i = 1
294  j= 1
295  n = 0
296 
297  x(1) = max(x1(1), x2(1))
298  do while (i .le. n1 .and. j .le. n2)
299  if(x1(i+1) .lt. x2(j+1) ) then
300  n = n + 1
301  x(n+1) = x1(i+1)
302  c(n) = c1(i) * c2(j)
303  p(n) = p1(i) + p2(j)
304  i = i + 1
305 
306  elseif(x1(i+1) .eq. x2(j+1)) then
307  n = n + 1
308  x(n+1) = x1(i+1)
309  c(n) = c1(i) * c2(j)
310  p(n) = p1(i) + p2(j)
311  i = i + 1
312  j = j + 1
313  else
314  n = n + 1
315  x(n+1) = x2(j+1)
316  c(n) = c1(i) * c2(j)
317  p(n) = p1(i) + p2(j)
318  j = j + 1
319  endif
320  enddo
block data include Zlatfit h c fitting region data x1(1)/0.03/
nodes i
********************block data cblkHeavy ********************integer j data *HeavyG2symbol p
Definition: cblkHeavy.h:7
********************block data cblkHeavy ********************integer j data *HeavyG2symbol *data *HeavyG2code kiron data j
Definition: cblkHeavy.h:36
block data include Zlatfit h c fitting region data x2(1)/0.5/data x1(2)/0.3/
integer n
Definition: Zcinippxc.h:1
! 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
dE dx *! Nuc Int sampling table c
Definition: cblkMuInt.h:130
Here is the caller graph for this function:

◆ ksamppwx()

subroutine ksamppwx ( real*8, dimension(n c,
real*8, dimension(n p,
real*8, dimension(n x0,
integer  n,
real*8, dimension(n-1)  xp,
real*8, dimension(n c2 
)

Definition at line 240 of file ksampPw.f.

Referenced by cmbreme1(), cmbreme2(), and main().

240  implicit none
241 !
242  integer n ! input. number of segments
243  real*8 c(n) ! input. c of c(X/x0)**(-p) in each segment.
244  real*8 p(n) ! input. see above
245  real*8 x0(n) ! input. see above
246  real*8 xp(n-1) ! output. crossing point of adjacent segments
247  real*8 c2(n) ! output. c2 = c/x0**p. That is, c2*x**(-p) becomes
248  ! the spectrum
249 !
250  integer i
251 
252 ! get coeff.
253  do i = 1, n
254  c2(i) = c(i)* x0(i)**p(i)
255  enddo
256 ! get crossing point of i and i+1th segmment
257 ! c_i X**(-p_i) = c_i+1 X**(-p_i+1)
258 
259  do i = 1, n-1
260  xp(i) =
261  * (c2(i)/c2(i+1) )**(1./(p(i) - p(i+1)))
262  enddo
263 
nodes i
********************block data cblkHeavy ********************integer j data *HeavyG2symbol p
Definition: cblkHeavy.h:7
integer n
Definition: Zcinippxc.h:1
dE dx *! Nuc Int sampling table c
Definition: cblkMuInt.h:130
Here is the caller graph for this function: