[430] | 1 | module 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 | |
---|
| 30 | contains |
---|
| 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 |
---|
| 497 | end module beam_beam_ptc |
---|