COSMOS v7.655  COSMOSv7655
(AirShowerMC)
cstdatmos0.fOld.h
Go to the documentation of this file.
1 ! Units are in SI.
2 
3 ! below cthick2h1 and cthick2h2 are so made that they coincide at
4 ! t=230
5 
6 !
7 ! -------------------------------------
8  real*8 function cvh2den(z)
9 ! --------------------------vertical height to density
10 ! at z > 50 km, very bad. but ok for our purpose
11  implicit none
12 !---- include 'Zstdatmos.h'
13 #include "Zstdatmos.h"
14  real*8 z
15  real*8 zsave, ans
16  data zsave/-1.d30/
17  save zsave, ans
18 
19 
20  if(z .ne. zsave) then
21  if(z .ge. hc) then
22  ans = 10.d0*exp( (hb-z)/hn)/hn
23  else
24  ans = 10.d0*hlhmi*( (ha-z)/hl)**(hmi-1.d0)
25  endif
26  zsave = z
27  endif
28  cvh2den = ans
29  end
30 !
31 ! -------------------------------------------------------------
32 ! real*8 z, cvh2temp
33 ! do z = 0., 1250.d3, 1d3
34 ! write(*, *) sngl(z/1.d3), sngl(cvh2temp(z))
35 ! enddo
36 ! end
37 
38  real*8 function cvh2temp(z)
39 ! vettical height to temperatur (Kelvin)
40 
41  real*8 z ! input. vertical height in m
42 ! output is temperature of the atmospher in Kelvin
43 
44  real*8 ans, error, hx
45  integer j, m, loc, nh
46 
47  parameter(nh=17)
48 
49  real*8 h(nh), t(nh) ! h (km) vs t
50 
51  data h/91, 100, 110, 120, 130, 140, 160, 180, 200, 250,
52  * 300, 350, 400, 450, 500, 550, 600/
53  data t/186.87, 195.08, 240.0, 360.0, 469.27, 559.63,
54  * 696.29, 790.07, 854.56, 941.33, 976.01, 990.06, 995.83,
55  * 998.22, 999.24, 999.67, 999.85/
56  data m/2/
57 
58  if(z .lt. 11.1d3) then
59  cvh2temp = (216.65d0 - 288.15d0)/11.1d3 * z + 288.15d0
60  elseif(z .lt. 20.0d3) then
61  cvh2temp = 216.65d0
62  elseif(z .lt. 32.2d3) then
63  cvh2temp = (228.756d0 - 216.65d0)/(32.2d3-20.d3) * (z-20.d3)
64  * + 216.65d0
65  elseif(z .lt. 47.4d3) then
66  cvh2temp = (270.65d0-228.756d0)/(47.4d3-32.2d3) * (z-32.2d3)
67  * + 228.756d0
68  elseif(z .lt. 51.d3) then
69  cvh2temp = 270.650
70  elseif(z .lt. 72.d3) then
71  cvh2temp = (214.263d0 - 270.65d0)/(72.d3- 51.0d3) * (z-51.d3)
72  * + 270.65d0
73  elseif(z .lt. 86.d3) then
74  cvh2temp = (186.87d0-214.263d0)/(86.d3-72.d3)* (z -72.d3)
75  * + 214.263d0
76  elseif( z .lt. 91.d3) then
77  cvh2temp = 186.87d0
78  elseif( z .lt. 600d3) then
79  hx = z/1.d3 ! in km
80  call kdwhereis(hx, nh, h, 1, loc)
81  j = min(max(loc - (m-1)/2,1), nh+1-m) ! max of m points from j
82  call kpolintp(h(j), 1, t(j), 1, m, hx, ans, error)
83  cvh2temp = ans
84  else
85  cvh2temp = 1000.
86  endif
87  end
88 
89 !---------------------------------------------
90  real*8 function cthick2h(t)
91  implicit none
92  real*8 t
93 #include "Zstdatmos.h"
94  external cblkStdAtmos
95 
96 
97  if(t .lt. t0) then
98  cthick2h = hb - log(t/10.d0) * hn ! for t<t0 = 2300
99  else
100  cthick2h = ha - hl * (t/10.d0)**hm ! for t > t0
101  endif
102  end
103 ! -------------------------------------
104  real*8 function cthick2den(t)
105 ! -------------------------------------
106  implicit none
107 !---- include 'Zstdatmos.h'
108 #include "Zstdatmos.h"
109  real*8 t
110 
111  if(t .lt. t0) then
112  cthick2den = t/hn
113  else
114  cthick2den =10.d0*hlhmi*
115  * (t/10.d0)**(1.d0 - hm)
116  endif
117  end
118 ! -------------------------------------
119  real*8 function cvh2denp(z)
120 ! -------------------------------------
121 
122 ! d rho/dz
123  implicit none
124 !---- include 'Zstdatmos.h'
125 #include "Zstdatmos.h"
126  real*8 z
127  real*8 zsave, ans, cvh2den
128  data zsave /-1.d30/
129  save ans, zsave
130 
131  if(z .ne. zsave) then
132  if(z .ge. hc) then
133 ! cvh2denp = - 10.d0*exp( (hb-z)/hn)/hn/hn
134 ! cvh2denp = cnst1*exp( (hb-z)/hn)
135  ans = - cvh2den(z)/hn
136  else
137 ! cvh2denp =- 10.d0*hlhmi*(hmi-1.d0)/hl*
138 ! * ( (ha-z)/hl)**(hmi-2.d0)
139 ! cvh2denp =cnst2 *
140 ! * ( (ha-z)/hl)**cnst3
141  ans = -(hmi-1.0d0)/(ha-z)*cvh2den(z)
142  endif
143  zsave= z
144  endif
145  cvh2denp = ans
146  end
147 ! ----------------------------------
148  real*8 function cvh2scaleh(z)
149 ! ----------------------------------
150 ! get scale height defined by
151 ! - rho/ (d rho/d z). This has discontinuity at
152 ! z = 11.1 km
153  implicit none
154  real*8 z
155  real*8 cvh2den, cvh2denp
156 !
157  cvh2scaleh = - cvh2den(z)/ cvh2denp(z)
158  end
159 ! -------------------------------------
160  real*8 function cvh2den2p(z)
161 ! -------------------------------------
162 ! d(d rho/dz)/dz
163  implicit none
164 !---- include 'Zstdatmos.h'
165 #include "Zstdatmos.h"
166  real*8 z
167 ! logical firsttime/.true./
168 ! save firsttime
169 ! real*8 cnst1, cnst2, cnst3
170 ! save cnst1, cnst2, cnst3
171  real*8 zsave, ans, cvh2denp, cvh2den
172  data zsave/-1.d30/
173  save zsave
174 !
175 ! if(firsttime) then
176 ! cnst1 = 10.d0/hn/hn/hn
177 ! cnst2 = 10.d0*hlhmi*(hmi-1.d0)*(hmi-2.d0)
178 ! * /hl/hl
179 ! cnst3 = hmi -3.d0
180 ! firsttime = .false.
181 ! endif
182  if(z .ne. zsave) then
183  if(z .ge. hc) then
184 ! cvh2den2p = 10.d0*exp( (hb-z)/hn)/hn/hn/hn
185 ! cvh2den2p = cnst1*exp( (hb-z)/hn)
186  ans = cvh2den(z)/hn/hn
187  else
188 ! cvh2den2p = 10.d0*hlhmi*(hmi-1.d0)*(hmi-2.d0)
189 ! * /hl/hl *
190 ! * ( (ha-z)/hl)**(hmi-3.d0)
191 ! cvh2den2p = cnst2 * ( (ha-z)/hl )** cnst3
192  ans = -(hmi-2.d0)/(ha -z) * cvh2denp(z)
193  endif
194  zsave =z
195  endif
196  cvh2den2p = ans
197  end
198 ! -------------------------------------
199  real*8 function cvh2den3p(z)
200 ! -------------------------------------
201 ! rho'''(z)
202  implicit none
203 !---- include 'Zstdatmos.h'
204 #include "Zstdatmos.h"
205  real*8 z
206 ! logical firsttime/.true./
207 ! save firsttime
208 ! real*8 cnst1, cnst2, cnst3
209 ! save cnst1, cnst2, cnst3
210  real*8 zsave, ans, cvh2den2p, cvh2den
211  data zsave/-1.d30/
212  save zsave, ans
213 !
214 ! if(firsttime) then
215 ! cnst1 = - 10.d0/hn/hn/hn/hn
216 ! cnst2 =- 10.d0*hlhmi*(hmi-1.d0)*(hmi-2.d0)*(hmi-3.d0)
217 ! * /hl/hl/hl
218 ! cnst3 = hmi -4.d0
219 ! firsttime = .false.
220 ! endif
221  if(z .ne. zsave) then
222  if(z .ge. hc) then
223 ! cvh2den3p = - 10.d0*exp( (hb-z)/hn)/hn/hn/hn/hn
224 ! cvh2den3p = cnst1*exp( (hb-z)/hn)
225  ans = - cvh2den(z)/hn/hn/hn
226  else
227 ! cvh2den3p =- 10.d0*hlhmi*(hmi-1.d0)*(hmi-2.d0)*(hmi-3)
228 ! * /hl/hl/hl *
229 ! * ( (ha-z)/h1)**(hmi-4.d0)
230 ! cvh2den3p = cnst2 * ( (ha-z)/hl )** cnst3
231  ans = -(hmi-3.0d0)/(ha-z) * cvh2den2p(z)
232  endif
233  zsave =z
234  endif
235  cvh2den3p = ans
236  end
237 ! ---------------------------------------
238  real*8 function cvh2thick(z)
239 ! ---------------------------------------
240  implicit none
241 !---- include 'Zstdatmos.h'
242 #include "Zstdatmos.h"
243  real*8 z
244 
245  if(z .gt. hc) then
246  cvh2thick = 10.d0* exp( (hb-z)/hn)
247  else
248  cvh2thick = 10.d0* ((ha-z)/hl)**hmi
249 !
250 ! the next factor is -1.5d-11 and is neglected.
251 ! * - 10.d0* ((ha-hc)/hl)**hmi +
252 ! * 10.d0*exp( (hb-hc)/hn)
253 !
254  endif
255  end
256  block data cblkStdAtmos
257 #include "Zstdatmos.h"
258 ! main frame cosmos data.
259 ! d rho/dz is discontinuous at hc, but others are
260 ! continuous. Thickness values are almost the same
261 ! as the standard obained by cstdatmos0-multi-seg.f
262 !
263  data ha/4512224.7572830657d-2/, hb/4541933.9782793734d-2/,
264  * hl/1244541.6061892177d-2/, hm/.1854893358365053d0/,
265  * hn/6.3300000000746224d3/, t0/230.000458235099d1/,
266  * hc /10.996296495761545d3/
267  data hmi/5.3911455097420/, hlhmi/4.3318322850207d-04/
268 
269 ! coef. to have continuous density values for h=h(t) and its derivative
270 ! about t=t0. However, the absolute thickness is little bit larger
271 ! than the standard, so that we don't use this one.
272 !
273 ! data ha/25289.093750000d0/, hb/48490.145953733/,
274 ! * hl/2712.61d0/, hm/0.32397778012688d0/,
275 ! * hn/6800 .0/, hmi/3.0866314338235/, hc/11.1d3/,
276 ! * hlhmi/1.1378815278872D-03/, t0/2443.3759999999/
277  end
278 
279 
parameter
integer npitbl real *nx parameter(n=101, npitbl=46, nx=n-1) real *8 uconst
modpdgxs::m
real(8), parameter, public m
Definition: cpdgXs.f:13
z
nodes z
Definition: cstdatmos0.fNew2.h:93
real
float real
Definition: Zdef.h:2
integer
int integer
Definition: Zdef.h:3
cvh2denp
real *8 function cvh2denp(z)
Definition: ciniSegAtoms.f:201
hm
! common for std atmosphere real *AlmostVacH ! almost vacuum height real *AlmostVacT ! thickness there hm
Definition: Zstdatmos.h:7
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
cvh2den3p
real *8 function cvh2den3p(z)
Definition: vh2den.f:159
hb
! common for std atmosphere real *AlmostVacH ! almost vacuum height real *AlmostVacT ! thickness there hb
Definition: Zstdatmos.h:7
kpolintp
subroutine kpolintp(xa, xstep, ya, ystep, n, x, y, error)
Definition: kpolintp.f:163
d3
block data cblkIncident data *Za1ry *HeightOfInj d3
Definition: cblkIncident.h:7
kdwhereis
subroutine kdwhereis(x, in, a, step, loc)
Definition: kdwhereis.f:27
hc
! common for std atmosphere real *AlmostVacH ! almost vacuum height real *AlmostVacT ! thickness there hc
Definition: Zstdatmos.h:7
data
block data cblkMuInt integer i ! Pair total X sec data(MuPrTX(i), i=1, 38)/1 0.431132E-02
nrl_atmos::h
real(4), dimension(:), allocatable, save h
Definition: cNRLAtmos.f:28
ne
int ne[nl][nth]
Definition: Zprivate.h:11
as
! Parameters used for hadronic cascade shower is generated newline ! For you may give as as or em as(order/case/separator insensitive) is to generate EM-cascade and AS. \newline ! Generate
j
********************block data cblkHeavy ********************integer j data *HeavyG2symbol *data *HeavyG2code kiron data j
Definition: cblkHeavy.h:36
cthick2den
real(8) function cthick2den(t)
Definition: cNRLAtmos.f:591
to
integer maxbin nregion c minsave drx2 ! drawing region real maxdep integer maxpos integer to
Definition: Zfit.h:15
block
integer maxFiles ! actual logical file ! it will be used during mpi skeleton making time ! so when the program starts to normal ! it will not be used block
Definition: Zmpi.h:23
hl
! common for std atmosphere real *AlmostVacH ! almost vacuum height real *AlmostVacT ! thickness there hl
Definition: Zstdatmos.h:7
d0
block data cblkEvhnp ! currently usable models data RegMdls ad *special data *Cekaon d0
Definition: cblkEvhnp.h:5
However
! Parameters used for hadronic cascade shower is generated newline ! For you may give as as or em quick generation of AS for heavy primaries is tried See chookASbyH f character *Generate2 don t touch this for skeleton flesh use integer MagBrem no magnetic bremsstrahlung is considered newline ! if and Ee energy loss due to magnetic brems is considered newline ! if and Ee real sampling of gamma is performed WaitRatio ! must be made small so that WaitRatio *E0 sim MagBremEmin integer MagPair no magnetic pair creation is considered newline ! if and Eg real sampling is the LPM effect is considered when Ee LpmBremEmin for electrons and ! Eg LpmPairEmin for gamma rays real *MagBremEmin E magnetic bremsstrahlung by electrons may be considered However
Definition: Ztrackp.h:169
defined
! this is to treat the stange code generated in Lund Fritiof ! If you find that you can get reliable events by discarding ! events with such strange then you may define ! the machine definition here ! If you want to discard the strange code events but ! want to be define below ! if you want discard the strage code events without ! any define below If DEBUG_STRANGECODE is you should not give the below ! if none of above two is defined
Definition: ZstrangeCode.h:3
height
*Zfirst p fm *Zfirst p Zfirst p Zfirst p *Zfirst p *Zfirst pos xyz Zfirst pos xyz *Zfirst pos *Zfirst pos Zfirst pos height
Definition: ZavoidUnionMap.h:1
cvh2den2p
real *8 function cvh2den2p(z)
Definition: ciniSegAtoms.f:251
logical
int logical
Definition: Zdef.h:1
cthick2h
nodes a atmos atmos temp real * cthick2h
Definition: cstdatmos0.fNew2.h:107
hmi
! common for std atmosphere real *AlmostVacH ! almost vacuum height real *AlmostVacT ! thickness there hmi
Definition: Zstdatmos.h:7
d
dE dx *! Nuc Int sampling table d
Definition: cblkMuInt.h:130
cvh2scaleh
real(8) function cvh2scaleh(vh)
Definition: cNRLAtmos.f:635
if
integer lengid integer lengdir character *dir integer kgetenv2 character *numb character *execid character *msg logical takehist save do nsites tkarspec nsites if(histdep(i) .gt. 0) then takehist
ha
! common for std atmosphere real *AlmostVacH ! almost vacuum height real *AlmostVacT ! thickness there hlhmi common comstdatm ha
Definition: Zstdatmos.h:7
cvh2den
real *8 function cvh2den(z)
Definition: ciniSegAtoms.f:54
nh
int nh[nl][nth]
Definition: Zprivate.h:13
main
program main
Definition: ascii2bin.f:1
t0
! common for std atmosphere real *AlmostVacH ! almost vacuum height real *AlmostVacT ! thickness there t0
Definition: Zstdatmos.h:7
cstdatmos0
subroutine cstdatmos0
Definition: cstdatmos0-multi-seg.f:38
t
nodes t
Definition: cstdatmos0.fNew2.h:84
hn
! common for std atmosphere real *AlmostVacH ! almost vacuum height real *AlmostVacT ! thickness there hn
Definition: Zstdatmos.h:7
SI
! Zkcele h ! unit here is SI(m, s, kg) real *8 sidcor
from
integer maxbin nregion c minsave drx2 ! drawing region real maxdep integer maxpos integer from
Definition: Zfit.h:15
cvh2thick
real *8 function cvh2thick(z)
Definition: ciniSegAtoms.f:93
cvh2temp
real(8) function cvh2temp(vh)
Definition: cNRLAtmos.f:61
is
block data cblkIncident data *Za1ry is
Definition: cblkIncident.h:5