1 | !=================================================================================== |
---|
2 | subroutine photon(iele,rad,dlength,energy,ieave,iquasto,d1,d2) |
---|
3 | !SUBROUTINE photon (i_elem_type, rad_curv_m, length_curr_elem, & |
---|
4 | ! Energy_beam_MeV, mass_rest_MeV, d1, d2) |
---|
5 | !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
---|
6 | ! |
---|
7 | ! subroutine to determine the energy losses at entrance |
---|
8 | ! and exit of magnet (dipole or quadrupole) |
---|
9 | ! |
---|
10 | ! |
---|
11 | ! input: |
---|
12 | ! iele - Magnet Type, iele=2 means a Quadrupole and |
---|
13 | ! its Radiation can be switched off with iquasto flag |
---|
14 | ! set in the Initialisation |
---|
15 | ! rad - radius of curvature |
---|
16 | ! (rho for bend |
---|
17 | ! 1/(abs(k1)*sqrt(x**2+y**2)) for quadrupole |
---|
18 | ! |
---|
19 | ! dlength - length of element |
---|
20 | ! energy - beam energy in MeV |
---|
21 | ! output: |
---|
22 | ! d1 - momentum loss at entrance (normalized to beam energy) |
---|
23 | ! d2 - momentum loss at exit (normalized to beam energy) |
---|
24 | ! |
---|
25 | ! FZ/28.02.2001 |
---|
26 | !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
---|
27 | USE precision_constants, ONLY:CLIGHT,dhbar,pmae,dp,sp,qelect,c_1d6,c_1d_10,c_1d_3,c_0_3079,c_137,one,two,three,five |
---|
28 | implicit none |
---|
29 | include 'photoni.inc' |
---|
30 | |
---|
31 | INTEGER, INTENT(IN) :: iele, ieave, iquasto |
---|
32 | REAL(dp), INTENT(IN) :: dlength,energy |
---|
33 | REAL(dp), INTENT(OUT) :: d1,d2 |
---|
34 | REAL(dp), INTENT(INOUT) :: rad |
---|
35 | |
---|
36 | REAL(dp) :: gamma,theta, dnpho, dnphohalf,ec,eave,ep1,ep2, ra2, & |
---|
37 | ran,xl,en,gamfac |
---|
38 | |
---|
39 | REAL(dp), parameter :: el=qelect, vl=CLIGHT |
---|
40 | real(sp) amu |
---|
41 | Double Precision :: ran2 ! function |
---|
42 | ! integer i,j,nmax,n1,n2,npmax,iseed,idumy,ieave,iquasto,iele |
---|
43 | ! double precision pi,vl,dhbar,el,gamma,energy,xilog,sollog,b,sol, & |
---|
44 | ! &xil,sll,ak,bk,theta,dlength,rad,dnpho,dnphohalf,ec,eave,ep1,ep2, & |
---|
45 | ! &d1,d2,ra2,ran,ran2,xl,en,gamfac |
---|
46 | ! parameter (nmax=100) |
---|
47 | |
---|
48 | ! common / lookup / xil(nmax), sll(nmax) |
---|
49 | ! common / lookup2 / ak(nmax), bk(nmax) |
---|
50 | ! common / rann / iseed, idumy |
---|
51 | ! common / trick / pi, vl, dhbar, el, gamma |
---|
52 | ! common / eave1 / ieave,iquasto |
---|
53 | |
---|
54 | ! real amu, amax |
---|
55 | ! integer ierr,nphoton |
---|
56 | INTEGER :: ierr, nphoton, n1, n2, npmax, i,j, idumy |
---|
57 | external rnpset, rnpssn |
---|
58 | ierr=0 |
---|
59 | rad=abs(rad) |
---|
60 | ! print*,"RRRRRRRRRRRRRRRRRad in photon",iele,rad,dlength,energy |
---|
61 | |
---|
62 | ! for electrons |
---|
63 | gamma=energy/pmae*c_1d_3 |
---|
64 | |
---|
65 | if (gamma.le.1) then |
---|
66 | write(*,*) ' error in subroutine photon - no initialization ' |
---|
67 | return |
---|
68 | endif |
---|
69 | |
---|
70 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
---|
71 | ! |
---|
72 | ! Stop Quadrupole Radiation |
---|
73 | ! |
---|
74 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
---|
75 | |
---|
76 | if (iquasto.eq.0.and.iele.eq.2) return |
---|
77 | ! it is check already before CALL this subr. |
---|
78 | |
---|
79 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
---|
80 | ! |
---|
81 | ! compute mean number of photons emitted |
---|
82 | ! |
---|
83 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
---|
84 | |
---|
85 | 989 continue |
---|
86 | |
---|
87 | ! trajectory bending angle |
---|
88 | theta = dlength/rad |
---|
89 | |
---|
90 | dnpho = five/two/sqrt(three)*gamma/c_137*theta |
---|
91 | ! |
---|
92 | dnphohalf = dnpho/two |
---|
93 | ! |
---|
94 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
---|
95 | ! |
---|
96 | ! compute critical energy and average (half) energy loss |
---|
97 | ! |
---|
98 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
---|
99 | ! |
---|
100 | ! critical photon energy as fraction of beam energy |
---|
101 | |
---|
102 | ec = three/two*vl*gamma**3/rad*dhbar/(energy*el*c_1d6) |
---|
103 | ! |
---|
104 | eave = dnphohalf*c_0_3079*ec |
---|
105 | ! |
---|
106 | !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
---|
107 | ! |
---|
108 | ! compute number of photons n1 and n2 actually emitted |
---|
109 | ! |
---|
110 | !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
---|
111 | ! |
---|
112 | amu = dnphohalf |
---|
113 | ! write(*,*) ' call rnpssn 1 ',amu,ierr |
---|
114 | call rnpssn(amu,nphoton,ierr) |
---|
115 | ! write(*,*) ' call rnpssn 2 ',amu,ierr,nphoton |
---|
116 | |
---|
117 | n1 = nphoton |
---|
118 | call rnpssn(amu,nphoton,ierr) |
---|
119 | ! write(*,*) ' after calling rnpssn ',amu,ierr,nphoton |
---|
120 | n2 = nphoton |
---|
121 | |
---|
122 | npmax = n1 + n2 |
---|
123 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
---|
124 | ! |
---|
125 | ! toss randum number to determine photon energy |
---|
126 | ! |
---|
127 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
---|
128 | ep1 = 0. |
---|
129 | ep2 = 0. |
---|
130 | |
---|
131 | do i = 1, npmax |
---|
132 | |
---|
133 | ra2 = ran2(idumy) |
---|
134 | if (ra2.lt.c_1d_10) goto 989 |
---|
135 | |
---|
136 | ran = log(ra2) |
---|
137 | do j = 1, nmax |
---|
138 | if (ran.lt.sll(j)) then |
---|
139 | if(j.eq.1) then |
---|
140 | xl = (ran-bk(1))/ak(1) |
---|
141 | goto 990 |
---|
142 | endif |
---|
143 | xl = (ran-bk(j))/ak(j) |
---|
144 | goto 990 |
---|
145 | endif |
---|
146 | enddo |
---|
147 | xl = (ran-bk(nmax))/ak(nmax) |
---|
148 | 990 continue |
---|
149 | |
---|
150 | ! |
---|
151 | ! note: xl is logarithm of photon energy |
---|
152 | ! divided by critical energy |
---|
153 | ! |
---|
154 | |
---|
155 | if (xl.lt.(-46.).or.(xl.gt.(2.3))) goto 989 |
---|
156 | |
---|
157 | ! photon energy as fraction of beam energy |
---|
158 | en = -exp(xl)*ec |
---|
159 | |
---|
160 | if (i.le.n1) then |
---|
161 | ep1 = ep1 + en |
---|
162 | else |
---|
163 | ep2 = ep2 + en |
---|
164 | endif |
---|
165 | |
---|
166 | enddo |
---|
167 | |
---|
168 | ! compute total relative photon energy (average loss is |
---|
169 | ! added for compensation) |
---|
170 | ! |
---|
171 | gamfac = gamma*gamma/(gamma*gamma-one) |
---|
172 | if( ieave.eq.1 ) then |
---|
173 | d1 = (ep1 + eave) * gamfac |
---|
174 | d2 = (ep2 + eave) * gamfac |
---|
175 | else |
---|
176 | d1 = ep1 * gamfac |
---|
177 | d2 = ep2 * gamfac |
---|
178 | endif |
---|
179 | |
---|
180 | return |
---|
181 | end subroutine photon |
---|
182 | !=================================================================================== |
---|
183 | |
---|
184 | |
---|
185 | !=================================================================================== |
---|
186 | !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
---|
187 | !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
---|
188 | double precision function ran2(idum) |
---|
189 | !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
---|
190 | ! |
---|
191 | ! The Function Returns a Uniform Random Number between 0.0 and 1.0. |
---|
192 | ! |
---|
193 | ! The Function uses a 'Subtractive Method'!!! |
---|
194 | ! |
---|
195 | ! Set 'idum' to a negative value to initialize or reinitialize the Sequence. |
---|
196 | ! |
---|
197 | ! (Numerical Recipies ran3(idumy), pg. 273.) |
---|
198 | ! |
---|
199 | !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
---|
200 | USE precision_constants, ONLY:one |
---|
201 | implicit none |
---|
202 | integer idum |
---|
203 | integer mbig, mseed, mz |
---|
204 | double precision fac |
---|
205 | parameter(mbig=1000000000,mseed=161803398,mz=0,fac=one/mbig) |
---|
206 | integer i, iff, ii, inext, inextp, k |
---|
207 | integer mj, mk, ma(55) |
---|
208 | save iff, inext, inextp, ma |
---|
209 | data iff /0/ |
---|
210 | |
---|
211 | ! Initialization: |
---|
212 | if(idum.lt.0.or.iff.eq.0) then |
---|
213 | iff = 1 |
---|
214 | mj = mseed-iabs(idum) |
---|
215 | mj = mod(mj,mbig) |
---|
216 | mk = 1 |
---|
217 | do i=1,54 |
---|
218 | ii = mod(2*i,55) |
---|
219 | ma(ii) = mk |
---|
220 | mk = mj - mk |
---|
221 | if(mk.lt.mz) mk = mk + mbig |
---|
222 | mj = ma(ii) |
---|
223 | enddo |
---|
224 | |
---|
225 | do k=1,4 |
---|
226 | do i=1,55 |
---|
227 | ma(i) = ma(i) - ma(1+mod(i+30,55)) |
---|
228 | if(ma(i).lt.mz) ma(i) = ma(i) + mbig |
---|
229 | enddo |
---|
230 | enddo |
---|
231 | inext = 0 |
---|
232 | inextp = 31 |
---|
233 | idum = 1 |
---|
234 | endif |
---|
235 | |
---|
236 | inext = inext + 1 |
---|
237 | if(inext.eq.56) inext = 1 |
---|
238 | inextp = inextp + 1 |
---|
239 | if(inextp.eq.56) inextp = 1 |
---|
240 | mj = ma(inext) - ma(inextp) |
---|
241 | if(mj.lt.mz) mj = mj + mbig |
---|
242 | ma(inext) = mj |
---|
243 | ran2 = mj * fac |
---|
244 | |
---|
245 | return |
---|
246 | END function ran2 |
---|
247 | !======================================================================================== |
---|