| 1 | subroutine emit(deltap, tol, orbit0, disp0, rt, u0, emit_v, & |
|---|
| 2 | nemit_v, bmax, gmax, dismax, tunes, sig_v, pdamp) |
|---|
| 3 | use bbfi |
|---|
| 4 | use twiss0fi |
|---|
| 5 | use emitfi |
|---|
| 6 | implicit none |
|---|
| 7 | |
|---|
| 8 | |
|---|
| 9 | !----------------------------------------------------------------------* |
|---|
| 10 | ! Purpose: * |
|---|
| 11 | ! Compute emittances by A. Chao's method. * |
|---|
| 12 | ! Input: * |
|---|
| 13 | ! DELTAP (real) Average energy error desired. * |
|---|
| 14 | ! tol (real) tolerance for distinction static/dynamic: |
|---|
| 15 | ! |e_val_5| < tol & |e-val_6| < tol: static, |
|---|
| 16 | ! default: 1.000001 |
|---|
| 17 | ! orbit0 (real) closed orbit from TMCLOR |
|---|
| 18 | ! disp0 (real) dispersion from TMCLOR |
|---|
| 19 | ! rt (real) one_turn matrix TMCLOR |
|---|
| 20 | ! Output: |
|---|
| 21 | ! u0 (real) Radiation loss per turn in GeV |
|---|
| 22 | ! emit_v (real) ex, ey, et (emittances) |
|---|
| 23 | ! nemit_v (real) exn, eyn, etn (normalised emitt., MAD style) |
|---|
| 24 | ! bmax (real) Maximum extents of modes. |
|---|
| 25 | ! gmax (real) Maximum divergences of modes. |
|---|
| 26 | ! dismax (real) Maximum dispersion. |
|---|
| 27 | ! tunes (real) qx, qy, qs |
|---|
| 28 | ! pdamp (real) damping partition numbers |
|---|
| 29 | ! sig_v (real) sigx, sigy, sigt, sige |
|---|
| 30 | !----------------------------------------------------------------------* |
|---|
| 31 | !---- Communication area for radiation damping. |
|---|
| 32 | double precision orbit0(6), orbit(6), orbit2(6), em(6,6), rd(6,6), reval(6) |
|---|
| 33 | double precision aival(6), rt(6,6), orbit1(6), ek(6) |
|---|
| 34 | double precision emit_v(3), nemit_v(3), tunes(3), sig_v(4) |
|---|
| 35 | double precision u0, pdamp(3) |
|---|
| 36 | double precision bmax(3,3), gmax(3,3), dismax(4) |
|---|
| 37 | double precision em2(6,6), disp(6), disp0(6),al_errors(align_max) |
|---|
| 38 | double precision tol, deltap, get_value, arad, suml, gammas, el |
|---|
| 39 | double precision betas, bx, gx, re(6,6), te(6,6,6), tt(6,6,6) |
|---|
| 40 | double precision get_variable, node_value |
|---|
| 41 | integer i, j, j1, j2, k, k1, k2, eflag, restart_sequ, n_align |
|---|
| 42 | integer node_al_errors, code, advance_node |
|---|
| 43 | logical m66sta, fmap, stabx, staby, stabt, frad |
|---|
| 44 | double precision zero, one, three, twopi |
|---|
| 45 | parameter (zero = 0.0d0, one = 1.d0, three = 3.0d0) |
|---|
| 46 | |
|---|
| 47 | twopi = get_variable('twopi ') |
|---|
| 48 | do i = 1, 6 |
|---|
| 49 | orbit(i) = orbit0(i) |
|---|
| 50 | disp(i) = disp0(i) |
|---|
| 51 | enddo |
|---|
| 52 | call dzero(emit_v, 3) |
|---|
| 53 | call dzero(nemit_v, 3) |
|---|
| 54 | call dzero(bmax, 9) |
|---|
| 55 | call dzero(gmax, 9) |
|---|
| 56 | call dzero(dismax, 4) |
|---|
| 57 | call dzero(tunes, 3) |
|---|
| 58 | call dzero(sig_v, 4) |
|---|
| 59 | call dzero(pdamp, 3) |
|---|
| 60 | u0 = zero |
|---|
| 61 | !---- Find eigenvectors at initial position. |
|---|
| 62 | if (m66sta(rt)) then |
|---|
| 63 | call laseig(rt, reval, aival, em) |
|---|
| 64 | stabt = .false. |
|---|
| 65 | ! print '('' Static map, eigenvalues:'',(/1X,2E15.8))', |
|---|
| 66 | ! + (reval(i), aival(i), i = 1, 4) |
|---|
| 67 | else |
|---|
| 68 | call ladeig(rt, reval, aival, em) |
|---|
| 69 | stabt = reval(5)**2 + aival(5)**2 .le. tol .and. & |
|---|
| 70 | reval(6)**2 + aival(6)**2 .le. tol |
|---|
| 71 | ! print '('' Static map, eigenvalues:'',(/1X,2E15.8))', |
|---|
| 72 | ! + (reval(i), aival(i), i = 1, 6) |
|---|
| 73 | endif |
|---|
| 74 | stabx = reval(1)**2 + aival(1)**2 .le. tol .and. & |
|---|
| 75 | reval(2)**2 + aival(2)**2 .le. tol |
|---|
| 76 | staby = reval(3)**2 + aival(3)**2 .le. tol .and. & |
|---|
| 77 | reval(4)**2 + aival(4)**2 .le. tol |
|---|
| 78 | |
|---|
| 79 | !---- Maximum extents. |
|---|
| 80 | do j = 1, 3 |
|---|
| 81 | j1 = 2 * j -1 |
|---|
| 82 | j2 = 2 * j |
|---|
| 83 | do k = 1, 3 |
|---|
| 84 | k1 = 2 * k - 1 |
|---|
| 85 | k2 = 2 * k |
|---|
| 86 | bmax(j,k) = em(j1,k1) * em(j1,k1) + em(j1,k2) * em(j1,k2) |
|---|
| 87 | gmax(j,k) = em(j2,k1) * em(j2,k1) + em(j2,k2) * em(j2,k2) |
|---|
| 88 | enddo |
|---|
| 89 | enddo |
|---|
| 90 | arad = get_value('probe ','arad ') |
|---|
| 91 | betas = get_value('probe ','beta ') |
|---|
| 92 | gammas = get_value('probe ','gamma ') |
|---|
| 93 | cg = arad * gammas**3 / three |
|---|
| 94 | frad = get_value('probe ','radiate ') .ne. zero |
|---|
| 95 | !---- Initialize damping calculation. |
|---|
| 96 | if (frad .and. stabt) then |
|---|
| 97 | sum(1) = zero |
|---|
| 98 | sum(2) = zero |
|---|
| 99 | sum(3) = zero |
|---|
| 100 | sumu0 = zero |
|---|
| 101 | call m66one(rd) |
|---|
| 102 | endif |
|---|
| 103 | call dzero(tt,216) |
|---|
| 104 | call m66one(rt) |
|---|
| 105 | eflag = 0 |
|---|
| 106 | suml = zero |
|---|
| 107 | bbd_cnt=0 |
|---|
| 108 | bbd_flag=1 |
|---|
| 109 | |
|---|
| 110 | i = restart_sequ() |
|---|
| 111 | 10 continue |
|---|
| 112 | bbd_pos=i |
|---|
| 113 | code = node_value('mad8_type ') |
|---|
| 114 | if(code.eq.39) code=15 |
|---|
| 115 | if(code.eq.38) code=24 |
|---|
| 116 | el = node_value('l ') |
|---|
| 117 | n_align = node_al_errors(al_errors) |
|---|
| 118 | if (n_align.ne.0) then |
|---|
| 119 | call dcopy(orbit, orbit2, 6) |
|---|
| 120 | call tmali1(orbit2,al_errors,betas,gammas,orbit,re) |
|---|
| 121 | if (.not. stabt) call m66byv(re, disp, disp) |
|---|
| 122 | call m66mpy(re, em, em) |
|---|
| 123 | if (frad .and. stabt) call m66mpy(re, rd, rd) |
|---|
| 124 | endif |
|---|
| 125 | !*---- Keep orbit at entrance. |
|---|
| 126 | call dcopy(orbit, orbit1, 6) |
|---|
| 127 | !---- Element matrix and length. |
|---|
| 128 | call tmmap(code,.true.,.true.,orbit,fmap,ek,re,te) |
|---|
| 129 | if (fmap) then |
|---|
| 130 | !---- Advance dispersion. |
|---|
| 131 | if (.not. stabt) then |
|---|
| 132 | call m66byv(re, disp, disp) |
|---|
| 133 | do j = 1, 4 |
|---|
| 134 | dismax(j) = max(abs(disp(j)),dismax(j)) |
|---|
| 135 | enddo |
|---|
| 136 | endif |
|---|
| 137 | !---- Radiation damping. |
|---|
| 138 | call m66mpy(re, em, em2) |
|---|
| 139 | if (frad .and. stabt) then |
|---|
| 140 | call emdamp(code, deltap, em, em2, orbit1, orbit, re) |
|---|
| 141 | call m66mpy(re, rd, rd) |
|---|
| 142 | endif |
|---|
| 143 | call dcopy(em2, em, 36) |
|---|
| 144 | |
|---|
| 145 | !---- Compute beta and gamma. |
|---|
| 146 | do j = 1, 3 |
|---|
| 147 | j1 = 2 * j -1 |
|---|
| 148 | j2 = 2 * j |
|---|
| 149 | do k = 1, 3 |
|---|
| 150 | k1 = 2 * k - 1 |
|---|
| 151 | k2 = 2 * k |
|---|
| 152 | bx = em(j1,k1) * em(j1,k1) + em(j1,k2) * em(j1,k2) |
|---|
| 153 | bmax(j,k) = max(bmax(j,k),bx) |
|---|
| 154 | gx = em(j2,k1) * em(j2,k1) + em(j2,k2) * em(j2,k2) |
|---|
| 155 | gmax(j,k) = max(gmax(j,k),gx) |
|---|
| 156 | enddo |
|---|
| 157 | enddo |
|---|
| 158 | suml = suml + el |
|---|
| 159 | endif |
|---|
| 160 | if (n_align.ne.0) then |
|---|
| 161 | call dcopy(orbit, orbit2, 6) |
|---|
| 162 | call tmali2(el,orbit2,al_errors,betas,gammas,orbit,re) |
|---|
| 163 | if (.not. stabt) call m66byv(re, disp, disp) |
|---|
| 164 | call m66mpy(re, em, em) |
|---|
| 165 | if (frad .and. stabt) call m66mpy(re, rd, rd) |
|---|
| 166 | endif |
|---|
| 167 | if (advance_node().ne.0) then |
|---|
| 168 | i=i+1 |
|---|
| 169 | goto 10 |
|---|
| 170 | endif |
|---|
| 171 | bbd_flag=0 |
|---|
| 172 | !---- Undamped tunes and beam extents. |
|---|
| 173 | qx = atan2(aival(1), reval(1)) / twopi |
|---|
| 174 | if (qx .lt. zero) qx = qx + one |
|---|
| 175 | qy = atan2(aival(3), reval(3)) / twopi |
|---|
| 176 | if (qy .lt. zero) qy = qy + one |
|---|
| 177 | qs = atan2(aival(5), reval(5)) / twopi |
|---|
| 178 | if (qs .lt. zero) qs = - qs |
|---|
| 179 | !---- Summary output. |
|---|
| 180 | call emsumm(rd,em,bmax,gmax,stabt,frad,u0,emit_v,nemit_v,tunes, & |
|---|
| 181 | sig_v,pdamp) |
|---|
| 182 | end subroutine emit |
|---|
| 183 | subroutine emsumm(rd,em,bmax,gmax,stabt,frad,u0,emit_v,nemit_v, & |
|---|
| 184 | tunes,sig_v,pdamp) |
|---|
| 185 | use emitfi |
|---|
| 186 | implicit none |
|---|
| 187 | |
|---|
| 188 | |
|---|
| 189 | !----------------------------------------------------------------------* |
|---|
| 190 | ! Purpose: * |
|---|
| 191 | ! Finish radiation damping calculation and print summary. * |
|---|
| 192 | ! Input: * |
|---|
| 193 | ! RD(6,6) (real) Damped one-turn transfer matrix. * |
|---|
| 194 | ! EM(6,6) (real) Undamped eigenvectors. * |
|---|
| 195 | ! BMAX(3,3) (real) Maximum extents of modes. * |
|---|
| 196 | ! GMAX(3,3) (real) Maximum divergences of modes. * |
|---|
| 197 | ! STABT (logical) anti-STATIC flag |
|---|
| 198 | ! FRAD (logical) radiation flag |
|---|
| 199 | ! Output: |
|---|
| 200 | ! u0 (real) Radiation loss per turn in GeV |
|---|
| 201 | ! emit_v (real) ex, ey, et (emittances) |
|---|
| 202 | ! nemit_v (real) exn, eyn, etn (normalised emitt., MAD style) |
|---|
| 203 | ! tunes (real) qx, qy, qs |
|---|
| 204 | ! pdamp (real) damping partition numbers |
|---|
| 205 | ! sig_v (real) sigx, sigy, sigt, sige |
|---|
| 206 | !----------------------------------------------------------------------* |
|---|
| 207 | integer j,j1,j2,k,k1,k2,iqpr2 |
|---|
| 208 | double precision arad,gammas,clg,etpr,expr,eypr,f0,sal,en0 |
|---|
| 209 | double precision amass, clight, hbar, freq0, u0, betas |
|---|
| 210 | double precision ten3p,tenp6,tenp9,three,twopi,one,two,zero,four |
|---|
| 211 | double precision ex,ey,et,exn,eyn,sigx,sigy,sige,sigt |
|---|
| 212 | double precision rd(6,6), em(6,6), bmax(3,3), gmax(3,3),pdamp(3) |
|---|
| 213 | double precision get_variable,get_value |
|---|
| 214 | double precision emit_v(3), nemit_v(3), tunes(3), sig_v(4) |
|---|
| 215 | logical stabt, frad |
|---|
| 216 | double precision reval(6), aival(6), alj(3), tau(3), tune(3) |
|---|
| 217 | double precision sigma(6,6), bstar(3,3), gstar(3,3), dummy(6,6) |
|---|
| 218 | parameter (zero = 0.d0, one = 1.d0, two = 2.d0, three = 3.0d0) |
|---|
| 219 | parameter (iqpr2 = 6, four = 4.d0) |
|---|
| 220 | parameter (ten3p = 1.0d3, tenp6 = 1.0d6, tenp9 = 1.0d9) |
|---|
| 221 | |
|---|
| 222 | call dzero(sigma,36) |
|---|
| 223 | ex=zero |
|---|
| 224 | ey=zero |
|---|
| 225 | et=zero |
|---|
| 226 | twopi=get_variable('twopi ') |
|---|
| 227 | clight = get_variable('clight ') |
|---|
| 228 | hbar = get_variable('hbar ') |
|---|
| 229 | arad = get_value('probe ','arad ') |
|---|
| 230 | betas = get_value('probe ','beta ') |
|---|
| 231 | gammas = get_value('probe ','gamma ') |
|---|
| 232 | amass = get_value('probe ','mass ') |
|---|
| 233 | freq0 = get_value('probe ','freq0 ') |
|---|
| 234 | !---- Synchrotron energy loss [GeV]. |
|---|
| 235 | if (stabt .and. frad) then |
|---|
| 236 | u0 = sumu0 |
|---|
| 237 | |
|---|
| 238 | !---- Tunes. |
|---|
| 239 | call ladeig(rd, reval, aival, dummy) |
|---|
| 240 | tune(1) = atan2(aival(1), reval(1)) / twopi |
|---|
| 241 | if (tune(1) .lt. zero) tune(1) = tune(1) + one |
|---|
| 242 | tune(2) = atan2(aival(3), reval(3)) / twopi |
|---|
| 243 | if (tune(2) .lt. zero) tune(2) = tune(2) + one |
|---|
| 244 | tune(3) = atan2(aival(5), reval(5)) / twopi |
|---|
| 245 | if (tune(3) .lt. zero) tune(3) = - tune(3) |
|---|
| 246 | |
|---|
| 247 | !---- Damping constants per turn. |
|---|
| 248 | alj(1) = - log(reval(1)**2 + aival(1)**2) / two |
|---|
| 249 | alj(2) = - log(reval(3)**2 + aival(3)**2) / two |
|---|
| 250 | alj(3) = - log(reval(5)**2 + aival(5)**2) / two |
|---|
| 251 | |
|---|
| 252 | !---- Damping partition numbers. |
|---|
| 253 | en0 = get_value('probe ','energy ') |
|---|
| 254 | sal = two * en0 / u0 |
|---|
| 255 | pdamp(1) = alj(1) * sal |
|---|
| 256 | pdamp(2) = alj(2) * sal |
|---|
| 257 | pdamp(3) = alj(3) * sal |
|---|
| 258 | !---- Emittances. |
|---|
| 259 | clg = ((55.d0 * hbar * clight) / (96.d0 * sqrt(three))) & |
|---|
| 260 | * ((arad * gammas**5) / amass) |
|---|
| 261 | ex = clg * sum(1) / alj(1) |
|---|
| 262 | ey = clg * sum(2) / alj(2) |
|---|
| 263 | et = clg * sum(3) / alj(3) |
|---|
| 264 | |
|---|
| 265 | !---- Damping constants per second and damping times. |
|---|
| 266 | f0 = freq0 * tenp6 |
|---|
| 267 | alj(1) = abs(alj(1) * f0) |
|---|
| 268 | alj(2) = abs(alj(2) * f0) |
|---|
| 269 | alj(3) = abs(alj(3) * f0) |
|---|
| 270 | tau(1) = one / alj(1) |
|---|
| 271 | tau(2) = one / alj(2) |
|---|
| 272 | tau(3) = one / alj(3) |
|---|
| 273 | endif |
|---|
| 274 | |
|---|
| 275 | !---- TRANSPORT sigma matrix. |
|---|
| 276 | call emce2i(stabt, em, ex, ey, et, sigma) |
|---|
| 277 | |
|---|
| 278 | !---- Extents at interaction point. |
|---|
| 279 | do j = 1, 3 |
|---|
| 280 | j1 = 2 * j -1 |
|---|
| 281 | j2 = 2 * j |
|---|
| 282 | do k = 1, 3 |
|---|
| 283 | k1 = 2 * k - 1 |
|---|
| 284 | k2 = 2 * k |
|---|
| 285 | bstar(j,k) = em(j1,k1) * em(j1,k1) + em(j1,k2) * em(j1,k2) |
|---|
| 286 | gstar(j,k) = em(j2,k1) * em(j2,k1) + em(j2,k2) * em(j2,k2) |
|---|
| 287 | enddo |
|---|
| 288 | enddo |
|---|
| 289 | |
|---|
| 290 | exn = ex * four * betas * gammas |
|---|
| 291 | eyn = ey * four * betas * gammas |
|---|
| 292 | sigx = sqrt(abs(sigma(1,1))) |
|---|
| 293 | sigy = sqrt(abs(sigma(3,3))) |
|---|
| 294 | if (sigma(5,5) .gt. zero .or. sigma(6,6) .gt. zero) then |
|---|
| 295 | sigt = sqrt(abs(sigma(5,5))) |
|---|
| 296 | sige = sqrt(abs(sigma(6,6))) |
|---|
| 297 | else |
|---|
| 298 | sigt = zero |
|---|
| 299 | sige = zero |
|---|
| 300 | endif |
|---|
| 301 | tunes(1) = qx |
|---|
| 302 | tunes(2) = qy |
|---|
| 303 | tunes(3) = qs |
|---|
| 304 | emit_v(1) = ex |
|---|
| 305 | emit_v(2) = ey |
|---|
| 306 | emit_v(3) = et |
|---|
| 307 | nemit_v(1) = exn |
|---|
| 308 | nemit_v(2) = eyn |
|---|
| 309 | sig_v(1) = sigx |
|---|
| 310 | sig_v(2) = sigy |
|---|
| 311 | sig_v(3) = sigt |
|---|
| 312 | sig_v(4) = sige |
|---|
| 313 | !---- Summary output; header and global parameters. |
|---|
| 314 | !---- Dynamic case. |
|---|
| 315 | expr = ex * tenp6 |
|---|
| 316 | eypr = ey * tenp6 |
|---|
| 317 | etpr = et * tenp6 |
|---|
| 318 | if (stabt) then |
|---|
| 319 | if (frad) write (iqpr2, 910) ten3p * u0 |
|---|
| 320 | write (iqpr2, 920) 1, 2, 3 |
|---|
| 321 | write (iqpr2, 930) qx, qy, qs |
|---|
| 322 | if (frad) write (iqpr2, 940) tune |
|---|
| 323 | write (iqpr2, 950) ((bstar(j,k), j = 1, 3), k = 1, 3), & |
|---|
| 324 | ((gstar(j,k), j = 1, 3), k = 1, 3), & |
|---|
| 325 | ((bmax(j,k), j = 1, 3), k = 1, 3), & |
|---|
| 326 | ((gmax(j,k), j = 1, 3), k = 1, 3) |
|---|
| 327 | if (frad) then |
|---|
| 328 | write (iqpr2, 960) pdamp, alj, (tau(j), j = 1, 3), & |
|---|
| 329 | expr, eypr, etpr |
|---|
| 330 | endif |
|---|
| 331 | else |
|---|
| 332 | write (iqpr2, 920) 1, 2 |
|---|
| 333 | write (iqpr2, 930) qx, qy |
|---|
| 334 | write (iqpr2, 970) ((bstar(j,k), j = 1, 2), k = 1, 2), & |
|---|
| 335 | ((gstar(j,k), j = 1, 2), k = 1, 2), & |
|---|
| 336 | ((bmax(j,k), j = 1, 2), k = 1, 2), & |
|---|
| 337 | ((gmax(j,k), j = 1, 2), k = 1, 2) |
|---|
| 338 | endif |
|---|
| 339 | |
|---|
| 340 | !---- RF system. |
|---|
| 341 | ! call enprrf |
|---|
| 342 | |
|---|
| 343 | 910 format(t6,'U0',t16,f14.6,' [MeV/turn]') |
|---|
| 344 | 920 format(' '/' ',t42,3(9x,'M o d e',3x,i1:)) |
|---|
| 345 | 930 format(' Fractional tunes',t30,'undamped',t42,3f20.8) |
|---|
| 346 | 940 format(' ',t30,'damped',t42,3f20.8) |
|---|
| 347 | 950 format(' '/' beta* [m]',t30,'x',t42,3e20.8/t30,'y',t42,3e20.8/ & |
|---|
| 348 | t30,'t',t42,3e20.8/ & |
|---|
| 349 | ' '/' gamma* [1/m]',t30,'px',t42,3e20.8/t30,'py',t42,3e20.8/ & |
|---|
| 350 | t30,'pt',t42,3e20.8/ & |
|---|
| 351 | ' '/' beta(max) [m]',t30,'x',t42,3e20.8/t30,'y',t42,3e20.8/ & |
|---|
| 352 | t30,'t',t42,3e20.8/ & |
|---|
| 353 | ' '/' gamma(max) [1/m]',t30,'px',t42,3e20.8/t30,'py',t42,3e20.8/ & |
|---|
| 354 | t30,'pt',t42,3e20.8) |
|---|
| 355 | 960 format(' '/' Damping partition numbers',t42,3f20.8/ & |
|---|
| 356 | ' Damping constants [1/s]',t46,3e20.8/ & |
|---|
| 357 | ' Damping times [s]',t46,3e20.8/ & |
|---|
| 358 | ' Emittances [pi micro m]',t42,3e20.8) |
|---|
| 359 | 970 format(' '/' beta* [m]',t30,'x',t42,2e20.8/t30,'y',t42,2e20.8/ & |
|---|
| 360 | ' '/' gamma* [1/m]',t30,'px',t42,2e20.8/t30,'py',t42,2e20.8/ & |
|---|
| 361 | ' '/' beta(max) [m]',t30,'x',t42,2e20.8/t30,'y',t42,2e20.8/ & |
|---|
| 362 | ' '/' gamma(max) [1/m]',t30,'px',t42,2e20.8/t30,'py',t42,2e20.8) |
|---|
| 363 | |
|---|
| 364 | end subroutine emsumm |
|---|
| 365 | subroutine emce2i(stabt, em, ex, ey, et, sigma) |
|---|
| 366 | implicit none |
|---|
| 367 | !----------------------------------------------------------------------* |
|---|
| 368 | ! Purpose: * |
|---|
| 369 | ! Convert eigenvectors to internal sigma matrix form. * |
|---|
| 370 | ! Input: * |
|---|
| 371 | ! EM(6,6) (real) Eigenvector matrix. * |
|---|
| 372 | ! EX (real) Horizontal emittance. * |
|---|
| 373 | ! EY (real) Vertical emittance. * |
|---|
| 374 | ! ET (real) Longitudinal emittance. * |
|---|
| 375 | ! Output: * |
|---|
| 376 | ! SIGMA(6,6)(real) Beam matrix in internal form. * |
|---|
| 377 | !----------------------------------------------------------------------* |
|---|
| 378 | integer j,k |
|---|
| 379 | double precision et, ex, ey, em(6,6), sigma(6,6) |
|---|
| 380 | logical stabt |
|---|
| 381 | |
|---|
| 382 | do j = 1, 6 |
|---|
| 383 | do k = 1, 6 |
|---|
| 384 | sigma(j,k) = & |
|---|
| 385 | ex * (em(j,1) * em(k,1) + em(j,2) * em(k,2)) + & |
|---|
| 386 | ey * (em(j,3) * em(k,3) + em(j,4) * em(k,4)) |
|---|
| 387 | if (stabt) then |
|---|
| 388 | sigma(j,k) = sigma(j,k) + & |
|---|
| 389 | et * (em(j,5) * em(k,5) + em(j,6) * em(k,6)) |
|---|
| 390 | endif |
|---|
| 391 | enddo |
|---|
| 392 | enddo |
|---|
| 393 | end subroutine emce2i |
|---|
| 394 | subroutine emdamp(code, deltap, em1, em2, orb1, orb2, re) |
|---|
| 395 | use twiss0fi |
|---|
| 396 | use emitfi |
|---|
| 397 | use twtrrfi |
|---|
| 398 | implicit none |
|---|
| 399 | |
|---|
| 400 | |
|---|
| 401 | !---------------------------------------------------------------------* |
|---|
| 402 | ! Purpose: * |
|---|
| 403 | ! Deal with radiation damping in an element. * |
|---|
| 404 | ! Input: * |
|---|
| 405 | ! code (int) MAD-8 element code * |
|---|
| 406 | ! deltap (real) momentum error * |
|---|
| 407 | ! EM1(6,6) (real) Matrix of eigenvectors at entrance. * |
|---|
| 408 | ! EM2(6,6) (real) Matrix of eigenvectors at exit. * |
|---|
| 409 | ! ORB1(6) (real) Orbit position at entrance. * |
|---|
| 410 | ! ORB2(6) (real) Orbit position at exit. * |
|---|
| 411 | ! Input/output: * |
|---|
| 412 | ! RE(6,6) (real) Transfer matrix for the element; changed on * |
|---|
| 413 | ! output to contain damping. * |
|---|
| 414 | !---------------------------------------------------------------------* |
|---|
| 415 | integer code |
|---|
| 416 | double precision deltap |
|---|
| 417 | double precision em1(6,6), em2(6,6), orb1(6), orb2(6), re(6,6) |
|---|
| 418 | double precision ten3m, ten6p, zero, half, one, two, three, four |
|---|
| 419 | double precision six, twelve |
|---|
| 420 | parameter (ten3m = 1.0d-3, ten6p = 1.0d+6) |
|---|
| 421 | parameter (zero = 0.0d0, half = 0.5d0) |
|---|
| 422 | parameter (one = 1.0d0, two = 2.0d0) |
|---|
| 423 | parameter (three = 3.0d0, twelve = 12.d0) |
|---|
| 424 | parameter (four = 4.0d0, six = 6.0d0) |
|---|
| 425 | integer i, j, ir, ii, n, n_ferr, iord, nn, ns, nd, nord |
|---|
| 426 | integer node_fd_errors |
|---|
| 427 | double precision rw(6,6), tw(6,6,6), ferror(2), rw0(6,6) |
|---|
| 428 | double precision normal(0:maxmul), skew(0:maxmul) |
|---|
| 429 | double precision vals(2,0:maxmul), field(2,0:maxmul) |
|---|
| 430 | double precision f_errors(0:50) |
|---|
| 431 | double precision o1(6), e1(6,6), o2(6), e2(6,6) |
|---|
| 432 | double precision x1, y1, t1, x2, y2, t2, px1, py1, pt1 |
|---|
| 433 | double precision px2, py2, pt2 |
|---|
| 434 | equivalence (x1, o1(1)), (px1, o1(2)) |
|---|
| 435 | equivalence (y1, o1(3)), (py1, o1(4)) |
|---|
| 436 | equivalence (t1, o1(3)), (pt1, o1(4)) |
|---|
| 437 | equivalence (x2, o2(1)), (px2, o2(2)) |
|---|
| 438 | equivalence (y2, o2(3)), (py2, o2(4)) |
|---|
| 439 | equivalence (t2, o2(3)), (pt2, o2(4)) |
|---|
| 440 | double precision el, node_value, tilt, bvk |
|---|
| 441 | double precision edg1, edg2, sk1, sk2, hgap, fint, sks, h, ct |
|---|
| 442 | double precision corr, hx, hy, hxx, hxy, hyy, h1, hcb1, hcbs1 |
|---|
| 443 | double precision tedg1, fact1, fact1x, rfac1, rfac1x, rfac1y |
|---|
| 444 | double precision h2, hcb2, tedg2, fact2, fact2x, rfac2 |
|---|
| 445 | double precision rfac2x, rfac2y, bi2gi2, betas, gammas |
|---|
| 446 | double precision get_value, e5sq1, e5sq2, e5sqs1, e5sqs2, x, y |
|---|
| 447 | double precision f1, f2, f1s, f2s, twon, str, st, pi, clight |
|---|
| 448 | double precision r1sq, r2sq, fh1, fh2, dr, di, drt, hcb |
|---|
| 449 | double precision rfv, rff, rfl, rfvlt, rffrq, rflag, time |
|---|
| 450 | double precision xkick, ykick, dpx, dpy, an, hyx, hcbs2,hbi |
|---|
| 451 | double precision sk3, rfac, rfacx, rfacy |
|---|
| 452 | double precision get_variable, fh |
|---|
| 453 | ! double precision sk0, sk0s |
|---|
| 454 | if (code .eq. 8) then |
|---|
| 455 | !--- multipole |
|---|
| 456 | el = node_value('lrad ') |
|---|
| 457 | else |
|---|
| 458 | el = node_value('l ') |
|---|
| 459 | endif |
|---|
| 460 | if (el .eq. zero.and.code .ne. 10) goto 500 |
|---|
| 461 | betas = get_value('probe ','beta ') |
|---|
| 462 | gammas = get_value('probe ','gamma ') |
|---|
| 463 | !---- Prepare data. |
|---|
| 464 | bvk = node_value('other_bv ') |
|---|
| 465 | !---- Switch on element type. |
|---|
| 466 | go to (500, 20, 30, 500, 50, 60, 70, 80, 500, 100, & |
|---|
| 467 | 500, 500, 500, 140, 150, 160, 500, 500, 500, 500, & |
|---|
| 468 | 500, 500, 500, 500, 500, 500, 500, 500, 500, 500, & |
|---|
| 469 | 500, 500, 500, 500, 500, 500, 500, 500, 150, 500), code |
|---|
| 470 | go to 500 |
|---|
| 471 | |
|---|
| 472 | !---- Dipole. |
|---|
| 473 | 20 continue |
|---|
| 474 | 30 continue |
|---|
| 475 | |
|---|
| 476 | ! FRS 16.11.2004 This is still in the intermediate spirit of k0&k0s |
|---|
| 477 | ! an = node_value('angle ') |
|---|
| 478 | ! sk0 = an / el |
|---|
| 479 | ! sk0s = node_value('k0s ') |
|---|
| 480 | ! if (sk0s .eq. zero) then |
|---|
| 481 | ! tilt = zero |
|---|
| 482 | ! else |
|---|
| 483 | ! tilt = atan2(sk0s, sk0) |
|---|
| 484 | ! sk0 = sqrt(sk0**2 + sk0s**2) |
|---|
| 485 | ! an = sk0 * el |
|---|
| 486 | ! endif |
|---|
| 487 | ! an = bvk * an |
|---|
| 488 | ! sk0 = bvk * sk0 |
|---|
| 489 | ! sk0s = bvk * sk0s |
|---|
| 490 | ! FRS 16.11.2004 here the correction |
|---|
| 491 | an = bvk * node_value('angle ') * el/node_value('l ') |
|---|
| 492 | tilt = -node_value('tilt ') |
|---|
| 493 | edg1 = bvk * node_value('e1 ') |
|---|
| 494 | edg2 = bvk * node_value('e2 ') |
|---|
| 495 | sk1 = bvk * node_value('k1 ') |
|---|
| 496 | sk2 = bvk * node_value('k2 ') |
|---|
| 497 | hgap = node_value('hgap ') |
|---|
| 498 | fint = node_value('fint ') |
|---|
| 499 | sks = zero |
|---|
| 500 | h = an / el |
|---|
| 501 | |
|---|
| 502 | !---- Refer orbit and eigenvectors to magnet midplane. |
|---|
| 503 | ct = cos(tilt) |
|---|
| 504 | st = sin(tilt) |
|---|
| 505 | o1(1) = ct * orb1(1) + st * orb1(3) |
|---|
| 506 | o1(2) = ct * orb1(2) + st * orb1(4) |
|---|
| 507 | o1(3) = - st * orb1(1) + ct * orb1(3) |
|---|
| 508 | o1(4) = - st * orb1(2) + ct * orb1(4) |
|---|
| 509 | o1(5) = orb1(5) |
|---|
| 510 | o1(6) = orb1(6) |
|---|
| 511 | o2(1) = ct * orb2(1) + st * orb2(3) |
|---|
| 512 | o2(2) = ct * orb2(2) + st * orb2(4) |
|---|
| 513 | o2(3) = - st * orb2(1) + ct * orb2(3) |
|---|
| 514 | o2(4) = - st * orb2(2) + ct * orb2(4) |
|---|
| 515 | o2(5) = orb2(5) |
|---|
| 516 | o2(6) = orb2(6) |
|---|
| 517 | do i = 1, 6 |
|---|
| 518 | e1(1,i) = ct * em1(1,i) + st * em1(3,i) |
|---|
| 519 | e1(2,i) = ct * em1(2,i) + st * em1(4,i) |
|---|
| 520 | e1(3,i) = - st * em1(1,i) + ct * em1(3,i) |
|---|
| 521 | e1(4,i) = - st * em1(2,i) + ct * em1(4,i) |
|---|
| 522 | e1(5,i) = em1(5,i) |
|---|
| 523 | e1(6,i) = em1(6,i) |
|---|
| 524 | e2(1,i) = ct * em2(1,i) + st * em2(3,i) |
|---|
| 525 | e2(2,i) = ct * em2(2,i) + st * em2(4,i) |
|---|
| 526 | e2(3,i) = - st * em2(1,i) + ct * em2(3,i) |
|---|
| 527 | e2(4,i) = - st * em2(2,i) + ct * em2(4,i) |
|---|
| 528 | e2(5,i) = em2(5,i) |
|---|
| 529 | e2(6,i) = em2(6,i) |
|---|
| 530 | enddo |
|---|
| 531 | !---- Move through orbit through fringing field; |
|---|
| 532 | ! Requested components of eigenvectors are not affected. |
|---|
| 533 | corr = (h + h) * hgap * fint |
|---|
| 534 | call m66one(rw) |
|---|
| 535 | call dzero(tw,216) |
|---|
| 536 | call tmfrng(.false.,h,sk1,edg1,zero,+one,corr,rw,tw) |
|---|
| 537 | call m66byv(rw,o1,o1) |
|---|
| 538 | call m66one(rw) |
|---|
| 539 | call dzero(tw,216) |
|---|
| 540 | call tmfrng(.false.,h,sk1,edg2,zero,-one,corr,rw,tw) |
|---|
| 541 | call dcopy(rw,rw0,36) |
|---|
| 542 | call m66inv(rw0,rw) |
|---|
| 543 | call m66byv(rw,o2,o2) |
|---|
| 544 | |
|---|
| 545 | !---- Local curvature and its derivatives, |
|---|
| 546 | ! Coefficients for damping matrix. |
|---|
| 547 | hx = sk1*x1 + sks*y1 + h + half*sk2 * (x1**2 - y1**2) |
|---|
| 548 | hy = sks*x1 - sk1*y1 - sk2*x1*y1 |
|---|
| 549 | hxx = sk1 + sk2*x1 |
|---|
| 550 | hxy = sks - sk2*y1 |
|---|
| 551 | hyx = hxy |
|---|
| 552 | hyy = - hxx |
|---|
| 553 | h1 = sqrt(hx**2 + hy**2) |
|---|
| 554 | hcb1 = h1**3 |
|---|
| 555 | hcbs1 = three*h1 * & |
|---|
| 556 | (hx * (hxx*px1 + hxy*py1) + hy * (hxy*px1 + hyy*py1)) |
|---|
| 557 | |
|---|
| 558 | tedg1 = tan(edg1) |
|---|
| 559 | fact1 = (one + h*x1) * (one - tedg1*x1) |
|---|
| 560 | fact1x = h - tedg1 - 2.0*h*tedg1*x1 |
|---|
| 561 | rfac1 = cg*el*h1**2*fact1 |
|---|
| 562 | rfac1x = cg*el * (two*(hx*hxx+hy*hyx)*fact1 + h1**2*fact1x) |
|---|
| 563 | rfac1y = cg*el * two*(hx*hxy+hy*hyy)*fact1 |
|---|
| 564 | |
|---|
| 565 | hx = sk1*x2 + sks*y2 + h + half*sk2 * (x2**2 - y2**2) |
|---|
| 566 | hy = sks*x2 - sk1*y2 - sk2*x2*y2 |
|---|
| 567 | hxx = sk1 + sk2*x2 |
|---|
| 568 | hxy = sks - sk2*y2 |
|---|
| 569 | hyx = hxy |
|---|
| 570 | hyy = - hxx |
|---|
| 571 | h2 = sqrt(hx**2 + hy**2) |
|---|
| 572 | hcb2 = h2**3 |
|---|
| 573 | hcbs2 = three*h2 * & |
|---|
| 574 | (hx * (hxx*px2 + hxy*py2) + hy * (hxy*px2 + hyy*py2)) |
|---|
| 575 | |
|---|
| 576 | tedg2 = tan(edg2) |
|---|
| 577 | fact2 = (one + h*x2) * (one - tedg2*x2) |
|---|
| 578 | fact2x = h - tedg2 - 2.0*h*tedg2*x2 |
|---|
| 579 | rfac2 = cg*el*h2**2*fact2 |
|---|
| 580 | rfac2x = cg*el * (two*(hx*hxx+hy*hyx)*fact2 + h2**2*fact2x) |
|---|
| 581 | rfac2y = cg*el * two*(hx*hxy+hy*hyy)*fact2 |
|---|
| 582 | |
|---|
| 583 | !---- Cubic integration over h**3 * E(i,5) * conjg(E(i,5)). |
|---|
| 584 | bi2gi2 = one / (betas * gammas)**2 |
|---|
| 585 | hbi = h / betas |
|---|
| 586 | do i = 1, 3 |
|---|
| 587 | ir = 2 * i - 1 |
|---|
| 588 | ii = 2 * i |
|---|
| 589 | |
|---|
| 590 | !---- E(i,5) * conjg(E(i,5)) and its derivative w.r.t. S. |
|---|
| 591 | e5sq1 = e1(5,ir)**2 + e1(5,ii)**2 |
|---|
| 592 | e5sq2 = e2(5,ir)**2 + e2(5,ii)**2 |
|---|
| 593 | e5sqs1 = two * (e1(5,ir) * (bi2gi2*e1(6,ir) - hbi*e1(1,ir)) & |
|---|
| 594 | + e1(5,ii) * (bi2gi2*e1(6,ii) - hbi*e1(1,ii))) |
|---|
| 595 | e5sqs2 = two * (e2(5,ir) * (bi2gi2*e2(6,ir) - hbi*e2(1,ir)) & |
|---|
| 596 | + e2(5,ii) * (bi2gi2*e2(6,ii) - hbi*e2(1,ii))) |
|---|
| 597 | |
|---|
| 598 | !---- Integrand and its derivative w.r.t. S. |
|---|
| 599 | f1 = hcb1 * e5sq1 |
|---|
| 600 | f2 = hcb2 * e5sq2 |
|---|
| 601 | f1s = hcbs1 * e5sq1 + hcb1 * e5sqs1 |
|---|
| 602 | f2s = hcbs2 * e5sq2 + hcb2 * e5sqs2 |
|---|
| 603 | |
|---|
| 604 | !---- Actual integration. |
|---|
| 605 | sum(i) = sum(i) + half * el * (f1 + f2) - & |
|---|
| 606 | el**2 * (f2s - f1s) / twelve |
|---|
| 607 | enddo |
|---|
| 608 | go to 77 |
|---|
| 609 | |
|---|
| 610 | !---- Quadrupole. |
|---|
| 611 | 50 continue |
|---|
| 612 | bvk = node_value('other_bv ') |
|---|
| 613 | sk1 = bvk * node_value('k1 ') |
|---|
| 614 | str = sk1 |
|---|
| 615 | n = 1 |
|---|
| 616 | twon = two |
|---|
| 617 | go to 75 |
|---|
| 618 | |
|---|
| 619 | !---- Sextupole. |
|---|
| 620 | 60 continue |
|---|
| 621 | bvk = node_value('other_bv ') |
|---|
| 622 | sk2 = bvk * node_value('k2 ') |
|---|
| 623 | str = sk2 / two |
|---|
| 624 | n = 2 |
|---|
| 625 | twon = four |
|---|
| 626 | go to 75 |
|---|
| 627 | |
|---|
| 628 | !---- Octupole. |
|---|
| 629 | 70 continue |
|---|
| 630 | bvk = node_value('other_bv ') |
|---|
| 631 | sk3 = bvk * node_value('k2 ') |
|---|
| 632 | str = sk3 / six |
|---|
| 633 | n = 3 |
|---|
| 634 | twon = six |
|---|
| 635 | |
|---|
| 636 | !---- Common to all pure multipoles. |
|---|
| 637 | 75 continue |
|---|
| 638 | call dcopy(orb1, o1, 6) |
|---|
| 639 | call dcopy(orb2, o2, 6) |
|---|
| 640 | call dcopy(em1, e1, 36) |
|---|
| 641 | call dcopy(em2, e2, 36) |
|---|
| 642 | |
|---|
| 643 | !---- Local curvature. |
|---|
| 644 | r1sq = orb1(1)**2 + orb1(3)**2 |
|---|
| 645 | r2sq = orb2(1)**2 + orb2(3)**2 |
|---|
| 646 | h1 = abs(str) * sqrt(r1sq)**n |
|---|
| 647 | h2 = abs(str) * sqrt(r2sq)**n |
|---|
| 648 | rfac = cg * str**2 * el |
|---|
| 649 | rfac1 = rfac * r1sq**n |
|---|
| 650 | rfac2 = rfac * r2sq**n |
|---|
| 651 | rfac1x = twon * rfac * r1sq**(n-1) * x1 |
|---|
| 652 | rfac2x = twon * rfac * r1sq**(n-1) * x2 |
|---|
| 653 | rfac1y = twon * rfac * r1sq**(n-1) * y1 |
|---|
| 654 | rfac2y = twon * rfac * r1sq**(n-1) * y2 |
|---|
| 655 | |
|---|
| 656 | !---- Trapezoidal integration over h**3 * E(k,5) * conjg(E(k,5)). |
|---|
| 657 | fh1 = half * el * h1**3 |
|---|
| 658 | fh2 = half * el * h2**3 |
|---|
| 659 | sum(1) = sum(1) + fh1 * (e1(5,1)**2 + e1(5,2)**2) & |
|---|
| 660 | + fh2 * (e2(5,1)**2 + e2(5,2)**2) |
|---|
| 661 | sum(2) = sum(2) + fh1 * (e1(5,3)**2 + e1(5,4)**2) & |
|---|
| 662 | + fh2 * (e2(5,3)**2 + e2(5,4)**2) |
|---|
| 663 | sum(3) = sum(3) + fh1 * (e1(5,5)**2 + e1(5,6)**2) & |
|---|
| 664 | + fh2 * (e2(5,5)**2 + e2(5,6)**2) |
|---|
| 665 | |
|---|
| 666 | !---- Damping matrices. |
|---|
| 667 | ! Code common to bending magnet and pure multipoles. |
|---|
| 668 | 77 call m66one(rw) |
|---|
| 669 | rw(2,1) = - rfac1x * (one + pt1) * px1 |
|---|
| 670 | rw(2,2) = one - rfac1 * (one + pt1) |
|---|
| 671 | rw(2,3) = - rfac1y * (one + pt1) * px1 |
|---|
| 672 | rw(2,6) = - rfac1 * px1 |
|---|
| 673 | rw(4,1) = - rfac1x * (one + pt1) * py1 |
|---|
| 674 | rw(4,3) = - rfac1y * (one + pt1) * py1 |
|---|
| 675 | rw(4,4) = one - rfac1 * (one + pt1) |
|---|
| 676 | rw(4,6) = - rfac1 * py1 |
|---|
| 677 | rw(6,1) = - rfac1x * (one + pt1)**2 |
|---|
| 678 | rw(6,3) = - rfac1y * (one + pt1)**2 |
|---|
| 679 | rw(6,6) = one - two * rfac1 * (one + pt1) |
|---|
| 680 | call m66mpy(re, rw, re) |
|---|
| 681 | |
|---|
| 682 | call m66one(rw) |
|---|
| 683 | rw(2,1) = - rfac2x * (one + pt2) * px2 |
|---|
| 684 | rw(2,2) = one - rfac2 * (one + pt2) |
|---|
| 685 | rw(2,3) = - rfac2y * (one + pt2) * px2 |
|---|
| 686 | rw(2,6) = - rfac2 * px2 |
|---|
| 687 | rw(4,1) = - rfac2x * (one + pt2) * py2 |
|---|
| 688 | rw(4,3) = - rfac2y * (one + pt2) * py2 |
|---|
| 689 | rw(4,4) = one - rfac2 * (one + pt2) |
|---|
| 690 | rw(4,6) = - rfac2 * py2 |
|---|
| 691 | rw(6,1) = - rfac2x * (one + pt2)**2 |
|---|
| 692 | rw(6,3) = - rfac2y * (one + pt2)**2 |
|---|
| 693 | rw(6,6) = one - two * rfac2 * (one + pt2) |
|---|
| 694 | call m66mpy(rw, re, re) |
|---|
| 695 | go to 500 |
|---|
| 696 | |
|---|
| 697 | !---- Thin multipoles, EL is the fictitious length for radiation. |
|---|
| 698 | 80 continue |
|---|
| 699 | !---- Multipole components. |
|---|
| 700 | call dzero(f_errors,maxferr+1) |
|---|
| 701 | n_ferr = node_fd_errors(f_errors) |
|---|
| 702 | bvk = node_value('other_bv ') |
|---|
| 703 | call dzero(normal,maxmul+1) |
|---|
| 704 | call dzero(skew,maxmul+1) |
|---|
| 705 | call get_node_vector('knl ',nn,normal) |
|---|
| 706 | call get_node_vector('ksl ',ns,skew) |
|---|
| 707 | call dzero(vals,2*(maxmul+1)) |
|---|
| 708 | do iord = 0, nn |
|---|
| 709 | vals(1,iord) = normal(iord) |
|---|
| 710 | enddo |
|---|
| 711 | do iord = 0, ns |
|---|
| 712 | vals(2,iord) = skew(iord) |
|---|
| 713 | enddo |
|---|
| 714 | |
|---|
| 715 | !---- Field error vals. |
|---|
| 716 | call dzero(field,2*(maxmul+1)) |
|---|
| 717 | if (n_ferr .gt. 0) then |
|---|
| 718 | call dcopy(f_errors,field,n_ferr) |
|---|
| 719 | endif |
|---|
| 720 | nd = 2 * max(nn, ns, n_ferr/2-1) |
|---|
| 721 | |
|---|
| 722 | !---- Other components and errors. |
|---|
| 723 | nord = 0 |
|---|
| 724 | do iord = 0, nd/2 |
|---|
| 725 | do j = 1, 2 |
|---|
| 726 | field(j,iord) = bvk * (vals(j,iord) + field(j,iord)) & |
|---|
| 727 | / (one + deltap) |
|---|
| 728 | if (field(j,iord) .ne. zero) nord = iord |
|---|
| 729 | enddo |
|---|
| 730 | enddo |
|---|
| 731 | |
|---|
| 732 | !---- Track orbit. |
|---|
| 733 | x = orb1(1) |
|---|
| 734 | y = orb1(3) |
|---|
| 735 | |
|---|
| 736 | !---- Multipole kick. |
|---|
| 737 | dr = zero |
|---|
| 738 | di = zero |
|---|
| 739 | do iord = nord, 0, -1 |
|---|
| 740 | drt = (dr * x - di * y) / float(iord+1) + field(1,iord) |
|---|
| 741 | di = (dr * y + di * x) / float(iord+1) + field(2,iord) |
|---|
| 742 | dr = drt |
|---|
| 743 | enddo |
|---|
| 744 | |
|---|
| 745 | !---- H is local "curvature" due to multipole kick. |
|---|
| 746 | h = sqrt(dr**2 + di**2) / el |
|---|
| 747 | hcb = half * el * h**3 |
|---|
| 748 | sum(1) = sum(1) + hcb * & |
|---|
| 749 | (em1(5,1)**2 + em1(5,2)**2 + em2(5,1)**2 + em2(5,2)**2) |
|---|
| 750 | sum(2) = sum(2) + hcb * & |
|---|
| 751 | (em1(5,3)**2 + em1(5,4)**2 + em2(5,3)**2 + em2(5,4)**2) |
|---|
| 752 | sum(3) = sum(3) + hcb * & |
|---|
| 753 | (em1(5,5)**2 + em1(5,6)**2 + em2(5,5)**2 + em2(5,6)**2) |
|---|
| 754 | |
|---|
| 755 | !---- Damping matrix, is the same at both ends. |
|---|
| 756 | rfac = cg * (dr**2 + di**2) / el |
|---|
| 757 | rfacx = cg * (- dr * re(2,1) + di * re(4,1)) / el |
|---|
| 758 | rfacy = cg * (- dr * re(2,3) + di * re(4,3)) / el |
|---|
| 759 | |
|---|
| 760 | call m66one(rw) |
|---|
| 761 | rw(2,1) = - rfacx * (one + orb1(6)) * orb1(2) |
|---|
| 762 | rw(2,2) = one - rfac * (one + orb1(6)) |
|---|
| 763 | rw(2,3) = - rfacy * (one + orb1(6)) * orb1(2) |
|---|
| 764 | rw(2,6) = - rfac * orb1(2) |
|---|
| 765 | rw(4,1) = - rfacx * (one + orb1(6)) * orb1(4) |
|---|
| 766 | rw(4,3) = - rfacy * (one + orb1(6)) * orb1(4) |
|---|
| 767 | rw(4,4) = one - rfac * (one + orb1(6)) |
|---|
| 768 | rw(4,6) = - rfac * orb1(4) |
|---|
| 769 | rw(6,1) = - rfacx * (one + orb1(6)) |
|---|
| 770 | rw(6,3) = - rfacy * (one + orb1(6)) |
|---|
| 771 | rw(6,6) = one - two * rfac * (one + orb1(6)) |
|---|
| 772 | call m66mpy(re, rw, re) |
|---|
| 773 | call m66mpy(rw, re, re) |
|---|
| 774 | go to 500 |
|---|
| 775 | |
|---|
| 776 | !---- RF cavities. |
|---|
| 777 | 100 continue |
|---|
| 778 | rfv = node_value('volt ') |
|---|
| 779 | rff = node_value('freq ') |
|---|
| 780 | rfl = node_value('lag ') |
|---|
| 781 | pi = get_variable('pi ') |
|---|
| 782 | clight = get_variable('clight ') |
|---|
| 783 | rfvlt = ten3m * rfv |
|---|
| 784 | rffrq = rff * (ten6p * two * pi / clight) |
|---|
| 785 | rflag = two * pi * rfl |
|---|
| 786 | time = half * (orb1(5) + orb2(5)) |
|---|
| 787 | sumu0 = sumu0 + rfvlt * sin(rflag - rffrq * time) |
|---|
| 788 | go to 500 |
|---|
| 789 | |
|---|
| 790 | !---- Orbit correctors. |
|---|
| 791 | 140 continue |
|---|
| 792 | 150 continue |
|---|
| 793 | 160 continue |
|---|
| 794 | n_ferr = node_fd_errors(f_errors) |
|---|
| 795 | do i = 1, 2 |
|---|
| 796 | ferror(i) = zero |
|---|
| 797 | enddo |
|---|
| 798 | if (n_ferr .gt. 0) call dcopy(f_errors, ferror, min(2, n_ferr)) |
|---|
| 799 | if(code.eq.14) then |
|---|
| 800 | xkick=bvk*(node_value('kick ')+node_value('chkick ')+ & |
|---|
| 801 | ferror(1)) |
|---|
| 802 | ykick=zero |
|---|
| 803 | else if(code.eq.15.or.code.eq.39) then |
|---|
| 804 | xkick=bvk*(node_value('hkick ')+node_value('chkick ')+ & |
|---|
| 805 | ferror(1)) |
|---|
| 806 | ykick=bvk*(node_value('vkick ')+node_value('cvkick ')+ & |
|---|
| 807 | ferror(2)) |
|---|
| 808 | else if(code.eq.16) then |
|---|
| 809 | xkick=zero |
|---|
| 810 | ykick=bvk*(node_value('kick ')+node_value('cvkick ')+ & |
|---|
| 811 | ferror(2)) |
|---|
| 812 | else |
|---|
| 813 | xkick=zero |
|---|
| 814 | ykick=zero |
|---|
| 815 | endif |
|---|
| 816 | !---- Sum up total kicks. |
|---|
| 817 | dpx = xkick / (one + deltap) |
|---|
| 818 | dpy = ykick / (one + deltap) |
|---|
| 819 | |
|---|
| 820 | !---- Local curvature. |
|---|
| 821 | hx = abs(dpx) / el |
|---|
| 822 | hy = abs(dpy) / el |
|---|
| 823 | rfac = cg * (hx**2 + hx**2) * el |
|---|
| 824 | |
|---|
| 825 | !---- Trapezoidal integration over h**3*E(k,5)*E*(k,5). |
|---|
| 826 | fh = half * el * sqrt(hx**2 + hy**2)**3 |
|---|
| 827 | sum(1) = sum(1) + fh * & |
|---|
| 828 | (em1(5,1)**2 + em1(5,2)**2 + em2(5,1)**2 + em2(5,2)**2) |
|---|
| 829 | sum(2) = sum(2) + fh * & |
|---|
| 830 | (em1(5,3)**2 + em1(5,4)**2 + em2(5,3)**2 + em2(5,4)**2) |
|---|
| 831 | sum(3) = sum(3) + fh * & |
|---|
| 832 | (em1(5,5)**2 + em1(5,6)**2 + em2(5,5)**2 + em2(5,6)**2) |
|---|
| 833 | |
|---|
| 834 | !---- Damping matrices. |
|---|
| 835 | call m66one(rw) |
|---|
| 836 | rw(2,2) = one - rfac * (one + orb1(6)) |
|---|
| 837 | rw(2,6) = - rfac * orb1(2) |
|---|
| 838 | rw(4,4) = one - rfac * (one + orb1(6)) |
|---|
| 839 | rw(4,6) = - rfac * orb1(4) |
|---|
| 840 | rw(6,6) = one - two * rfac * (one + orb1(6)) |
|---|
| 841 | call m66mpy(re, rw, re) |
|---|
| 842 | |
|---|
| 843 | call m66one(rw) |
|---|
| 844 | rw(2,2) = one - rfac * (one + orb2(6)) |
|---|
| 845 | rw(2,6) = - rfac * orb2(2) |
|---|
| 846 | rw(4,4) = one - rfac * (one + orb2(6)) |
|---|
| 847 | rw(4,6) = - rfac * orb2(4) |
|---|
| 848 | rw(6,6) = one - two * rfac * (one + orb2(6)) |
|---|
| 849 | call m66mpy(rw, re, re) |
|---|
| 850 | 500 continue |
|---|
| 851 | |
|---|
| 852 | end subroutine emdamp |
|---|
![(please configure the [header_logo] section in trac.ini)](/PSPA/chrome/site/pspa.gif)
