[430] | 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 | !======================================================================================== |
---|