COSMOS v7.655  COSMOSv7655
(AirShowerMC)
cmkPrimSTbl.f
Go to the documentation of this file.
1 ! *******************************************************
2  subroutine cmkprimstbl(each, icon)
3 ! make primary energy sampling table for a given
4 ! component.
5 ! *******************************************************
6  implicit none
7 
8 #include "Zptcl.h"
9 #include "Zprimary.h"
10 #include "Zprimaryc.h"
11  type(component):: each
12  integer icon
13 !
14  integer k, nseg, l
15 
16  character*150 msg
17  integer firstseg, i, lastseg
18  real*8 cut, cut2
19 !
20  nseg = each%no_of_seg ! # of segmets = datapoints -1
21  cut = each%cut
22  cut2 =each%cut2
23 ! first convert flux into normal expression
24 ! energy unit is kept as input
25 ! original each flux = (dI/dE)*E**flatterer
26 !
27 ! check energy order.
28  do l = 1, nseg
29  if(each%energy(l) .ge. each%energy(l+1)) then
30  write(*,*) each%label,'-th component=',each%symb,
31  * ' energy is not in ascending order.'
32  write(*,*) ' It is in the ', l, '-th segment'
33  endif
34  enddo
35 ! flux without multiplication of E**flatterer
36  do l = 1, nseg+1
37  each%flux(l) = each%flux(l)/each%energy(l)**each%flatterer
38  enddo
39  firstseg = 1
40  lastseg = nseg
41 
42 ! now flux = dI/dE = const *E**(-beta). get beta next
43  if(each%diff_or_inte .eq. 'd') then
44 ! get beta in E**(-beta)dE in each segment
45  do k =1, nseg
46  each%beta(k) = - log(each%flux(k)/each%flux(k+1))
47  * /
48  * log(each%energy(k)/each%energy(k+1))
49  enddo
50 ! assume last segment is the end of spectrum
51  if(nseg .ge. 1) each%beta(nseg+1) = 100.
52 !
53 ! -----------------
54 ! find segment i <= cut < i+1, if cut > 0
55 ! and adjust each.emin.
56 !
57  if(cut .gt. 0. .and. nseg .ge. 1 ) then
58  if(cut .ge. each%energy(nseg+1)) then
59  call cerrormsg(
60  * 'the primary minmum cut >= table max',0)
61  endif
62  do i = nseg, 1, -1
63  if(cut .ge. each%energy(i)) then
64  firstseg = i
65  each%emin = cut
66  goto 100
67  endif
68  enddo
69  each%emin = each%energy(1)
70  firstseg = 1
71 ! nothing to do
72  goto 200
73  100 continue
74 ! flux at cut
75  k = firstseg
76  each%flux(k) = each%flux(k)*
77  * ( cut/each%energy(k))**(-each%beta(k))
78  each%energy(k) = cut
79  endif
80  200 continue
81 ! ------------------
82 ! find segment i < cut2 <= i+1, if cut2 > 0
83 ! and adjust each.emax,
84 !
85 
86 
87  if(cut2 .gt. 0. .and. nseg .ge. 1 ) then
88  if(cut2 .le. each%energy(1)) then
89  call cerrormsg(
90  * 'the primary max cut <= table min',0)
91  endif
92  do i = 1, nseg
93  if(cut2 .le. each%energy(i+1)) then
94  lastseg = i
95  each%emax= cut2
96  goto 300
97  endif
98  enddo
99  each%emax = each%energy(nseg+1)
100  lastseg = nseg
101 ! nothing to do
102  goto 400
103  300 continue
104 ! flux at cut2
105  k = lastseg
106  each%flux(k+1) = each%flux(k)*
107  * ( cut2/each%energy(k))**(-each%beta(k))
108  each%energy(k+1) = cut2
109  endif
110  400 continue
111 
112 ! get integral distribution
113  if(nseg .ge. 1) then
114 ! assume spectrum is cut at each.emax
115  each%norm_inte(lastseg+1)= 0.
116  do k = nseg+1, lastseg+1, -1
117  each%norm_inte(k) = 0.
118  enddo
119  do k = lastseg, firstseg, -1
120  if(abs(each%beta(k)-1.d0) .gt. 1.d-6) then
121  each%norm_inte(k) = each%norm_inte(k+1) +
122  * ( each%energy(k+1)*each%flux(k+1)
123  * -each%energy(k)*each%flux(k) )/(1.d0-each%beta(k))
124  else
125  each%norm_inte(k) = each%norm_inte(k+1) +
126  * each%flux(k)* each%energy(k)*
127  * log(each%energy(k+1)/each%energy(k))
128  endif
129  enddo
130  endif
131  elseif(each%diff_or_inte .eq. 'i') then
132 ! integral. first make power of integral spectrum
133  do k =1, nseg
134  each%beta(k) = - log(each%flux(k)/
135  * each%flux(k+1))
136  * /
137  * log(each%energy(k)/each%energy(k+1))
138  enddo
139  if(nseg .ge. 1) then
140  each%beta(nseg + 1) = each%beta(nseg)
141  endif
142 ! move flux into norm_inte
143  do k = 1, nseg+1
144  each%norm_inte(k) = each%flux(k)
145  enddo
146  else
147  write(*,*) ' error specification of diff/integral=',
148  * each%diff_or_inte
149  icon = 1
150  goto 900
151  endif
152 
153 ! total integral value
154  if(nseg .ge. 1) then
155  each%inte_value = each%norm_inte(firstseg)
156  else
157 ! delta funcition
158  each%inte_value = each%flux(1)
159  endif
160 ! normalize the flux
161  do k = firstseg, lastseg
162  each%norm_inte(k) = each%norm_inte(k)/each%inte_value
163  enddo
164  900 continue
165  end
166 
167 
168 
169 
170 
cerrormsg
subroutine cerrormsg(msg, needrtn)
Definition: cerrorMsg.f:4
d0
block data cblkEvhnp ! currently usable models data RegMdls ad *special data *Cekaon d0
Definition: cblkEvhnp.h:5
cmkprimstbl
subroutine cmkprimstbl(each, icon)
Definition: cmkPrimSTbl.f:3
d
dE dx *! Nuc Int sampling table d
Definition: cblkMuInt.h:130