COSMOS v7.655  COSMOSv7655
(AirShowerMC)
cstdatmos0.fNew2.h
Go to the documentation of this file.
1 ! This is to express atmosphere using linear interpolation.
2 ! cspline is not used because it generates some odd behaviour
3 !
4 !
5 ! -------------------------------------
6  real*8 function cvh2den(z)
7 ! --------------------------vertical height to density
8  implicit none
9 #include "Zatmos.h"
10  real*8 z ! input. m
11 ! real*8 zsave
12  real*8 ans
13 ! data zsave/-1.d30/
14 ! save zsave, ans
15  integer i
16  real*8 a
17 
18 ! if(z .ne. zsave) then
19  if( z .gt. atmos%z(atmos%nodes) ) then
20  ans = atmos%rho(atmos%nodes)*
21  * exp(-(z-atmos%z(atmos%nodes))/Hinf)
22  elseif(z .lt. atmos%z(1)) then
23  ans = atmos%rho(1)*
24  * exp( (atmos%z(1)-z)/atmos%H(1) )
25  else
26  call kdwhereis(z, atmos%nodes, atmos%z, 1, i)
27  a = atmos%a(i)
28  if(a .ne. 0.d0) then
29  ans = atmos%rho(i)*
30  * (1+ a*(z-atmos%z(i))/atmos%H(i))**(-1.0d0-1.d0/a)
31  else
32  ans =
33  * atmos%rho(i) * exp(- (z-atmos%z(i))/atmos%H(i))
34  endif
35  endif
36 ! zsave = z
37 ! endif
38  cvh2den = ans
39  end
40 ! ----------------------------------
41  real*8 function cvh2temp(z)
42  implicit none
43 #include "Zatmos.h"
44 ! vettical height to temperatur (Kelvin)
45 
46  real*8 z ! input. vertical height in m
47 ! output is temperature of the atmospher in Kelvin
48 
49  real*8 ans
50  integer i
51 
52  if( z .gt. atmos%z(atmos%nodes) ) then
53  ans = atmos%T(atmos%nodes)
54  elseif(z .lt. atmos%z(1)) then
55  ans = atmos%T(1) + atmos%b(1)*(z - atmos%z(1))
56  else
57  call kdwhereis(z, atmos%nodes, atmos%z, 1, i)
58  ans = atmos%T(i) + atmos%b(i)*(z-atmos%z(i))
59  endif
60  cvh2temp = ans
61  end
62 
63 !---------------------------------------------
64  real*8 function cthick2h(t)
65  implicit none
66 #include "Zatmos.h"
67  real*8 t ! input. air thickness in kg/m^2
68 
69  real*8 logt, ans
70  integer i
71  real*8 dod0, fd, a
72 
73  logt = log(t)
74  if(t .ge. atmos%cumd(1) ) then
75  ans = atmos%z(1) -
76  * (logt - atmos%logcumd(1) )*atmos%H(1)
77  elseif(t .le. atmos%cumd(atmos%nodes)) then
78  ans = atmos%z(atmos%nodes) -
79  * Hinf* log(t/atmos%cumd(atmos%nodes))
80  else
81  call kdwhereis(t, atmos%nodes, atmos%cumd, 1, i)
82 ! i is such that X(i) > x >= x(i+1);
83 ! i=1, 2, .. nodes-1
84  dod0 =( atmos%cumd(i) - t )/ atmos%d0(i)
85  a = atmos%a(i)
86  fd = 1. - dod0
87  if(a .ne. 0.) then
88  ans = (fd**(-a)- 1.0d0)*atmos%H(i)/a + atmos%z(i)
89  else
90  ans = -log(fd)* atmos%H(i) + atmos%z(i)
91  endif
92 !
93 ! write(0,*) ' t=',t, ' z=', ans, ' i=',i, ' dod0=', dod0, ' a=',a
94 ! write(0,*) ' atmos.H,z=', atmos.H(i), atmos.z(i)
95 !
96  endif
97 
98  cthick2h = ans
99  end
100 
101 !---------------------------------------------
102  real*8 function cthick2den(t)
103  implicit none
104 #include "Zatmos.h"
105  real*8 t ! input. air thickness in kg/m^2
106 
107  real*8 logt, ans, temp
108  real*8 cthick2h, cvh2den
109 
110 
111  logt = log(t)
112  if(t .gt. atmos%cumd(1) ) then
113  temp = cthick2h(t)
114  ans = cvh2den( temp )
115  elseif(t .lt. atmos%cumd(atmos%nodes)) then
116  temp = cthick2h(t)
117  ans = cvh2den( temp )
118  else
119  temp = cthick2h(t)
120  ans = cvh2den( temp )
121  endif
122  cthick2den = ans
123  end
124 ! -------------------------------------
125  real*8 function cvh2denp(z)
126 ! -------------------------------------
127 ! d rho/dz
128  implicit none
129 #include "Zatmos.h"
130  real*8 z
131 
132  real*8 ans
133  integer i
134  real*8 a
135 
136  if( z .gt. atmos%z(atmos%nodes) ) then
137  ans = 0.
138  elseif(z .lt. atmos%z(1)) then
139  ans =- atmos%rho(1)/atmos%H(1)*
140  * exp( (atmos%z(1)-z)/atmos%H(1) )
141  else
142  call kdwhereis(z, atmos%nodes, atmos%z, 1, i)
143  a = atmos%a(i)
144  if(a .ne. 0.d0) then
145  ans = atmos%rho(i)*a/atmos%H(i) *
146  * (1.0d0 + a*(z-atmos%z(i))/atmos%H(i) )**(-2.d0-1.0d0/a)
147  else
148  ans =- atmos%rho(i)/atmos%H(i)*
149  * exp(- (z-atmos%z(i)/atmos%H(i)))
150  endif
151  endif
152  cvh2denp = ans
153  end
154 ! ----------------------------------
155  real*8 function cvh2scaleh(z)
156 ! ----------------------------------
157  implicit none
158 #include "Zatmos.h"
159  real*8 z
160  real*8 ans
161  integer i
162 
163  if( z .gt. atmos%z(atmos%nodes-1) ) then
164  ans = atmos%H(atmos%nodes-1)
165  elseif(z .lt. atmos%z(1)) then
166  ans = atmos%H(1)
167  else
168  call kdwhereis(z, atmos%nodes, atmos%z, 1, i)
169  ans = atmos%H(i) + atmos%a(i)*(z-atmos%z(i))
170  endif
171  cvh2scaleh = ans
172  end
173 ! -------------------------------------
174  real*8 function cvh2den2p(z)
175 ! -------------------------------------
176 ! d(d rho/dz)/dz
177  implicit none
178 #include "Zatmos.h"
179  real*8 z
180  real*8 ans
181  integer i
182  real*8 a
183 
184  if( z .gt. atmos%z(atmos%nodes) ) then
185  ans = 0.
186  elseif(z .lt. atmos%z(1)) then
187  ans = atmos%rho(1)/atmos%H(1)/atmos%H(1)*
188  * exp( (atmos%z(1)-z)/atmos%H(1) )
189  else
190  call kdwhereis(z, atmos%nodes, atmos%z, 1, i)
191  a = atmos%a(i)
192  if(a .ne. 0.) then
193  ans = atmos%rho(i)* (a/atmos%H(i))**2 *
194  * (1. + a*(z-atmos%z(i))/atmos%H(i))**(-3.-1./a)
195  else
196  ans = atmos%rho(i)/atmos%H(i)**2 *
197  * exp(-(z-atmos%z(i))/atmos%H(i))
198  endif
199  endif
200  cvh2den2p = ans
201  end
202 ! ---------------------------------------
203  real*8 function cvh2thick(z)
204 ! ---------------------------------------
205  implicit none
206 #include "Zatmos.h"
207  real*8 z
208 
209  real*8 ans
210 
211  integer i
212  real*8 a
213 
214  if( z .gt. atmos%z(atmos%nodes) ) then
215  ans = atmos%cumd(atmos%nodes) *
216  * exp((atmos%z(atmos%nodes) - z)/Hinf )
217  elseif(z .lt. atmos%z(1)) then
218  ans = atmos%cumd(1)*
219  * exp( (atmos%z(1)-z)/atmos%H(1) )
220  else
221  call kdwhereis(z, atmos%nodes, atmos%z, 1, i)
222  a = atmos%a(i)
223  if(a .ne. 0.) then
224  ans = atmos%cumd(i) - atmos%d0(i)*(1.-
225  * (1+ a*(z-atmos%z(i))/atmos%H(i))**(-1/a))
226  else
227  ans = atmos%cumd(i) - atmos%d0(i)*(1.-
228  * exp(-(z-atmos%z(i))/atmos%H(i)))
229  endif
230  endif
231  cvh2thick = ans
232  end
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
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
kdwhereis
subroutine kdwhereis(x, in, a, step, loc)
Definition: kdwhereis.f:27
dod0
nodes dod0
Definition: cstdatmos0.fNew2.h:84
data
block data cblkMuInt integer i ! Pair total X sec data(MuPrTX(i), i=1, 38)/1 0.431132E-02
ne
int ne[nl][nth]
Definition: Zprivate.h:11
nrl_atmos::temp
real(4), dimension(:), allocatable, save temp
Definition: cNRLAtmos.f:29
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
d0
block data cblkEvhnp ! currently usable models data RegMdls ad *special data *Cekaon d0
Definition: cblkEvhnp.h:5
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
cthick2h
nodes a atmos atmos temp real * cthick2h
Definition: cstdatmos0.fNew2.h:107
d
dE dx *! Nuc Int sampling table d
Definition: cblkMuInt.h:130
T
! type magfield sequence ! Note that position vector where the magnetic field is given ! is not included here ! unit of field strength is in T(1 gauss=10 **-4 T) union map real *8 x
a
nodes a
Definition: cstdatmos0.fNew2.h:93
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
cvh2den
real *8 function cvh2den(z)
Definition: ciniSegAtoms.f:54
b
dE dx *! Nuc Int sampling table b
Definition: cblkMuInt.h:130
t
nodes t
Definition: cstdatmos0.fNew2.h:84
H
nodes a atmos H(i)
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
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