1 | subroutine ttbb(track,ktrack,parvec) |
---|
2 | |
---|
3 | implicit none |
---|
4 | |
---|
5 | !----------------------------------------------------------------------* |
---|
6 | ! purpose: * |
---|
7 | ! track a set of particle through a beam-beam interaction region. * |
---|
8 | ! see mad physicist's manual for the formulas used. * |
---|
9 | !input: * |
---|
10 | ! input/output: * |
---|
11 | ! track(6,*)(double) track coordinates: (x, px, y, py, t, pt). * |
---|
12 | ! ktrack (integer) number of tracks. * |
---|
13 | !----------------------------------------------------------------------* |
---|
14 | integer ktrack,beamshape,b_dir_int |
---|
15 | double precision track(6,*),parvec(*),fk,dp |
---|
16 | double precision gamma0,beta0,beta_dp,ptot,b_dir |
---|
17 | double precision q,q_prime,get_value,node_value,get_variable |
---|
18 | double precision zero,one,two |
---|
19 | logical first |
---|
20 | save first |
---|
21 | data first / .true. / |
---|
22 | parameter(zero=0.0d0,one=1.0d0,two=2.0d0) |
---|
23 | ! if x > explim, exp(-x) is outside machine limits. |
---|
24 | |
---|
25 | !---- Calculate momentum deviation and according changes |
---|
26 | ! of the relativistic factor beta0 |
---|
27 | dp = get_variable('track_deltap ') |
---|
28 | q = get_value('probe ','charge ') |
---|
29 | q_prime = node_value('charge ') |
---|
30 | gamma0 = parvec(7) |
---|
31 | beta0 = sqrt(one-one/gamma0**2) |
---|
32 | ptot = beta0*gamma0*(one+dp) |
---|
33 | beta_dp = ptot / sqrt(one + ptot**2) |
---|
34 | b_dir_int = node_value('bbdir ') |
---|
35 | b_dir=dble(b_dir_int) |
---|
36 | b_dir = b_dir/sqrt(b_dir*b_dir + 1.0d-32) |
---|
37 | !---- pre-factor, if zero, anything else does not need to be calculated |
---|
38 | fk = two*parvec(5)*parvec(6)/parvec(7)/beta0/(one+dp)/q* & |
---|
39 | (one-beta0*beta_dp*b_dir)/(beta_dp+0.5*(b_dir-one)*b_dir*beta0) |
---|
40 | ! |
---|
41 | if (fk .eq. zero) return |
---|
42 | !---- choose beamshape: 1-Gaussian, 2-flattop=trapezoidal, 3-hollow-parabolic |
---|
43 | beamshape = node_value('bbshape ') |
---|
44 | if(beamshape.lt.1.or.beamshape.gt.3) then |
---|
45 | beamshape=1 |
---|
46 | if(first) then |
---|
47 | first = .false. |
---|
48 | call aawarn('TTBB: ', & |
---|
49 | 'beamshape out of range, set to default=1') |
---|
50 | endif |
---|
51 | endif |
---|
52 | if(beamshape.eq.1) call ttbb_gauss(track,ktrack,fk) |
---|
53 | if(beamshape.eq.2) call ttbb_flattop(track,ktrack,fk) |
---|
54 | if(beamshape.eq.3) call ttbb_hollowparabolic(track,ktrack,fk) |
---|
55 | |
---|
56 | end subroutine ttbb |
---|
57 | subroutine ttbb_gauss(track,ktrack,fk) |
---|
58 | |
---|
59 | use bbfi |
---|
60 | implicit none |
---|
61 | |
---|
62 | ! ---------------------------------------------------------------------* |
---|
63 | ! purpose: kicks the particles of the beam considered with a beam * |
---|
64 | ! having a Gaussian perpendicular shape * |
---|
65 | ! input and output as those of subroutine ttbb * |
---|
66 | ! ---------------------------------------------------------------------* |
---|
67 | logical bborbit |
---|
68 | integer ktrack,itrack,ipos,get_option |
---|
69 | double precision track(6,*),pi,sx,sy,xm,ym,sx2,sy2,xs, & |
---|
70 | ys,rho2,fk,tk,phix,phiy,rk,xb,yb,crx,cry,xr,yr,r,r2,cbx,cby, & |
---|
71 | get_variable,node_value,zero,one,two,three,ten3m,explim |
---|
72 | parameter(zero=0d0,one=1d0,two=2d0,three=3d0,ten3m=1d-3, & |
---|
73 | explim=150d0) |
---|
74 | ! if x > explim, exp(-x) is outside machine limits. |
---|
75 | |
---|
76 | !---- initialize. |
---|
77 | bborbit = get_option('bborbit ') .ne. 0 |
---|
78 | pi=get_variable('pi ') |
---|
79 | sx = node_value('sigx ') |
---|
80 | sy = node_value('sigy ') |
---|
81 | xm = node_value('xma ') |
---|
82 | ym = node_value('yma ') |
---|
83 | ipos = 0 |
---|
84 | if (.not. bborbit) then |
---|
85 | !--- find position of closed orbit bb_kick |
---|
86 | do ipos = 1, bbd_cnt |
---|
87 | if (bbd_loc(ipos) .eq. bbd_pos) goto 1 |
---|
88 | enddo |
---|
89 | ipos = 0 |
---|
90 | 1 continue |
---|
91 | endif |
---|
92 | sx2 = sx*sx |
---|
93 | sy2 = sy*sy |
---|
94 | !---- limit formulae for sigma(x) = sigma(y). |
---|
95 | if (abs(sx2 - sy2) .le. ten3m * (sx2 + sy2)) then |
---|
96 | do itrack = 1, ktrack |
---|
97 | xs = track(1,itrack) - xm |
---|
98 | ys = track(3,itrack) - ym |
---|
99 | rho2 = xs * xs + ys * ys |
---|
100 | tk = rho2 / (two * sx2) |
---|
101 | if (tk .gt. explim) then |
---|
102 | phix = xs * fk / rho2 |
---|
103 | phiy = ys * fk / rho2 |
---|
104 | else if (rho2 .ne. zero) then |
---|
105 | phix = xs * fk / rho2 * (one - exp(-tk) ) |
---|
106 | phiy = ys * fk / rho2 * (one - exp(-tk) ) |
---|
107 | else |
---|
108 | phix = zero |
---|
109 | phiy = zero |
---|
110 | endif |
---|
111 | if (ipos .ne. 0) then |
---|
112 | !--- subtract closed orbit kick |
---|
113 | phix = phix - bb_kick(1,ipos) |
---|
114 | phiy = phiy - bb_kick(2,ipos) |
---|
115 | endif |
---|
116 | track(2,itrack) = track(2,itrack) + phix |
---|
117 | track(4,itrack) = track(4,itrack) + phiy |
---|
118 | enddo |
---|
119 | |
---|
120 | !---- case sigma(x) > sigma(y). |
---|
121 | else if (sx2 .gt. sy2) then |
---|
122 | r2 = two * (sx2 - sy2) |
---|
123 | r = sqrt(r2) |
---|
124 | rk = fk * sqrt(pi) / r |
---|
125 | do itrack = 1, ktrack |
---|
126 | xs = track(1,itrack) - xm |
---|
127 | ys = track(3,itrack) - ym |
---|
128 | xr = abs(xs) / r |
---|
129 | yr = abs(ys) / r |
---|
130 | call ccperrf(xr, yr, crx, cry) |
---|
131 | tk = (xs * xs / sx2 + ys * ys / sy2) / two |
---|
132 | if (tk .gt. explim) then |
---|
133 | phix = rk * cry |
---|
134 | phiy = rk * crx |
---|
135 | else |
---|
136 | xb = (sy / sx) * xr |
---|
137 | yb = (sx / sy) * yr |
---|
138 | call ccperrf(xb, yb, cbx, cby) |
---|
139 | phix = rk * (cry - exp(-tk) * cby) |
---|
140 | phiy = rk * (crx - exp(-tk) * cbx) |
---|
141 | endif |
---|
142 | track(2,itrack) = track(2,itrack) + phix * sign(one,xs) |
---|
143 | track(4,itrack) = track(4,itrack) + phiy * sign(one,ys) |
---|
144 | if (ipos .ne. 0) then |
---|
145 | !--- subtract closed orbit kick |
---|
146 | track(2,itrack) = track(2,itrack) - bb_kick(1,ipos) |
---|
147 | track(4,itrack) = track(4,itrack) - bb_kick(2,ipos) |
---|
148 | endif |
---|
149 | enddo |
---|
150 | |
---|
151 | !---- case sigma(x) < sigma(y). |
---|
152 | else |
---|
153 | r2 = two * (sy2 - sx2) |
---|
154 | r = sqrt(r2) |
---|
155 | rk = fk * sqrt(pi) / r |
---|
156 | do itrack = 1, ktrack |
---|
157 | xs = track(1,itrack) - xm |
---|
158 | ys = track(3,itrack) - ym |
---|
159 | xr = abs(xs) / r |
---|
160 | yr = abs(ys) / r |
---|
161 | call ccperrf(yr, xr, cry, crx) |
---|
162 | tk = (xs * xs / sx2 + ys * ys / sy2) / two |
---|
163 | if (tk .gt. explim) then |
---|
164 | phix = rk * cry |
---|
165 | phiy = rk * crx |
---|
166 | else |
---|
167 | xb = (sy / sx) * xr |
---|
168 | yb = (sx / sy) * yr |
---|
169 | call ccperrf(yb, xb, cby, cbx) |
---|
170 | phix = rk * (cry - exp(-tk) * cby) |
---|
171 | phiy = rk * (crx - exp(-tk) * cbx) |
---|
172 | endif |
---|
173 | track(2,itrack) = track(2,itrack) + phix * sign(one,xs) |
---|
174 | track(4,itrack) = track(4,itrack) + phiy * sign(one,ys) |
---|
175 | if (ipos .ne. 0) then |
---|
176 | !--- subtract closed orbit kick |
---|
177 | track(2,itrack) = track(2,itrack) - bb_kick(1,ipos) |
---|
178 | track(4,itrack) = track(4,itrack) - bb_kick(2,ipos) |
---|
179 | endif |
---|
180 | enddo |
---|
181 | endif |
---|
182 | end subroutine ttbb_gauss |
---|
183 | subroutine ttbb_flattop(track,ktrack,fk) |
---|
184 | |
---|
185 | use bbfi |
---|
186 | implicit none |
---|
187 | |
---|
188 | ! ---------------------------------------------------------------------* |
---|
189 | ! purpose: kicks the particles of the beam considered with a beam * |
---|
190 | ! having an trapezoidal and, so flat top radial profile * |
---|
191 | ! input and output as those of subroutine ttbb * |
---|
192 | ! ---------------------------------------------------------------------* |
---|
193 | logical bborbit,first |
---|
194 | integer ktrack,itrack,ipos,get_option |
---|
195 | double precision track(6,*),pi,r0x,r0y,wi,wx,wy,xm,ym, & |
---|
196 | r0x2,r0y2,xs,ys,rho,rho2,fk,phir,phix,phiy,get_variable, & |
---|
197 | node_value,zero,one,two,three,ten3m,explim,norm,r1,zz |
---|
198 | parameter(zero=0d0,one=1d0,two=2d0,three=3d0,ten3m=1d-3, & |
---|
199 | explim=150d0) |
---|
200 | save first |
---|
201 | data first / .true. / |
---|
202 | ! if x > explim, exp(-x) is outside machine limits. |
---|
203 | |
---|
204 | !---- initialize. |
---|
205 | bborbit = get_option('bborbit ') .ne. 0 |
---|
206 | pi=get_variable('pi ') |
---|
207 | ! mean radii of the is given via variables sigx and sigy |
---|
208 | r0x = node_value('sigx ') |
---|
209 | r0y = node_value('sigy ') |
---|
210 | wi = node_value('width ') |
---|
211 | xm = node_value('xma ') |
---|
212 | ym = node_value('yma ') |
---|
213 | ipos = 0 |
---|
214 | if (.not. bborbit) then |
---|
215 | !--- find position of closed orbit bb_kick |
---|
216 | do ipos = 1, bbd_cnt |
---|
217 | if (bbd_loc(ipos) .eq. bbd_pos) goto 1 |
---|
218 | enddo |
---|
219 | ipos = 0 |
---|
220 | 1 continue |
---|
221 | endif |
---|
222 | r0x2 = r0x*r0x |
---|
223 | r0y2 = r0y*r0y |
---|
224 | wx = r0x*wi |
---|
225 | wy = r0y*wi |
---|
226 | !---- limit formulae for mean radius(x) = mean radius(y), |
---|
227 | !----- preliminary the only case considered. |
---|
228 | ! |
---|
229 | if (abs(r0x2 - r0y2) .gt. ten3m * (r0x2 + r0y2)) then |
---|
230 | zz = 0.5*(r0x + r0y) |
---|
231 | r0x=zz |
---|
232 | r0y=zz |
---|
233 | r0x2=r0x*r0x |
---|
234 | r0y2=r0y*r0y |
---|
235 | if(first) then |
---|
236 | first=.false. |
---|
237 | call aawarn('TTBB_FLATTOP: ','beam is assumed to be circular') |
---|
238 | endif |
---|
239 | endif |
---|
240 | norm = (12.0*r0x**2 + wx**2)/24.0 |
---|
241 | r1=r0x-wx/2.0 |
---|
242 | do itrack = 1, ktrack |
---|
243 | xs = track(1,itrack) - xm |
---|
244 | ys = track(3,itrack) - ym |
---|
245 | rho2 = xs * xs + ys * ys |
---|
246 | rho = sqrt(rho2) |
---|
247 | if(rho.le.r1) then |
---|
248 | phir = 0.5/norm |
---|
249 | phix = phir*xs |
---|
250 | phiy = phir*ys |
---|
251 | else if(rho.gt.r1.and.rho.lt.r1+wx) then |
---|
252 | phir = ((r0x**2/4.0 - r0x**3/6.0/wx - r0x*wx/8.0 + & |
---|
253 | wx**2/48.0)/rho2 + 0.25 + 0.5*r0x/wx - & |
---|
254 | rho/3.0/wx)/norm |
---|
255 | phix = phir*xs |
---|
256 | phiy = phir*ys |
---|
257 | else if(rho.ge.r1+wx) then |
---|
258 | phir = 1.0/rho2 |
---|
259 | phix = xs*phir |
---|
260 | phiy = ys*phir |
---|
261 | endif |
---|
262 | track(2,itrack) = track(2,itrack)+phix*fk |
---|
263 | track(4,itrack) = track(4,itrack)+phiy*fk |
---|
264 | end do |
---|
265 | |
---|
266 | end subroutine ttbb_flattop |
---|
267 | subroutine ttbb_hollowparabolic(track,ktrack,fk) |
---|
268 | |
---|
269 | use bbfi |
---|
270 | implicit none |
---|
271 | |
---|
272 | ! ---------------------------------------------------------------------* |
---|
273 | ! purpose: kicks the particles of the beam considered with a beam * |
---|
274 | ! having a hollow-parabolic perpendicular shape * |
---|
275 | ! input and output as those of subroutine ttbb * |
---|
276 | ! ---------------------------------------------------------------------* |
---|
277 | logical bborbit,first |
---|
278 | integer ktrack,itrack,ipos,get_option |
---|
279 | double precision track(6,*),pi,r0x,r0y,wi,wx,wy,xm,ym, & |
---|
280 | r0x2,r0y2,xs,ys,rho,rho2,fk,phir,phix,phiy,get_variable, & |
---|
281 | node_value,zero,one,two,three,ten3m,explim,zz |
---|
282 | parameter(zero=0d0,one=1d0,two=2d0,three=3d0,ten3m=1d-3, & |
---|
283 | explim=150d0) |
---|
284 | save first |
---|
285 | data first / .true. / |
---|
286 | ! if x > explim, exp(-x) is outside machine limits. |
---|
287 | |
---|
288 | !---- initialize. |
---|
289 | bborbit = get_option('bborbit ') .ne. 0 |
---|
290 | pi=get_variable('pi ') |
---|
291 | ! mean radii of the is given via variables sigx and sigy |
---|
292 | r0x = node_value('sigx ') |
---|
293 | r0y = node_value('sigy ') |
---|
294 | wi = node_value('width ') |
---|
295 | ! width is given as FWHM of parabolic density profile, but formulas were |
---|
296 | ! derived with half width at the bottom of the parabolic density profile |
---|
297 | wi = wi/sqrt(2.0) |
---|
298 | xm = node_value('xma ') |
---|
299 | ym = node_value('yma ') |
---|
300 | ipos = 0 |
---|
301 | if (.not. bborbit) then |
---|
302 | !--- find position of closed orbit bb_kick |
---|
303 | do ipos = 1, bbd_cnt |
---|
304 | if (bbd_loc(ipos) .eq. bbd_pos) goto 1 |
---|
305 | enddo |
---|
306 | ipos = 0 |
---|
307 | 1 continue |
---|
308 | endif |
---|
309 | r0x2 = r0x*r0x |
---|
310 | r0y2 = r0y*r0y |
---|
311 | wx = wi*r0x |
---|
312 | wy = wi*r0y |
---|
313 | !---- limit formulae for mean radius(x) = mean radius(y), |
---|
314 | !----- preliminary the only case considered. |
---|
315 | ! |
---|
316 | if (abs(r0x2 - r0y2) .gt. ten3m * (r0x2 + r0y2)) then |
---|
317 | zz = 0.5*(r0x + r0y) |
---|
318 | r0x=zz |
---|
319 | r0y=zz |
---|
320 | r0x2=r0x*r0x |
---|
321 | r0y2=r0y*r0y |
---|
322 | if(first) then |
---|
323 | first=.false. |
---|
324 | call aawarn('TTBB_HOLLOWPARABOLIC: ', & |
---|
325 | 'beam is assumed to be circular') |
---|
326 | endif |
---|
327 | endif |
---|
328 | do itrack = 1, ktrack |
---|
329 | xs = track(1,itrack) - xm |
---|
330 | ys = track(3,itrack) - ym |
---|
331 | rho2 = xs * xs + ys * ys |
---|
332 | rho = sqrt(rho2) |
---|
333 | if(rho.le.r0x-wx) then |
---|
334 | phix = zero |
---|
335 | phiy = zero |
---|
336 | else if(rho.gt.r0x-wx.and.rho.lt.r0x+wx) then |
---|
337 | phir=0.75/wx/r0x/rho2*(r0x**4/12.0/wx**2 - r0x**2/2.0 + & |
---|
338 | 2.0*r0x*wx/3.0 - wx**2/4.0 + rho2/2.0*(1.0 - & |
---|
339 | r0x**2/wx**2) + rho**3/3.0*2.0*r0x/wx**2 - & |
---|
340 | rho**4/4.0/wx**2) |
---|
341 | phix = phir*xs |
---|
342 | phiy = phir*ys |
---|
343 | else |
---|
344 | phir = 1.0/rho2 |
---|
345 | phix = xs*phir |
---|
346 | phiy = ys*phir |
---|
347 | endif |
---|
348 | track(2,itrack) = track(2,itrack)+phix*fk |
---|
349 | track(4,itrack) = track(4,itrack)+phiy*fk |
---|
350 | end do |
---|
351 | |
---|
352 | end subroutine ttbb_hollowparabolic |
---|
353 | !!$subroutine ttbb_old(track,ktrack,parvec) |
---|
354 | !!$ |
---|
355 | !!$ use bbfi |
---|
356 | !!$ implicit none |
---|
357 | !!$ |
---|
358 | !!$ !----------------------------------------------------------------------* |
---|
359 | !!$ ! purpose: * |
---|
360 | !!$ ! track a set of particle through a beam-beam interaction region. * |
---|
361 | !!$ ! see mad physicist's manual for the formulas used. * |
---|
362 | !!$ !input: * |
---|
363 | !!$ ! input/output: * |
---|
364 | !!$ ! track(6,*)(double) track coordinates: (x, px, y, py, t, pt). * |
---|
365 | !!$ ! ktrack (integer) number of tracks. * |
---|
366 | !!$ !----------------------------------------------------------------------* |
---|
367 | !!$ logical bborbit |
---|
368 | !!$ integer ktrack,itrack,ipos,get_option |
---|
369 | !!$ double precision track(6,*),parvec(*),pi,sx,sy,xm,ym,sx2,sy2,xs, & |
---|
370 | !!$ ys,rho2,fk,tk,phix,phiy,rk,xb,yb,crx,cry,xr,yr,r,r2,cbx,cby, & |
---|
371 | !!$ get_variable,node_value,zero,one,two,three,ten3m,explim |
---|
372 | !!$ parameter(zero=0d0,one=1d0,two=2d0,three=3d0,ten3m=1d-3, & |
---|
373 | !!$ explim=150d0) |
---|
374 | !!$ ! if x > explim, exp(-x) is outside machine limits. |
---|
375 | !!$ |
---|
376 | !!$ !---- initialize. |
---|
377 | !!$ bborbit = get_option('bborbit ') .ne. 0 |
---|
378 | !!$ pi=get_variable('pi ') |
---|
379 | !!$ sx = node_value('sigx ') |
---|
380 | !!$ sy = node_value('sigy ') |
---|
381 | !!$ xm = node_value('xma ') |
---|
382 | !!$ ym = node_value('yma ') |
---|
383 | !!$ fk = two * parvec(5) * parvec(6) / parvec(7) |
---|
384 | !!$ if (fk .eq. zero) return |
---|
385 | !!$ ipos = 0 |
---|
386 | !!$ if (.not. bborbit) then |
---|
387 | !!$ !--- find position of closed orbit bb_kick |
---|
388 | !!$ do ipos = 1, bbd_cnt |
---|
389 | !!$ if (bbd_loc(ipos) .eq. bbd_pos) goto 1 |
---|
390 | !!$ enddo |
---|
391 | !!$ ipos = 0 |
---|
392 | !!$1 continue |
---|
393 | !!$ endif |
---|
394 | !!$ sx2 = sx*sx |
---|
395 | !!$ sy2 = sy*sy |
---|
396 | !!$ !---- limit formulae for sigma(x) = sigma(y). |
---|
397 | !!$ if (abs(sx2 - sy2) .le. ten3m * (sx2 + sy2)) then |
---|
398 | !!$ do itrack = 1, ktrack |
---|
399 | !!$ xs = track(1,itrack) - xm |
---|
400 | !!$ ys = track(3,itrack) - ym |
---|
401 | !!$ rho2 = xs * xs + ys * ys |
---|
402 | !!$ tk = rho2 / (two * sx2) |
---|
403 | !!$ if (tk .gt. explim) then |
---|
404 | !!$ phix = xs * fk / rho2 |
---|
405 | !!$ phiy = ys * fk / rho2 |
---|
406 | !!$ else if (rho2 .ne. zero) then |
---|
407 | !!$ phix = xs * fk / rho2 * (one - exp(-tk) ) |
---|
408 | !!$ phiy = ys * fk / rho2 * (one - exp(-tk) ) |
---|
409 | !!$ else |
---|
410 | !!$ phix = zero |
---|
411 | !!$ phiy = zero |
---|
412 | !!$ endif |
---|
413 | !!$ if (ipos .ne. 0) then |
---|
414 | !!$ !--- subtract closed orbit kick |
---|
415 | !!$ phix = phix - bb_kick(1,ipos) |
---|
416 | !!$ phiy = phiy - bb_kick(2,ipos) |
---|
417 | !!$ endif |
---|
418 | !!$ track(2,itrack) = track(2,itrack) + phix |
---|
419 | !!$ track(4,itrack) = track(4,itrack) + phiy |
---|
420 | !!$ enddo |
---|
421 | !!$ |
---|
422 | !!$ !---- case sigma(x) > sigma(y). |
---|
423 | !!$ else if (sx2 .gt. sy2) then |
---|
424 | !!$ r2 = two * (sx2 - sy2) |
---|
425 | !!$ r = sqrt(r2) |
---|
426 | !!$ rk = fk * sqrt(pi) / r |
---|
427 | !!$ do itrack = 1, ktrack |
---|
428 | !!$ xs = track(1,itrack) - xm |
---|
429 | !!$ ys = track(3,itrack) - ym |
---|
430 | !!$ xr = abs(xs) / r |
---|
431 | !!$ yr = abs(ys) / r |
---|
432 | !!$ call ccperrf(xr, yr, crx, cry) |
---|
433 | !!$ tk = (xs * xs / sx2 + ys * ys / sy2) / two |
---|
434 | !!$ if (tk .gt. explim) then |
---|
435 | !!$ phix = rk * cry |
---|
436 | !!$ phiy = rk * crx |
---|
437 | !!$ else |
---|
438 | !!$ xb = (sy / sx) * xr |
---|
439 | !!$ yb = (sx / sy) * yr |
---|
440 | !!$ call ccperrf(xb, yb, cbx, cby) |
---|
441 | !!$ phix = rk * (cry - exp(-tk) * cby) |
---|
442 | !!$ phiy = rk * (crx - exp(-tk) * cbx) |
---|
443 | !!$ endif |
---|
444 | !!$ track(2,itrack) = track(2,itrack) + phix * sign(one,xs) |
---|
445 | !!$ track(4,itrack) = track(4,itrack) + phiy * sign(one,ys) |
---|
446 | !!$ if (ipos .ne. 0) then |
---|
447 | !!$ !--- subtract closed orbit kick |
---|
448 | !!$ track(2,itrack) = track(2,itrack) - bb_kick(1,ipos) |
---|
449 | !!$ track(4,itrack) = track(4,itrack) - bb_kick(2,ipos) |
---|
450 | !!$ endif |
---|
451 | !!$ enddo |
---|
452 | !!$ |
---|
453 | !!$ !---- case sigma(x) < sigma(y). |
---|
454 | !!$ else |
---|
455 | !!$ r2 = two * (sy2 - sx2) |
---|
456 | !!$ r = sqrt(r2) |
---|
457 | !!$ rk = fk * sqrt(pi) / r |
---|
458 | !!$ do itrack = 1, ktrack |
---|
459 | !!$ xs = track(1,itrack) - xm |
---|
460 | !!$ ys = track(3,itrack) - ym |
---|
461 | !!$ xr = abs(xs) / r |
---|
462 | !!$ yr = abs(ys) / r |
---|
463 | !!$ call ccperrf(yr, xr, cry, crx) |
---|
464 | !!$ tk = (xs * xs / sx2 + ys * ys / sy2) / two |
---|
465 | !!$ if (tk .gt. explim) then |
---|
466 | !!$ phix = rk * cry |
---|
467 | !!$ phiy = rk * crx |
---|
468 | !!$ else |
---|
469 | !!$ xb = (sy / sx) * xr |
---|
470 | !!$ yb = (sx / sy) * yr |
---|
471 | !!$ call ccperrf(yb, xb, cby, cbx) |
---|
472 | !!$ phix = rk * (cry - exp(-tk) * cby) |
---|
473 | !!$ phiy = rk * (crx - exp(-tk) * cbx) |
---|
474 | !!$ endif |
---|
475 | !!$ track(2,itrack) = track(2,itrack) + phix * sign(one,xs) |
---|
476 | !!$ track(4,itrack) = track(4,itrack) + phiy * sign(one,ys) |
---|
477 | !!$ if (ipos .ne. 0) then |
---|
478 | !!$ !--- subtract closed orbit kick |
---|
479 | !!$ track(2,itrack) = track(2,itrack) - bb_kick(1,ipos) |
---|
480 | !!$ track(4,itrack) = track(4,itrack) - bb_kick(2,ipos) |
---|
481 | !!$ endif |
---|
482 | !!$ enddo |
---|
483 | !!$ endif |
---|
484 | !!$end subroutine ttbb_old |
---|