source: PSPA/madxPSPA/src/trrun_phot.f90 @ 442

Last change on this file since 442 was 430, checked in by touze, 11 years ago

import madx-5.01.00
File size: 8.7 KB
Line 
1subroutine trphot(el,curv,rfac,deltap)
2
3  implicit none
4
5  !----------------------------------------------------------------------*
6  ! Purpose:                                                             *
7  !   Generate random energy loss for photons, using a look-up table to  *
8  !   invert the function Y.  Ultra-basic interpolation computed;        *
9  !   leads to an extrapolation outside the table using the two outmost  *
10  !   point on each side (low and high).                                 *
11  !   Assumes ultra-relativistic particles (beta = 1).                   *
12  ! Author: Ghislain Roy                                                 *
13  ! Input:                                                               *
14  !   EL     (double)       Element length.                              *
15  !   CURV   (double)       Local curvature of orbit.                    *
16  ! Output:                                                              *
17  !   RFAC   (double)       Relative energy loss due to photon emissions.*
18  !----------------------------------------------------------------------*
19  !---- Generate pseudo-random integers in batches of NR.                *
20  !     The random integers are generated in the range [0, MAXRAN).      *
21  !---- Table definition: maxtab, taby(maxtab), tabxi(maxtab)            *
22  !----------------------------------------------------------------------*
23  integer i,ierror,j,nphot,nr,maxran,maxtab
24  parameter(nr=55,maxran=1000000000,maxtab=101)
25  double precision amean,curv,dlogr,el,frndm,rfac,scalen,scaleu,    &
26       slope,ucrit,xi,deltap,hbar,clight,arad,pc,gamma,amass,real_am,    &
27       get_value,get_variable,tabxi(maxtab),taby(maxtab),zero,one,two,   &
28       three,five,twelve,fac1
29  parameter(zero=0d0,one=1d0,two=2d0,three=3d0,five=5d0,twelve=12d0,&
30       fac1=3.256223d0)
31  character(20) text
32  data (taby(i), i = 1, 52)                                         &
33       / -1.14084005d0,  -0.903336763d0, -0.769135833d0, -0.601840854d0, &
34       -0.448812515d0, -0.345502228d0, -0.267485678d0, -0.204837948d0,   &
35       -0.107647471d0, -0.022640628d0,  0.044112321d0,  0.0842842236d0,  &
36       0.132941082d0,  0.169244036d0,  0.196492359d0,  0.230918407d0,    &
37       0.261785239d0,  0.289741248d0,  0.322174788d0,  0.351361096d0,    &
38       0.383441716d0,  0.412283719d0,  0.442963421d0,  0.472622454d0,    &
39       0.503019691d0,  0.53197819d0,   0.561058342d0,  0.588547111d0,    &
40       0.613393188d0,  0.636027336d0,  0.675921738d0,  0.710166812d0,    &
41       0.725589216d0,  0.753636241d0,  0.778558254d0,  0.811260045d0,    &
42       0.830520391d0,  0.856329501d0,  0.879087269d0,  0.905612588d0,    &
43       0.928626955d0,  0.948813677d0,  0.970829248d0,  0.989941061d0,    &
44       1.0097903d0,    1.02691281d0,   1.04411256d0,   1.06082714d0,     &
45       1.0750246d0,    1.08283985d0,   1.0899564d0,    1.09645379d0 /
46  data (taby(i), i = 53, 101)                                       &
47       /  1.10352755d0,   1.11475027d0,   1.12564385d0,   1.1306442d0,   &
48       1.13513422d0,   1.13971806d0,   1.14379156d0,   1.14741969d0,     &
49       1.15103698d0,   1.15455759d0,   1.15733826d0,   1.16005647d0,     &
50       1.16287541d0,   1.16509759d0,   1.16718769d0,   1.16911888d0,     &
51       1.17075884d0,   1.17225218d0,   1.17350936d0,   1.17428589d0,     &
52       1.17558432d0,   1.17660713d0,   1.17741513d0,   1.17805469d0,     &
53       1.17856193d0,   1.17896497d0,   1.17928565d0,   1.17954147d0,     &
54       1.17983139d0,   1.1799767d0,    1.18014216d0,   1.18026078d0,     &
55       1.18034601d0,   1.1804074d0,    1.18045175d0,   1.1804837d0,      &
56       1.18051291d0,   1.18053186d0,   1.18054426d0,   1.18055236d0,     &
57       1.18055761d0,   1.18056166d0,   1.18056381d0,   1.1805656d0,      &
58       1.18056655d0,   1.18056703d0,   1.18056726d0,   1.1805675d0,      &
59       1.18056762d0 /
60  data (tabxi(i), i = 1, 52)                                        &
61       / -7.60090017d0,  -6.90775537d0,  -6.50229025d0,  -5.99146461d0,  &
62       -5.52146101d0,  -5.20300722d0,  -4.96184492d0,  -4.76768923d0,    &
63       -4.46540833d0,  -4.19970512d0,  -3.98998451d0,  -3.86323285d0,    &
64       -3.70908213d0,  -3.59356928d0,  -3.50655794d0,  -3.39620972d0,    &
65       -3.29683733d0,  -3.20645332d0,  -3.10109282d0,  -3.0057826d0,     &
66       -2.9004221d0,   -2.80511189d0,  -2.70306253d0,  -2.60369015d0,    &
67       -2.50103593d0,  -2.4024055d0 ,  -2.30258512d0,  -2.20727491d0,    &
68       -2.12026358d0,  -2.04022098d0,  -1.89712d0   ,  -1.7719568d0,     &
69       -1.71479833d0,  -1.60943794d0,  -1.51412773d0,  -1.38629436d0,    &
70       -1.30933332d0,  -1.20397282d0,  -1.10866261d0,  -0.99425226d0,    &
71       -0.89159810d0,  -0.79850775d0,  -0.69314718d0,  -0.59783697d0,    &
72       -0.49429631d0,  -0.40047753d0,  -0.30110508d0,  -0.19845095d0,    &
73       -0.10536054d0,  -0.05129330d0,   0.0d0,          0.048790119d0 /
74  data (tabxi(i), i = 53, 101)                                      &
75       /  0.104360029d0,  0.198850885d0,  0.300104618d0,  0.350656837d0, &
76       0.398776114d0,  0.451075643d0,  0.500775278d0,  0.548121393d0,    &
77       0.598836541d0,  0.652325153d0,  0.69813472d0 ,  0.746687889d0,    &
78       0.802001595d0,  0.850150883d0,  0.900161386d0,  0.951657832d0,    &
79       1.00063193d0,   1.05082154d0,   1.09861231d0,   1.13140213d0,     &
80       1.1939224d0,    1.25276291d0,   1.3083328d0,    1.36097658d0,     &
81       1.4109869d0,    1.45861506d0,   1.50407743d0,   1.54756248d0,     &
82       1.60943794d0,   1.64865863d0,   1.70474803d0,   1.75785792d0,     &
83       1.80828881d0,   1.85629797d0,   1.90210748d0,   1.9459101d0,      &
84       2.0014801d0,    2.05412364d0,   2.10413408d0,   2.15176225d0,     &
85       2.19722462d0,   2.25129175d0,   2.29253483d0,   2.35137534d0,     &
86       2.40694523d0,   2.45100522d0,   2.501436d0,     2.60268974d0,     &
87       2.64617491d0 /
88
89  !Get constants
90  clight = get_variable('clight ')
91  hbar   = get_variable('hbar ')
92  arad   = get_value('probe ','arad ')
93  pc     = get_value('probe ','pc ')
94  amass  = get_value('probe ','mass ')
95  gamma  = get_value('probe ','gamma ')
96  deltap = get_value('probe ','deltap ')
97  scalen = five / (twelve * hbar * clight)
98  scaleu = hbar * three * clight / two
99
100  !---- AMEAN is the average number of photons emitted.,
101  !     NPHOT is the integer number generated from Poisson's law.
102  amean = scalen * abs(arad*pc*(one+deltap)*el*curv) * sqrt(three)
103  rfac = zero
104  real_am = amean
105  if (real_am .gt. zero) then
106     call dpoissn(real_am, nphot, ierror)
107
108     if (ierror .ne. 0) then
109        write(text, '(1p,d20.12)') amean
110        call aafail('TRPHOT: ','Fatal: Poisson input mean =' // text)
111     endif
112
113     !---- For all photons, sum the radiated photon energy,
114     !     in units of UCRIT (relative to total energy).
115
116     if (nphot .ne. 0) then
117        ucrit = scaleu * gamma**2 * abs(curv) / amass
118        xi = zero
119        do i = 1, nphot
120
121           !---- Find a uniform random number in the range [ 0,3.256223 ].
122           !     Note that the upper limit is not exactly 15*sqrt(3)/8
123           !     because of imprecision in the integration of F.
124           dlogr = log(fac1 * frndm())
125
126           !---- Now look for the energy of the photon in the table TABY/TABXI
127           do j = 2, maxtab
128              if (dlogr .le. taby(j) ) go to 20
129           enddo
130
131           !---- Perform linear interpolation and sum up energy lost.
13220         slope = (dlogr - taby(j-1)) / (taby(j) - taby(j-1))
133           xi = dexp(tabxi(j-1) + slope * (tabxi(j) - tabxi(j-1)))
134           rfac = rfac + ucrit * xi
135        enddo
136     endif
137  endif
138end subroutine trphot
139subroutine dpoissn (amu,n,ierror)
140
141  implicit none
142
143  !----------------------------------------------------------------------*
144  !    POISSON GENERATOR                                                 *
145  !    CODED FROM LOS ALAMOS REPORT      LA-5061-MS                      *
146  !    PROB(N)=EXP(-AMU)*AMU**N/FACT(N)                                  *
147  !        WHERE FACT(N) STANDS FOR FACTORIAL OF N                       *
148  !    ON RETURN IERROR.EQ.0 NORMALLY                                    *
149  !              IERROR.EQ.1 IF AMU.LE.0.                                *
150  !----------------------------------------------------------------------*
151  integer n, ierror
152  double precision amu,ran,pir,frndm,grndm,expma,amax,zero,one,half
153  parameter(zero=0d0,one=1d0,half=5d-1)
154  !    AMAX IS THE VALUE ABOVE WHICH THE NORMAL DISTRIBUTION MUST BE USED
155  data amax/88d0/
156
157  ierror= 0
158  if(amu.le.zero) then
159     !    MEAN SHOULD BE POSITIVE
160     ierror=1
161     n = 0
162     go to 999
163  endif
164  if(amu.gt.amax) then
165     !   NORMAL APPROXIMATION FOR AMU.GT.AMAX
166     ran = grndm()
167     n=ran*sqrt(amu)+amu+half
168     goto 999
169  endif
170  expma=exp(-amu)
171  pir=one
172  n=-1
17310 n=n+1
174  pir=pir*frndm()
175  if(pir.gt.expma) go to 10
176999 end subroutine dpoissn
Note: See TracBrowser for help on using the repository browser.