COSMOS v7.655  COSMOSv7655
(AirShowerMC)
cmBremI1.f
Go to the documentation of this file.
1 ! program testBremI1
2 ! implicit none
3 ! real*8 upsilon, xc, y, cmBremI1
4 ! xc = 0.
5 !c do while (xc .lt. 4.e-3)
6 ! upsilon = 100.
7 ! do while (upsilon .gt. 1.e-6)
8 ! y = cmBremI1(upsilon, xc)
9 ! write(*, *) sngl(upsilon), sngl(y)
10 ! upsilon = upsilon / 10.**0.1
11 ! if(y .lt. 1.e-5) upsilon = 1.e-9
12 ! enddo
13 ! write(*, *)
14 !c
15 !c enddo
16 ! end
17 ! **********************************************************
18  real*8 function cmbremi1(upsilon, xc)
19  implicit none
20  real*8 upsilon ! input. Upsilon value
21  real*8 xc ! input. fix the cut off, below which the synchroron
22  ! is treated as energy loss only.
23  ! one of 3.1627 x 10-3, 10-3, ... 10-6, or 0.0
24  ! if xc is not one of these, neareset one
25  ! is chosen.
26 
27 ! some values;
28 ! U 31.6 10 1 .1 .01
29 ! I1 1.99 3.6 4.85 5.19
30 ! I2 .19 .227 .15 1.76e-2 3.53e-4
31 !
32 ! Psi0: the integral x= Xc to 1 of 4 k(z)/z * 1/(2+ 3Uz)^2 dz
33 ! where z = 2/3/U* x/(1-x).
34 !
35 ! Xc (=Eg/Ee) must be one
36 ! 3.162x10-3, 10-3, 3.162x10-4, 10-4,... 10-6, 0.
37 !
38 ! The 'cross-section' for emitting fractional gamma ray energy x ~ x+dx,
39 ! in the unit distance is given by
40 ! p(Ee, x, U)dx = root(3)/2/pi (SyncConvR)/Ee* U * (1-x)/x K(2zeta)dx
41 ! (/meter, if Ee is in GeV). ( U is E/m * B/Bc ).
42 ! This is integrated from x = Xc to 1 to give the total 'crosssection'
43 ! The following tbl containes 11 coefficients for the polynomial
44 ! approximation of the part of the integral as a function
45 ! of Upsilon. (in log-log)
46 !
47 ! polynomial value times root(3)/2/pi SyncConvR/Ee *U
48 ! is the total cross-section (= probability/m).
49 !
50 ! tbl(*, i) i=1, 8 contains table of 11 coefficients of polynaomials
51 ! Integral (x =Xc to 1 ) of (1-x)/x K( 2/3/U x/(1-x))/ 4 dx
52 ! where U= upsilon value. Xc= 3.1627 x 10-3 to 10-6 step /3.1627
53 ! (i=1 is for 3.162x10-3, i=2 for 10-3, ...).
54 !
55  real*8 psi, x, xc0coef(8)
56  integer ixc, i
57 !
58  real*8 tbl(11, 8)
59  real*8 psi3_3(11), psi_3(11), psi3_4(11), psi_4(11)
60  real*8 psi3_5(11), psi_5(11), psi3_6(11), psi_6(11)
61  equivalence(psi3_3(1), tbl(1,1))
62  equivalence(psi_3(1), tbl(1,2))
63  equivalence(psi3_4(1), tbl(1,3))
64  equivalence(psi_4(1), tbl(1,4))
65  equivalence(psi3_5(1), tbl(1,5))
66  equivalence(psi_5(1), tbl(1,6))
67  equivalence(psi3_6(1), tbl(1,7))
68  equivalence(psi_6(1), tbl(1,8))
69 
70 
71 
72 
73  data (psi3_3(i), i= 1, 11)/
74  1 -0.36519777 , -0.16437488 , -0.42191193e-01,
75  2 0.45068651e-02, 0.37406044e-03, -0.11694408e-03,
76  3 -0.21981366e-04, 0.50330259e-05, 0.56416748e-06,
77  4 -0.59483304e-07, -0.10075934e-07
78  * /
79  data (psi_3(i), i= 1, 11)/
80  1 -0.27621962 , -0.17806948 , -0.36641237e-01,
81  2 0.34820468e-02, 0.23404694e-03, -0.11645243e-03,
82  3 -0.11166093e-05, 0.51644288e-05, 0.66054175e-07,
83  4 -0.94357089e-07, -0.73799132e-08
84  * /
85  data (psi3_4(i), i= 1, 11)/
86  1 -0.21981296 , -0.18463515 , -0.33139695e-01,
87  2 0.22214970e-02, 0.96012771e-04, -0.47834498e-04,
88  3 0.15951965e-04, 0.32866129e-05, -0.42901554e-06,
89  4 -0.10473240e-06, -0.53301833e-08
90  * /
91  data (psi_4(i), i= 1, 11)/
92  1 -0.18243710 , -0.18712415 , -0.31745835e-01,
93  2 0.90498800e-03, 0.12296262e-03, 0.46807678e-04,
94  3 0.22577025e-04, 0.42917421e-06, -0.75768289e-06,
95  4 -0.97835715e-07, -0.37327138e-08
96  * /
97 
98  data (psi3_5(i), i= 1, 11)/
99  1 -0.15645188 , -0.18758306 , -0.32142658e-01,
100  2 -0.27076403e-03, 0.30784765e-03, 0.14090033e-03,
101  3 0.21276677e-04, -0.25452860e-05, -0.96816286e-06,
102  4 -0.87676147e-07, -0.26691205e-08
103  * /
104  data (psi_5(i), i= 1, 11)/
105  1 -0.13939747 , -0.19218187 , -0.32177166e-01,
106  2 0.43072437e-03, 0.50884199e-03, 0.78017665e-04,
107  3 0.33119508e-05, -0.20977899e-05, -0.47432646e-06,
108  4 -0.37755362e-07, -0.10680872e-08
109  * /
110  data (psi3_6(i), i= 1, 11)/
111  1 -0.12864838 , -0.19609994 , -0.31528914e-01,
112  2 0.11811036e-02, 0.56385559e-03, 0.84414004e-05,
113  3 -0.73982020e-05, -0.91893847e-06, -0.85709610e-07,
114  4 -0.61186589e-08, -0.19199417e-09
115  * /
116  data (psi_6(i), i= 1, 11)/
117  1 -0.12129256 , -0.19786006 , -0.31152534e-01,
118  2 0.13900376e-02, 0.58057791e-03, -0.13423486e-04,
119  3 -0.10636349e-04, -0.54026892e-06, 0.34550944e-07,
120  4 0.35276755e-08, 0.72269160e-10
121  * /
122 ! for Xc=0. for U>3.e-3.
123  data ( xc0coef(i), i= 1, 8)/
124  1 -0.10535316 , -0.20019514 , -0.30586898e-01,
125  2 0.13258432e-02, 0.58454445e-03, -0.13505923e-04,
126  3 -0.97296726e-05, -0.53504562e-06
127  * /
128 
129  if(xc .eq. 0.) then
130  if(upsilon .lt. 3.e-3) then
131  psi = 0.269
132  else
133  psi = xc0coef(8)
134  x = log(upsilon)
135  do i = 7, 1, -1
136  psi = x * psi + xc0coef(i)
137  enddo
138  endif
139  else
140 !
141  ixc = log10(3.1627e-3/ xc)*2 + 1.1
142  if(abs(10.**(-2 - ixc/2.0) / xc -1.) .gt. 0.05) then
143  call cerrormsg('input xc is invalid: cmBremI1', 0)
144  endif
145  if(ixc .lt. 1 .or. ixc .gt. 8) then
146  call cerrormsg
147  * ('input xc is too large or too small: cmBremI1', 0)
148  endif
149 !
150 
151  psi = tbl(11, ixc)
152  x = log(upsilon)
153  do i = 10, 1, -1
154  psi = x * psi + tbl(i, ixc)
155  enddo
156  endif
157 
158  cmbremi1 =4* exp(psi)
159  end
160 
cerrormsg
subroutine cerrormsg(msg, needrtn)
Definition: cerrorMsg.f:4
cmbremi1
real *8 function cmbremi1(upsilon, xc)
Definition: cmBremI1.f:19
e
dE dx *! Nuc Int sampling table e
Definition: cblkMuInt.h:130