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

Go to the source code of this file.

Functions/Subroutines

subroutine cinippx
 
real *8 function cinteldndx (x)
 
subroutine cinippxn
 
real *8 function cinteldndxn (x)
 

Function/Subroutine Documentation

◆ cinippx()

subroutine cinippx ( )

Definition at line 2 of file cinippx.f.

References a, b, c, cerrormsg(), cinteldndx(), d, d0, dx, intendndx(), kbinchop(), ktrpzintegt(), n, ndndx, ppsx, and uconst.

Referenced by cs2lp().

2  implicit none
3 ! initialize leading particle x sampling.
4 !
5 
6 #include "Zcinippxc.h"
7 
8  external cinteldndx
9 
10  integer i, icon
11  real*8 s
12  real*8 x, x1, eps, ans
13  real*8 dndx(n) ! unnormalized dn/dx of leading particle at x= 0, 0.01,,
14 
15 ! to 1.0
16  data ( dndx(i), i= 1, n)/
17  1 0, 0.178704, 0.330057, 0.439452, 0.518544, 0.576176,
18  2 0.6189, 0.651, 0.677, 0.699, 0.719, 0.737, 0.754, 0.770,
19  3 0.7855, 0.800, 0.814, 0.845, 0.869, 0.894, 0.843, 0.804,
20  4 0.806, 0.9, 0.819, 0.878, 0.856, 0.831, 0.779, 0.811,
21  5 0.866, 0.831, 0.874, 0.804, 0.826, 0.781, 0.77, 0.758,
22  6 0.777, 0.746, 0.738, 0.734, 0.71, 0.704, 0.74, 0.666,
23  7 0.671, 0.669, 0.658, 0.653, 0.655, 0.631, 0.601, 0.583,
24  8 0.568, 0.546, 0.516, 0.482, 0.485, 0.455, 0.479, 0.476,
25  9 0.435, 0.419, 0.412, 0.348, 0.382, 0.347, 0.332, 0.337,
26  a 0.3, 0.335, 0.314, 0.296, 0.267, 0.284, 0.277, 0.25,
27  b 0.217, 0.253, 0.249, 0.198, 0.197, 0.194, 0.209, 0.189,
28  c 0.215, 0.193, 0.208, 0.206, 0.212, 0.193, 0.199, 0.214,
29  d 0.248, 0.254, 0.307, 0.371, 0.574, 1.914, 4.043/
30 
31  dx = 0.01d0
32  eps = 1.d-4
33 
34 ! integral of dn/dx from 0 to 1.
35  call ktrpzintegt(dndx, 1, n, 0.d0, dx, 1.d0, s)
36 !
37 ! make normalized table
38 !
39  do i = 1, n
40  ndndx(i) = dndx(i)/s
41  enddo
42 ! make table of intgral(0:x) of ndndx. for x = 0 to 1.0 step
43 ! 0.01
44  x = 0.
45  do i=2, n
46  x = x + dx
47  call ktrpzintegt(ndndx, 1, n, 0.d0, dx, x, intendndx(i) )
48  enddo
49  intendndx(1) = 0.
50  intendndx(n) = 1.
51 ! solve inte(0:x) of ndndx = u for u = 0 to 1.0 step 0.01
52  x = dx
53  uconst = 0.
54  x1 = 0.
55  do i = 2, n-1
56  uconst = uconst + dx
57  call kbinchop(cinteldndx, x1, 1.d0, x, eps, ans, icon)
58  if(icon .ne. 0) then
59  call cerrormsg(
60  * 'failed in making leading particle sampling table', 0)
61  endif
62  ppsx(i) = ans
63  x = ans
64  x1 = x
65  enddo
66  ppsx(1) = 0.
67  ppsx(n) = 1.
cerrormsg
subroutine cerrormsg(msg, needrtn)
Definition: cerrorMsg.f:4
x1
block data include Zlatfit h c fitting region data x1(1)/0.03/
ans
atmos%rho(atmos%nodes) **exp(-(z-atmos%z(atmos%nodes))/Hinf) elseif(z .lt. atmos%z(1)) then ans=atmos%rho(1) **exp((atmos%z(1) -z)/atmos%H(1)) else call kdwhereis(z, atmos%nodes, atmos%z, 1, i) a=atmos%a(i) if(a .ne. 0.d0) then ans=atmos%rho(i) **(1+a *(z-atmos%z(i))/atmos%H(i)) **(-1.0d0-1.d0/a) else ans=*atmos%rho(i) *exp(-(z-atmos%z(i))/atmos%H(i)) endif endif ! zsave=z ! endif cvh2den=ans end ! ---------------------------------- real *8 function cvh2temp(z) implicit none ! vettical height to temperatur(Kelvin) real *8 z ! input. vertical height in m ! output is temperature of the atmospher in Kelvin real *8 ans integer i if(z .gt. atmos%z(atmos%nodes)) then ans=atmos%T(atmos%nodes) elseif(z .lt. atmos%z(1)) then ans=atmos%T(1)+atmos%b(1) *(z - atmos%z(1)) else call kdwhereis(z, atmos%nodes, atmos%z, 1, i) ans=atmos%T(i)+atmos%b(i) *(z-atmos%z(i)) endif cvh2temp=ans end !--------------------------------------------- real *8 function cthick2h(t) implicit none real *8 t ! input. air thickness in kg/m^2 real *8 logt, ans integer i real *8 dod0, fd, a logt=log(t) if(t .ge. atmos%cumd(1)) then ans=atmos%z(1) - *(logt - atmos%logcumd(1)) *atmos%H(1) elseif(t .le. atmos%cumd(atmos%nodes)) then ans=atmos%z(atmos%nodes) - *Hinf *log(t/atmos%cumd(atmos%nodes)) else call kdwhereis(t, atmos%nodes, atmos%cumd, 1, i) ! i is such that X(i) > x >=x(i+1) ans
Definition: cstdatmos0.fNew2.h:82
i
nodes i
Definition: cstdatmos0.fNew2.h:93
dx
integer npitbl real *nx dx real dx
Definition: Zcinippxc.h:10
cinteldndx
real *8 function cinteldndx(x)
Definition: cinippx.f:70
ndndx
integer npitbl real *nx dx real ndndx
Definition: Zcinippxc.h:10
kbinchop
subroutine kbinchop(f, x1, x2, x, eps, ans, icon)
Definition: kbinChop.f:20
uconst
integer npitbl real *nx dx real uconst
Definition: Zcinippxc.h:10
intendndx
integer npitbl real *nx dx real intendndx(n) real *8 ndndxn(n)
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
ppsx
integer npitbl real *nx dx real ppsx
Definition: Zcinippxc.h:10
ktrpzintegt
subroutine ktrpzintegt(t, intv, n, x0, dx, x, ans)
Definition: ktrpzIntegT.f:3
nrl_atmos::a
real(4), save a
Definition: cNRLAtmos.f:20
nrl_atmos::b
real(4), save b
Definition: cNRLAtmos.f:21
n
integer n
Definition: Zcinippxc.h:1
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
c
dE dx *! Nuc Int sampling table c
Definition: cblkMuInt.h:130
Here is the call graph for this function:
Here is the caller graph for this function:

◆ cinippxn()

subroutine cinippxn ( )

Definition at line 81 of file cinippx.f.

References a, b, c, cerrormsg(), cinteldndxn(), d, d0, dx, intendndxn, kbinchop(), ktrpzintegt(), n, ndndxn, and uconst.

Referenced by cs2lp().

81  implicit none
82 ! initialize leading particle x sampling.
83 ! for p--> n case
84 !
85 
86 #include "Zcinippxc.h"
87 
88  external cinteldndxn
89 
90  integer i, icon
91  real*8 s
92  real*8 x, x1, eps, ans
93  real*8 dndx(n) ! unnormalized dn/dx of leading particle at x= 0, 0.01,,
94 
95 ! to 1.0
96  data ( dndx(i), i= 1, n)/
97  1 0, 0.1161, 0.214, 0.285, 0.337, 0.374,0.402, 0.423,
98  2 0.440, 0.4549, 0.467, 0.479, 0.490, 0.50, 0.510,
99  3 0.520, 0.529, 0.549, 0.564, 0.581, 0.58, 0.59, 0.59,
100  4 0.59, 0.58, 0.58, 0.57, 0.564, 0.543, 0.532, 0.521,
101  5 0.510, 0.505, 0.50, 0.50, 0.50, 0.50, 0.49, 0.48,
102  6 0.47, 0.46, 0.45, 0.45, 0.44, 0.44,0.43, 0.41, 0.397,
103  7 0.39, 0.38, 0.375, 0.37, 0.34, 0.33, 0.32, 0.315, 0.31,
104  8 0.30, 0.295, 0.29, 0.28, 0.27, 0.25, 0.24, 0.23, 0.22,
105  9 0.21, 0.20, 0.195, 0.19, 0.18, 0.172, 0.17, 0.16, 0.15,
106  a 0.145, 0.14, 0.13, 0.121, 0.115, 0.102, 0.095, 0.091,
107  b 0.0808, 0.076, 0.073, 0.070, 0.067, 0.062, 0.056, 0.054,
108  c 0.047, 0.041, 0.037, 0.034, 0.029, 0.025, 0.020, 0.015,
109  d 0.006, 0.001/
110 
111 
112  dx = 0.01d0
113  eps = 1.d-4
114 
115 ! integral of dn/dx from 0 to 1.
116  call ktrpzintegt(dndx, 1, n, 0.d0, dx, 1.d0, s)
117 !
118 ! make normalized table
119 !
120  do i = 1, n
121  ndndxn(i) = dndx(i)/s
122  enddo
123 ! make table of intgral(0:x) of ndndx. for x = 0 to 1.0 step
124 ! 0.01
125  x = 0.
126  do i=2, n
127  x = x + dx
128  call ktrpzintegt(ndndxn, 1, n, 0.d0, dx, x, intendndxn(i) )
129  enddo
130  intendndxn(1) = 0.
131  intendndxn(n) = 1.
132 ! solve inte(0:x) of ndndx = u for u = 0 to 1.0 step 0.01
133  x = dx
134  uconst = 0.
135  x1 = 0.
136  do i = 2, n-1
137  uconst = uconst + dx
138  call kbinchop(cinteldndxn, x1, 1.d0, x, eps, ans, icon)
139  if(icon .ne. 0) then
140  call cerrormsg(
141  * 'failed in making leading particle sampling table', 0)
142  endif
143  ppsxn(i) = ans
144  x = ans
145  x1 = x
146  enddo
147  ppsxn(1) = 0.
148  ppsxn(n) = 1.
cerrormsg
subroutine cerrormsg(msg, needrtn)
Definition: cerrorMsg.f:4
x1
block data include Zlatfit h c fitting region data x1(1)/0.03/
ans
atmos%rho(atmos%nodes) **exp(-(z-atmos%z(atmos%nodes))/Hinf) elseif(z .lt. atmos%z(1)) then ans=atmos%rho(1) **exp((atmos%z(1) -z)/atmos%H(1)) else call kdwhereis(z, atmos%nodes, atmos%z, 1, i) a=atmos%a(i) if(a .ne. 0.d0) then ans=atmos%rho(i) **(1+a *(z-atmos%z(i))/atmos%H(i)) **(-1.0d0-1.d0/a) else ans=*atmos%rho(i) *exp(-(z-atmos%z(i))/atmos%H(i)) endif endif ! zsave=z ! endif cvh2den=ans end ! ---------------------------------- real *8 function cvh2temp(z) implicit none ! vettical height to temperatur(Kelvin) real *8 z ! input. vertical height in m ! output is temperature of the atmospher in Kelvin real *8 ans integer i if(z .gt. atmos%z(atmos%nodes)) then ans=atmos%T(atmos%nodes) elseif(z .lt. atmos%z(1)) then ans=atmos%T(1)+atmos%b(1) *(z - atmos%z(1)) else call kdwhereis(z, atmos%nodes, atmos%z, 1, i) ans=atmos%T(i)+atmos%b(i) *(z-atmos%z(i)) endif cvh2temp=ans end !--------------------------------------------- real *8 function cthick2h(t) implicit none real *8 t ! input. air thickness in kg/m^2 real *8 logt, ans integer i real *8 dod0, fd, a logt=log(t) if(t .ge. atmos%cumd(1)) then ans=atmos%z(1) - *(logt - atmos%logcumd(1)) *atmos%H(1) elseif(t .le. atmos%cumd(atmos%nodes)) then ans=atmos%z(atmos%nodes) - *Hinf *log(t/atmos%cumd(atmos%nodes)) else call kdwhereis(t, atmos%nodes, atmos%cumd, 1, i) ! i is such that X(i) > x >=x(i+1) ans
Definition: cstdatmos0.fNew2.h:82
i
nodes i
Definition: cstdatmos0.fNew2.h:93
dx
integer npitbl real *nx dx real dx
Definition: Zcinippxc.h:10
kbinchop
subroutine kbinchop(f, x1, x2, x, eps, ans, icon)
Definition: kbinChop.f:20
uconst
integer npitbl real *nx dx real uconst
Definition: Zcinippxc.h:10
d0
block data cblkEvhnp ! currently usable models data RegMdls ad *special data *Cekaon d0
Definition: cblkEvhnp.h:5
intendndxn
integer npitbl real *nx dx real * intendndxn
Definition: Zcinippxc.h:10
d
dE dx *! Nuc Int sampling table d
Definition: cblkMuInt.h:130
ktrpzintegt
subroutine ktrpzintegt(t, intv, n, x0, dx, x, ans)
Definition: ktrpzIntegT.f:3
cinteldndxn
real *8 function cinteldndxn(x)
Definition: cinippx.f:151
nrl_atmos::a
real(4), save a
Definition: cNRLAtmos.f:20
nrl_atmos::b
real(4), save b
Definition: cNRLAtmos.f:21
n
integer n
Definition: Zcinippxc.h:1
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
ndndxn
integer npitbl real *nx dx real ndndxn
Definition: Zcinippxc.h:10
c
dE dx *! Nuc Int sampling table c
Definition: cblkMuInt.h:130
Here is the call graph for this function:
Here is the caller graph for this function:

◆ cinteldndx()

real*8 function cinteldndx ( real*8  x)

Definition at line 70 of file cinippx.f.

References d0, dx, ktrpzintegt(), n, ndndx, and uconst.

Referenced by cinippx().

70  implicit none
71  real*8 x
72 
73 #include "Zcinippxc.h"
74 
75  real*8 ans
76 
77  call ktrpzintegt(ndndx, 1, n, 0.d0, dx, x, ans)
78  cinteldndx = ans - uconst
ans
atmos%rho(atmos%nodes) **exp(-(z-atmos%z(atmos%nodes))/Hinf) elseif(z .lt. atmos%z(1)) then ans=atmos%rho(1) **exp((atmos%z(1) -z)/atmos%H(1)) else call kdwhereis(z, atmos%nodes, atmos%z, 1, i) a=atmos%a(i) if(a .ne. 0.d0) then ans=atmos%rho(i) **(1+a *(z-atmos%z(i))/atmos%H(i)) **(-1.0d0-1.d0/a) else ans=*atmos%rho(i) *exp(-(z-atmos%z(i))/atmos%H(i)) endif endif ! zsave=z ! endif cvh2den=ans end ! ---------------------------------- real *8 function cvh2temp(z) implicit none ! vettical height to temperatur(Kelvin) real *8 z ! input. vertical height in m ! output is temperature of the atmospher in Kelvin real *8 ans integer i if(z .gt. atmos%z(atmos%nodes)) then ans=atmos%T(atmos%nodes) elseif(z .lt. atmos%z(1)) then ans=atmos%T(1)+atmos%b(1) *(z - atmos%z(1)) else call kdwhereis(z, atmos%nodes, atmos%z, 1, i) ans=atmos%T(i)+atmos%b(i) *(z-atmos%z(i)) endif cvh2temp=ans end !--------------------------------------------- real *8 function cthick2h(t) implicit none real *8 t ! input. air thickness in kg/m^2 real *8 logt, ans integer i real *8 dod0, fd, a logt=log(t) if(t .ge. atmos%cumd(1)) then ans=atmos%z(1) - *(logt - atmos%logcumd(1)) *atmos%H(1) elseif(t .le. atmos%cumd(atmos%nodes)) then ans=atmos%z(atmos%nodes) - *Hinf *log(t/atmos%cumd(atmos%nodes)) else call kdwhereis(t, atmos%nodes, atmos%cumd, 1, i) ! i is such that X(i) > x >=x(i+1) ans
Definition: cstdatmos0.fNew2.h:82
dx
integer npitbl real *nx dx real dx
Definition: Zcinippxc.h:10
cinteldndx
real *8 function cinteldndx(x)
Definition: cinippx.f:70
ndndx
integer npitbl real *nx dx real ndndx
Definition: Zcinippxc.h:10
uconst
integer npitbl real *nx dx real uconst
Definition: Zcinippxc.h:10
d0
block data cblkEvhnp ! currently usable models data RegMdls ad *special data *Cekaon d0
Definition: cblkEvhnp.h:5
ktrpzintegt
subroutine ktrpzintegt(t, intv, n, x0, dx, x, ans)
Definition: ktrpzIntegT.f:3
n
integer n
Definition: Zcinippxc.h:1
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:

◆ cinteldndxn()

real*8 function cinteldndxn ( real*8  x)

Definition at line 151 of file cinippx.f.

References d0, dx, ktrpzintegt(), n, ndndxn, and uconst.

Referenced by cinippxn().

151  implicit none
152  real*8 x
153 
154 #include "Zcinippxc.h"
155 
156  real*8 ans
157 
158  call ktrpzintegt(ndndxn, 1, n, 0.d0, dx, x, ans)
159  cinteldndxn = ans - uconst
ans
atmos%rho(atmos%nodes) **exp(-(z-atmos%z(atmos%nodes))/Hinf) elseif(z .lt. atmos%z(1)) then ans=atmos%rho(1) **exp((atmos%z(1) -z)/atmos%H(1)) else call kdwhereis(z, atmos%nodes, atmos%z, 1, i) a=atmos%a(i) if(a .ne. 0.d0) then ans=atmos%rho(i) **(1+a *(z-atmos%z(i))/atmos%H(i)) **(-1.0d0-1.d0/a) else ans=*atmos%rho(i) *exp(-(z-atmos%z(i))/atmos%H(i)) endif endif ! zsave=z ! endif cvh2den=ans end ! ---------------------------------- real *8 function cvh2temp(z) implicit none ! vettical height to temperatur(Kelvin) real *8 z ! input. vertical height in m ! output is temperature of the atmospher in Kelvin real *8 ans integer i if(z .gt. atmos%z(atmos%nodes)) then ans=atmos%T(atmos%nodes) elseif(z .lt. atmos%z(1)) then ans=atmos%T(1)+atmos%b(1) *(z - atmos%z(1)) else call kdwhereis(z, atmos%nodes, atmos%z, 1, i) ans=atmos%T(i)+atmos%b(i) *(z-atmos%z(i)) endif cvh2temp=ans end !--------------------------------------------- real *8 function cthick2h(t) implicit none real *8 t ! input. air thickness in kg/m^2 real *8 logt, ans integer i real *8 dod0, fd, a logt=log(t) if(t .ge. atmos%cumd(1)) then ans=atmos%z(1) - *(logt - atmos%logcumd(1)) *atmos%H(1) elseif(t .le. atmos%cumd(atmos%nodes)) then ans=atmos%z(atmos%nodes) - *Hinf *log(t/atmos%cumd(atmos%nodes)) else call kdwhereis(t, atmos%nodes, atmos%cumd, 1, i) ! i is such that X(i) > x >=x(i+1) ans
Definition: cstdatmos0.fNew2.h:82
dx
integer npitbl real *nx dx real dx
Definition: Zcinippxc.h:10
uconst
integer npitbl real *nx dx real uconst
Definition: Zcinippxc.h:10
d0
block data cblkEvhnp ! currently usable models data RegMdls ad *special data *Cekaon d0
Definition: cblkEvhnp.h:5
ktrpzintegt
subroutine ktrpzintegt(t, intv, n, x0, dx, x, ans)
Definition: ktrpzIntegT.f:3
cinteldndxn
real *8 function cinteldndxn(x)
Definition: cinippx.f:151
n
integer n
Definition: Zcinippxc.h:1
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
ndndxn
integer npitbl real *nx dx real ndndxn
Definition: Zcinippxc.h:10
Here is the call graph for this function:
Here is the caller graph for this function: