source: PSPA/madxPSPA/libs/ptc/src/Sma0_beam_beam_ptc.f90 @ 430

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

import madx-5.01.00
File size: 15.0 KB
Line 
1module beam_beam_ptc
2  !  use orbit_ptc
3  use S_TRACKING
4  ! use madx_ptc_module
5  implicit none
6  public
7  private BBKICKP, BBKICKR,PATCH_BBR,PATCH_BBP
8  private ccperrfP, ccperrfr ,ccperrf
9  !  private TRACK_NODE_LAYOUT_FLAG_R,TRACK_NODE_LAYOUT_FLAG_P
10  private imax
11  integer ::imax=1000
12
13  INTERFACE ccperrf
14     MODULE PROCEDURE ccperrfP
15     MODULE PROCEDURE ccperrfr
16  END INTERFACE
17
18
19  INTERFACE BBKICK
20     MODULE PROCEDURE BBKICKR
21     MODULE PROCEDURE BBKICKP
22  END INTERFACE
23  INTERFACE PATCH_BB
24     MODULE PROCEDURE PATCH_BBR
25     MODULE PROCEDURE PATCH_BBP
26  END INTERFACE
27
28  logical(lp), target :: do_beam_beam= my_false
29
30contains
31
32  SUBROUTINE PATCH_BBR(B,X,k,BETA0,exact,ENTERING)
33    implicit none
34    ! MISALIGNS REAL FIBRES IN PTC ORDER FOR FORWARD AND BACKWARD FIBRES
35    TYPE(BEAM_BEAM_NODE),TARGET,INTENT(INOUT):: B
36    real(dp), INTENT(INOUT):: X(6)
37    logical(lp),INTENT(IN):: exact,ENTERING
38    REAL(DP),INTENT(IN):: BETA0
39    REAL(DP) A(3),D(3)
40    TYPE(INTERNAL_STATE) k !,OPTIONAL :: K
41
42    IF(ENTERING) THEN
43       X(3)=B%A_X1*X(3);X(4)=B%A_X1*X(4);
44       CALL ROT_YZ(B%A(1),X,BETA0,exact,k%TIME)
45       CALL ROT_XZ(B%A(2),X,BETA0,exact,k%TIME)
46       CALL ROT_XY(B%A(3),X)  !,exact)
47       CALL TRANS(B%D,X,BETA0,exact,k%TIME)
48       X(3)=B%A_X2*X(3);X(4)=B%A_X2*X(4);
49    ELSE
50       A=-B%A
51       D=-B%D
52       X(3)=B%A_X2*X(3);X(4)=B%A_X2*X(4);
53       CALL TRANS(D,X,BETA0,exact,k%TIME)
54       CALL ROT_XY(A(3),X)  !,exact)
55       CALL ROT_XZ(A(2),X,BETA0,exact,k%TIME)
56       CALL ROT_YZ(A(1),X,BETA0,exact,k%TIME)
57       X(3)=B%A_X1*X(3);X(4)=B%A_X1*X(4);
58    ENDIF
59
60
61  END SUBROUTINE PATCH_BBR
62
63  SUBROUTINE PATCH_BBP(B,X,k,BETA0,exact,ENTERING)
64    implicit none
65    ! MISALIGNS REAL FIBRES IN PTC ORDER FOR FORWARD AND BACKWARD FIBRES
66    TYPE(BEAM_BEAM_NODE),TARGET,INTENT(INOUT):: B
67    TYPE(REAL_8), INTENT(INOUT):: X(6)
68    logical(lp),INTENT(IN):: exact,ENTERING
69    REAL(DP),INTENT(IN):: BETA0
70    REAL(DP) A(3),D(3)
71    TYPE(INTERNAL_STATE) k !,OPTIONAL :: K
72
73    IF(ENTERING) THEN
74       X(3)=B%A_X1*X(3);X(4)=B%A_X1*X(4);
75       CALL ROT_YZ(B%A(1),X,BETA0,exact,k%TIME)
76       CALL ROT_XZ(B%A(2),X,BETA0,exact,k%TIME)
77       CALL ROT_XY(B%A(3),X)  !,exact)
78       CALL TRANS(B%D,X,BETA0,exact,k%TIME)
79       X(3)=B%A_X2*X(3);X(4)=B%A_X2*X(4);
80    ELSE
81       A=-B%A
82       D=-B%D
83       X(3)=B%A_X2*X(3);X(4)=B%A_X2*X(4);
84       CALL TRANS(D,X,BETA0,exact,k%TIME)
85       CALL ROT_XY(A(3),X)  !,exact)
86       CALL ROT_XZ(A(2),X,BETA0,exact,k%TIME)
87       CALL ROT_YZ(A(1),X,BETA0,exact,k%TIME)
88       X(3)=B%A_X1*X(3);X(4)=B%A_X1*X(4);
89    ENDIF
90
91
92  END SUBROUTINE PATCH_BBP
93
94
95  subroutine BBKICKR(BB,X)
96
97    implicit none
98    !----------------------------------------------------------------------*
99    ! purpose:                                                             *
100    !   track a set of particle through a beam-beam interaction region.    *
101    !   see mad physicist's manual for the formulas used.                  *
102    !input:                                                                *
103    ! input/output:                                                        *
104    !   b%x(:,6)(double)  track coordinates: (x, px, y, py,  pt,t).        *
105    !   b%n    (integer) number of tracks.                                 *
106    !----------------------------------------------------------------------*
107    real(dp) sx2,sy2,xs,ys,rho2,fk,tk,phix,phiy,rk,xb,yb,crx,cry,xr,yr,r,r2,&
108         cbx,cby,ten3m,explim,sx,sy,xm,ym
109    TYPE(BEAM_BEAM_NODE), INTENT(INOUT) ::BB
110    REAL(DP), INTENT(INOUT) :: X(6)
111    parameter(ten3m=1.0e-3_dp,explim=150.0_dp)
112    SX=BB%SX
113    SY=BB%SY
114    XM=BB%XM
115    YM=BB%YM
116    FK=BB%FK
117    !    write(6,*) "bb%FK = " ,bb%FK
118
119    if (fk == 0.0_dp)  return
120    sx2 = sx*sx
121    sy2 = sy*sy
122    !---- limit formulae for sigma(x) = sigma(y).
123    if (abs(sx2 - sy2) .le. ten3m * (sx2 + sy2)) then !ten3m = 1.0d-3
124       xs = x(1) - xm
125       ys = x(3) - ym
126       rho2 = xs * xs + ys * ys
127       tk = rho2 / (sx2 + sy2)
128       if (tk .gt. explim) then
129          phix = xs * fk / rho2
130          phiy = ys * fk / rho2
131       else if (rho2 .ne. 0.0_dp) then
132          phix = xs * fk / rho2 * (1.0_dp - exp(-tk) )
133          phiy = ys * fk / rho2 * (1.0_dp - exp(-tk) )
134       else
135          phix = 0.0_dp
136          phiy = 0.0_dp
137       endif
138       phix = phix - bb%bbk(1) ! subtract horizontal bb kick
139       phiy = phiy - bb%bbk(2) ! subtract vertical co
140       x(2) = x(2) + phix
141       x(4) = x(4) + phiy
142
143       !---- case sigma(x) > sigma(y).
144    else if (sx2 > sy2) then
145       r2 = 2.0_dp * (sx2 - sy2)
146       r  = sqrt(r2)
147       rk = fk * sqrt(pi) / r
148       xs = x(1) - xm
149       ys = x(3) - ym
150       xr = abs(xs) / r
151       yr = abs(ys) / r
152       call ccperrf(xr, yr, crx, cry)
153       tk = (xs * xs / sx2 + ys * ys / sy2) / 2.0_dp
154       if (tk .gt. explim) then
155          phix = rk * cry
156          phiy = rk * crx
157       else
158          xb = (sy / sx) * xr
159          yb = (sx / sy) * yr
160          call ccperrf(xb, yb, cbx, cby)
161          phix = rk * (cry - exp(-tk) * cby)
162          phiy = rk * (crx - exp(-tk) * cbx)
163       endif
164       x(2) = x(2) + phix * sign(1.0_dp,xs)
165       x(4) = x(4) + phiy * sign(1.0_dp,ys)
166       x(2) = x(2) - BB%bbk(1)
167       x(4) = x(4) - BB%bbk(2)
168
169       !---- case sigma(x) < sigma(y).
170    else
171       r2 = 2.0_dp * (sy2 - sx2)
172       r  = sqrt(r2)
173       rk = fk * sqrt(pi) / r
174       !       do itrack = 1, b%n
175       !         IF(B%U(itrack)) CYCLE
176       !          xs = track(1,itrack) - xm
177       !          ys = track(3,itrack) - ym
178       xs = x(1) - xm
179       ys = x(3) - ym
180       xr = abs(xs) / r
181       yr = abs(ys) / r
182       call ccperrf(yr, xr, cry, crx)
183       tk = (xs * xs / sx2 + ys * ys / sy2) / 2.0_dp
184       if (tk .gt. explim) then
185          phix = rk * cry
186          phiy = rk * crx
187       else
188          xb  = (sy / sx) * xr
189          yb  = (sx / sy) * yr
190          call ccperrf(yb, xb, cby, cbx)
191          phix = rk * (cry - exp(-tk) * cby)
192          phiy = rk * (crx - exp(-tk) * cbx)
193       endif
194
195       x(2) = x(2) + phix * sign(1.0_dp,xs)
196       x(4) = x(4) + phiy * sign(1.0_dp,ys)
197
198       x(2) = x(2) - BB%bbk(1)
199       x(4) = x(4) - BB%bbk(2)
200
201    endif
202    !    IF(DZ/=ZERO) CALL DRIFT_BEAM_BACK_TO_POSITION(th,DZ,B)
203
204  end subroutine BBKICKR
205
206  subroutine ccperrfr(xx, yy, wx, wy)
207    implicit none
208    !----------------------------------------------------------------------*
209    ! purpose:                                                             *
210    !   modification of wwerf, double precision complex error function,    *
211    !   written at cern by K. Koelbig.                                     *
212    ! input:                                                               *
213    !   xx, yy    (double)    real + imag argument                         *
214    ! output:                                                              *
215    !   wx, wy    (double)    real + imag function result                  *
216    !----------------------------------------------------------------------*
217    integer n,nc,nu
218    real(dp), INTENT(INOUT):: xx,yy,wx,wy
219    real(dp) x,y,q,h,xl,xh,yh,tx,ty,tn,sx,sy,saux,  &
220         rx(33),ry(33),cc,xlim,ylim,fac1,fac2,fac3
221    parameter(cc=1.12837916709551_dp,        &
222         xlim=5.33_dp,ylim=4.29_dp,fac1=3.2_dp,fac2=23.0_dp,fac3=21.0_dp)
223
224    x = abs(xx)
225    y = abs(yy)
226
227    if (y .lt. ylim  .and.  x .lt. xlim) then
228       q  = (1.0_dp - y / ylim) * sqrt(1.0_dp - (x/xlim)**2)
229       h  = 1.0_dp / (fac1 * q)
230       nc = 7 + int(fac2*q)
231       xl = h**(1 - nc)
232       xh = y + 0.5_dp/h
233       yh = x
234       nu = 10 + int(fac3*q)
235       rx(nu+1) = 0.0_dp
236       ry(nu+1) = 0.0_dp
237
238       do n = nu, 1, -1
239          tx = xh + n * rx(n+1)
240          ty = yh - n * ry(n+1)
241          tn = tx*tx + ty*ty
242          rx(n) = 0.5_dp * tx / tn
243          ry(n) = 0.5_dp * ty / tn
244       enddo
245
246       sx = 0.0_dp
247       sy = 0.0_dp
248
249       do n = nc, 1, -1
250          saux = sx + xl
251          sx = rx(n) * saux - ry(n) * sy
252          sy = rx(n) * sy + ry(n) * saux
253          xl = h * xl
254       enddo
255
256       wx = cc * sx
257       wy = cc * sy
258    else
259       xh = y
260       yh = x
261       rx(1) = 0.0_dp
262       ry(1) = 0.0_dp
263
264       do n = 9, 1, -1
265          tx = xh + n * rx(1)
266          ty = yh - n * ry(1)
267          tn = tx*tx + ty*ty
268          rx(1) = 0.5_dp * tx / tn
269          ry(1) = 0.5_dp * ty / tn
270       enddo
271
272       wx = cc * rx(1)
273       wy = cc * ry(1)
274    endif
275
276    !      if(y .eq. zero) wx = exp(-x**2)
277    if(yy .lt. 0.0_dp) then
278       wx =   2.0_dp * exp(y*y-x*x) * cos(2.0_dp*x*y) - wx
279       wy = - 2.0_dp * exp(y*y-x*x) * sin(2.0_dp*x*y) - wy
280       if(xx .gt. 0.0_dp) wy = -wy
281    else
282       if(xx .lt. 0.0_dp) wy = -wy
283    endif
284
285  end SUBROUTINE ccperrfr
286
287
288
289  subroutine BBKICKP(BB,X)
290
291    implicit none
292    !----------------------------------------------------------------------*
293    ! purpose:                                                             *
294    !   track a set of particle through a beam-beam interaction region.    *
295    !   see mad physicist's manual for the formulas used.                  *
296    !input:                                                                *
297    ! input/output:                                                        *
298    !   b%x(:,6)(double)  track coordinates: (x, px, y, py,  pt,t).      *
299    !   b%n    (integer) number of tracks.                              *
300    !----------------------------------------------------------------------*
301    integer it
302    TYPE(REAL_8) xs,ys,tk,phix,phiy,xb,yb,crx,cry
303    TYPE(REAL_8) xr,yr,cbx,cby,rho2
304    REAL(DP) sx2,sy2,sx,sy,xm,ym,fk,ten3m,explim,xn1,xn2,xs1,xs2,arglim,rk
305    REAL(DP) r,r2,n
306    TYPE(BEAM_BEAM_NODE), INTENT(INOUT) ::BB
307    TYPE(REAL_8), INTENT(INOUT) :: X(6)
308    parameter(ten3m=1.0e-3_dp,arglim=1.0e-2_dp,explim=150.0_dp)
309
310    if (BB%fk == 0.0_dp)  return
311
312    CALL ALLOC(xr,yr,cbx,cby,rho2)
313    CALL ALLOC(xs,ys,tk,phix,phiy)
314    CALL ALLOC(xb,yb,crx,cry)
315
316
317    SX=BB%SX
318    SY=BB%SY
319    XM=BB%XM
320    YM=BB%YM
321    FK=BB%FK
322
323    sx2 = sx*sx
324    sy2 = sy*sy
325    !---- limit formulae for sigma(x) = sigma(y).
326    xn1=abs(sx2 - sy2)
327    xn2= ten3m * (sx2 + sy2)
328    xs1=sx2
329    xs2= sy2
330    if (xn1 .le.xn2 ) then !ten3m = 1.0d-3
331       xs = x(1) - xm
332       ys = x(3) - ym
333       rho2 = xs * xs + ys * ys
334       tk = rho2 / (sx2 + sy2)
335       if (tk .gt. explim) then
336          phix = xs * fk / rho2
337          phiy = ys * fk / rho2
338       else if (tk > arglim) then
339          phix = xs * fk / rho2 * (1.0_dp - exp(-tk) )
340          phiy = ys * fk / rho2 * (1.0_dp - exp(-tk) )
341       else
342
343          xr=1.0_dp
344          yr=1.0_dp
345
346          n=mybig
347          do it=1,imax
348             xr=-xr*tk/(it+1)
349             yr=yr+xr
350             if(it>10)n=full_abs(xr)
351             if(n<=puny) exit
352          enddo
353          if(it>imax-2) then
354             write(6,*) it,n
355             write(6,*) " Stopped in Beam-Beam "
356          endif
357          phix = xs * fk / (2.0_dp * sx2) * YR ! fudge
358          phiY = Ys * fk / (2.0_dp * sx2) * YR ! fudge
359       endif
360
361       phix = phix - bb%bbk(1)
362       phiy = phiy - bb%bbk(2)
363
364       x(2) = x(2) + phix
365       x(4) = x(4) + phiy
366
367
368       !---- case sigma(x) > sigma(y).
369    else
370
371
372       xs = x(1) - xm
373       ys = x(3) - ym
374       tk = (xs * xs / sx2 + ys * ys / sy2) / 2.0_dp
375
376       if(xs1 > xs2) then
377          r2 = 2.0_dp * (sx2 - sy2)
378          r  = sqrt(r2)
379          rk = fk * sqrt(pi) / r
380          xr = xs / r   !
381          yr = ys / r   !
382          call ccperrf(xr, yr, crx, cry)
383          if (tk .gt. explim) then
384             phix = rk * cry
385             phiy = rk * crx
386          else
387             xb = (sy / sx) * xr
388             yb = (sx / sy) * yr
389             call ccperrf(xb, yb, cbx, cby)
390             phix = rk * (cry - exp(-tk) * cby)
391             phiy = rk * (crx - exp(-tk) * cbx)
392          endif
393          x(2) = x(2) + phix
394          x(4) = x(4) + phiy
395          !          if (.NOT.bborbit)  then
396          x(2) = x(2) - BB%bbk(1)
397          x(4) = x(4) - BB%bbk(2)
398          !          endif
399          !       enddo
400
401          !---- case sigma(x) < sigma(y).
402       else
403          r2 = 2.0_dp * (sy2 - sx2)
404          r  = sqrt(r2)
405          rk = fk * sqrt(pi) / r
406
407          !       xs = x(1) - xm
408          !       ys = x(3) - ym
409          xr = xs / r !abs
410          yr = ys / r !abs
411          call ccperrf(yr, xr, cry, crx)
412          !       tk = (xs * xs / sx2 + ys * ys / sy2) / two
413          if (tk .gt. explim) then
414             phix = rk * cry
415             phiy = rk * crx
416          else
417             xb  = (sy / sx) * xr
418             yb  = (sx / sy) * yr
419             call ccperrf(yb, xb, cby, cbx)
420             phix = rk * (cry - exp(-tk) * cby)
421             phiy = rk * (crx - exp(-tk) * cbx)
422          endif
423          !          track(2,itrack) = track(2,itrack) + phix * sign(one,xs)
424          !          track(4,itrack) = track(4,itrack) + phiy * sign(one,ys)
425          x(2) = x(2) + phix
426          x(4) = x(4) + phiy
427          !          if (.NOT.bborbit)  then
428          !--- subtract closed orbit kick
429          !            track(2,itrack) = track(2,itrack) - bb_kick(1,ipos)
430          !            track(4,itrack) = track(4,itrack) - bb_kick(2,ipos)
431          x(2) = x(2) - BB%bbk(1)
432          x(4) = x(4) - BB%bbk(2)
433          !          endif
434          !       enddo
435       endif
436    endif
437
438    CALL KILL(xr,yr,cbx,cby,rho2)
439    CALL KILL(xs,ys,tk,phix,phiy)
440    CALL KILL(xb,yb,crx,cry)
441
442  end subroutine BBKICKP
443
444  subroutine ccperrfP(xx, yy, wx, wy)
445    implicit none
446    !----------------------------------------------------------------------*
447    ! purpose:                                                             *
448    !   modification of wwerf, double precision complex error function,    *
449    !   written at cern by K. Koelbig.                                     *
450    ! input:                                                               *
451    !   xx, yy    (double)    real + imag argument                         *
452    ! output:                                                              *
453    !   wx, wy    (double)    real + imag function result                  *
454    !----------------------------------------------------------------------*
455    TYPE(REAL_8),INTENT(INOUT):: xx,yy,wx,wy
456    TYPE(double_complex) z,zt,w
457    complex(dp) z0,w0,w1,wt0
458    real(dp) xx0, yy0, wx0, wy0
459    integer i
460    call alloc( z,zt,w)
461
462    z=xx+i_*yy
463    z0=z
464    z=z-z0
465
466    xx0=real(z0)
467    yy0=aimag(z0)
468    call ccperrf(xx0, yy0, wx0, wy0)
469
470    w0=wx0+i_*wy0
471
472    w1=-2.0_dp*z0*w0+2.0_dp*i_/sqrt(pi)
473
474    w=w0+w1*z
475
476    zt=z
477
478    do i=2,c_%no
479
480       zt=z*zt
481       wt0=  -2.0_dp*(w0+z0*w1)/i
482
483       w=w+wt0*zt
484
485       w0=w1
486       w1=wt0
487
488    enddo
489
490    wx=real(w)
491    wy=aimag(w)
492
493
494    call kill( z,zt,w)
495
496  end SUBROUTINE ccperrfP
497end module  beam_beam_ptc
Note: See TracBrowser for help on using the repository browser.