| 1 | LOGICAL FUNCTION is_drift() |
|---|
| 2 | double precision node_value, code |
|---|
| 3 | code = node_value('mad8_type ') |
|---|
| 4 | is_drift = code.eq.1; |
|---|
| 5 | END FUNCTION is_drift |
|---|
| 6 | |
|---|
| 7 | LOGICAL FUNCTION is_matrix() |
|---|
| 8 | double precision node_value, code |
|---|
| 9 | code = node_value('mad8_type ') |
|---|
| 10 | is_matrix = code.eq.4; |
|---|
| 11 | END FUNCTION is_matrix |
|---|
| 12 | |
|---|
| 13 | LOGICAL FUNCTION is_quad() |
|---|
| 14 | double precision node_value, code |
|---|
| 15 | code = node_value('mad8_type ') |
|---|
| 16 | is_quad = code.eq.5; |
|---|
| 17 | END FUNCTION is_quad |
|---|
| 18 | |
|---|
| 19 | LOGICAL FUNCTION is_thin() |
|---|
| 20 | double precision node_value, el, zero |
|---|
| 21 | parameter(zero=0d0) |
|---|
| 22 | el = node_value('l ') |
|---|
| 23 | is_thin = el.eq.zero; |
|---|
| 24 | END FUNCTION is_thin |
|---|
| 25 | |
|---|
| 26 | subroutine trrun(switch,turns,orbit0,rt,part_id,last_turn, & |
|---|
| 27 | last_pos,z,dxt,dyt,last_orbit,eigen,coords,e_flag,code_buf,l_buf) |
|---|
| 28 | |
|---|
| 29 | use bbfi |
|---|
| 30 | use twiss0fi |
|---|
| 31 | use name_lenfi |
|---|
| 32 | use trackfi |
|---|
| 33 | implicit none |
|---|
| 34 | |
|---|
| 35 | !----------------------------------------------------------------------* |
|---|
| 36 | ! Purpose: * |
|---|
| 37 | ! Interface RUN and DYNAP command to tracking routine * |
|---|
| 38 | ! * |
|---|
| 39 | !-- Input: * |
|---|
| 40 | ! switch (int) 1: RUN, 2: DYNAP fastune * |
|---|
| 41 | ! turns (int) number of turns to track * |
|---|
| 42 | ! orbit0 (dp. array) start of closed orbit * |
|---|
| 43 | ! rt (dp. matrix) one-turn matrix * |
|---|
| 44 | !-- Output: * |
|---|
| 45 | ! eigen dp(6,6) eigenvector matrix * |
|---|
| 46 | ! coords dp(6,0:turns,npart) (only switch > 1) particle coords. * |
|---|
| 47 | ! e_flag (int) (only switch > 1) 0: OK, else part. lost * |
|---|
| 48 | !-- Buffers: * |
|---|
| 49 | ! part_id int(npart) * |
|---|
| 50 | ! last_turn int(npart) * |
|---|
| 51 | ! last_pos dp(npart) * |
|---|
| 52 | ! z dp(6,npart) * |
|---|
| 53 | ! last_orbit dp(6,npart) * |
|---|
| 54 | ! code_buf int(nelem) local mad-8 code storage * |
|---|
| 55 | ! l_buf dp(nelem) local length storage * |
|---|
| 56 | !----------------------------------------------------------------------* |
|---|
| 57 | logical onepass,onetable,last_out,info,aperflag,doupdate |
|---|
| 58 | integer j,code,restart_sequ,advance_node,node_al_errors,n_align, & |
|---|
| 59 | nlm,jmax,j_tot,turn,turns,i,k,get_option,ffile,SWITCH,nint,ndble, & |
|---|
| 60 | nchar,part_id(*),last_turn(*),char_l,segment, e_flag, nobs,lobs, & |
|---|
| 61 | int_arr(1),tot_segm,code_buf(*) |
|---|
| 62 | double precision tmp_d,orbit0(6),orbit(6),el,re(6,6),rt(6,6), & |
|---|
| 63 | al_errors(align_max),z(6,*),zz(6),dxt(*),dyt(*),eigen(6,6),sum, & |
|---|
| 64 | node_value,one, & |
|---|
| 65 | get_variable,last_pos(*),last_orbit(6,*),maxaper(6),get_value, & |
|---|
| 66 | zero,obs_orb(6),coords(6,0:turns,*),l_buf(*),deltap |
|---|
| 67 | parameter(zero=0d0,one=1d0) |
|---|
| 68 | character(12) tol_a, char_a |
|---|
| 69 | !hbu |
|---|
| 70 | double precision spos |
|---|
| 71 | !hbu |
|---|
| 72 | character(4) vec_names(7) |
|---|
| 73 | !hbu |
|---|
| 74 | character(name_len) el_name |
|---|
| 75 | character(120) msg |
|---|
| 76 | data tol_a,char_a / 'maxaper ', ' ' / |
|---|
| 77 | !hbu |
|---|
| 78 | data vec_names / 'x', 'px', 'y', 'py', 't', 'pt','s' / |
|---|
| 79 | |
|---|
| 80 | logical is_drift, is_thin, is_quad, is_matrix |
|---|
| 81 | |
|---|
| 82 | !-------added by Yipeng SUN 01-12-2008-------------- |
|---|
| 83 | deltap = get_value('probe ','deltap ') |
|---|
| 84 | |
|---|
| 85 | if(deltap.eq.0) then |
|---|
| 86 | onepass = get_option('onepass ') .ne. 0 |
|---|
| 87 | if(onepass) then |
|---|
| 88 | else |
|---|
| 89 | call trclor(orbit0) |
|---|
| 90 | endif |
|---|
| 91 | else |
|---|
| 92 | endif |
|---|
| 93 | !-------added by Yipeng SUN 01-12-2008-------------- |
|---|
| 94 | |
|---|
| 95 | !---- AK 2006 04 23 |
|---|
| 96 | !---- This version of trrun.F gets rid of all problems concerning delta_p |
|---|
| 97 | !---- by eliminating any delta_p dependence and using full 6D formulae only!!! |
|---|
| 98 | !---- Only the parts of the code that deal with radiation effects still use |
|---|
| 99 | !---- the quantity delta_p |
|---|
| 100 | |
|---|
| 101 | ! print *,"madX::trrun.F" |
|---|
| 102 | ! print *," " |
|---|
| 103 | ! print *," AK special version 2006/04/23" |
|---|
| 104 | ! print *," =============================" |
|---|
| 105 | ! print *," Full 6D formulae internally only." |
|---|
| 106 | |
|---|
| 107 | |
|---|
| 108 | if(fsecarb) then |
|---|
| 109 | write (msg,*) 'Second order terms of arbitrary Matrix not '// & |
|---|
| 110 | 'allowed for tracking.' |
|---|
| 111 | call aawarn('trrun: ',msg) |
|---|
| 112 | return |
|---|
| 113 | endif |
|---|
| 114 | aperflag = .false. |
|---|
| 115 | e_flag = 0 |
|---|
| 116 | ! flag to avoid double entry of last line |
|---|
| 117 | last_out = .false. |
|---|
| 118 | onepass = get_option('onepass ') .ne. 0 |
|---|
| 119 | onetable = get_option('onetable ') .ne. 0 |
|---|
| 120 | info = get_option('info ') * get_option('warn ') .gt. 0 |
|---|
| 121 | if(onepass) call m66one(rt) |
|---|
| 122 | call trbegn(rt,eigen) |
|---|
| 123 | if (switch .eq. 1) then |
|---|
| 124 | ffile = get_value('run ', 'ffile ') |
|---|
| 125 | else |
|---|
| 126 | ffile = 1 |
|---|
| 127 | endif |
|---|
| 128 | ! for one_table case |
|---|
| 129 | segment = 0 |
|---|
| 130 | tot_segm = turns / ffile + 1 |
|---|
| 131 | if (mod(turns, ffile) .ne. 0) tot_segm = tot_segm + 1 |
|---|
| 132 | call trinicmd(switch,orbit0,eigen,jmax,z,turns,coords) |
|---|
| 133 | !--- jmax may be reduced by particle loss - keep number in j_tot |
|---|
| 134 | j_tot = jmax |
|---|
| 135 | !--- get vector of six coordinate maxapers (both RUN and DYNAP) |
|---|
| 136 | call comm_para(tol_a, nint, ndble, nchar, int_arr, maxaper, & |
|---|
| 137 | char_a, char_l) |
|---|
| 138 | !--- set particle id |
|---|
| 139 | do k=1,jmax |
|---|
| 140 | part_id(k) = k |
|---|
| 141 | enddo |
|---|
| 142 | !hbu--- init info for tables initial s position is 0 |
|---|
| 143 | !hbu initial s position is 0 |
|---|
| 144 | spos=0 |
|---|
| 145 | !hbu start of line, element 0 |
|---|
| 146 | nlm=0 |
|---|
| 147 | !hbu |
|---|
| 148 | el_name='start ' |
|---|
| 149 | !--- enter start coordinates in summary table |
|---|
| 150 | do i = 1,j_tot |
|---|
| 151 | tmp_d = i |
|---|
| 152 | call double_to_table_curr('tracksumm ', 'number ', tmp_d) |
|---|
| 153 | tmp_d = 0 |
|---|
| 154 | call double_to_table_curr('tracksumm ', 'turn ', tmp_d) |
|---|
| 155 | do j = 1, 6 |
|---|
| 156 | tmp_d = z(j,i) - orbit0(j) |
|---|
| 157 | call double_to_table_curr('tracksumm ', vec_names(j), tmp_d) |
|---|
| 158 | enddo |
|---|
| 159 | !hbu add s |
|---|
| 160 | call double_to_table_curr('tracksumm ',vec_names(7),spos) |
|---|
| 161 | call augment_count('tracksumm ') |
|---|
| 162 | enddo |
|---|
| 163 | !--- enter first turn, and possibly eigen in tables |
|---|
| 164 | if (switch .eq. 1) then |
|---|
| 165 | if (onetable) then |
|---|
| 166 | call track_pteigen(eigen) |
|---|
| 167 | !hbu add s, node id and name |
|---|
| 168 | call tt_putone(jmax, 0, tot_segm, segment, part_id, & |
|---|
| 169 | z, orbit0,spos,nlm,el_name) |
|---|
| 170 | else |
|---|
| 171 | do i = 1, jmax |
|---|
| 172 | !hbu |
|---|
| 173 | call tt_puttab(part_id(i), 0, 1, z(1,i), orbit0,spos) |
|---|
| 174 | enddo |
|---|
| 175 | endif |
|---|
| 176 | endif |
|---|
| 177 | !---- Initialize kinematics and orbit |
|---|
| 178 | bet0 = get_value('beam ','beta ') |
|---|
| 179 | betas = get_value('probe ','beta ') |
|---|
| 180 | gammas= get_value('probe ','gamma ') |
|---|
| 181 | bet0i = one / bet0 |
|---|
| 182 | beti = one / betas |
|---|
| 183 | dtbyds = get_value('probe ','dtbyds ') |
|---|
| 184 | deltas = get_variable('track_deltap ') |
|---|
| 185 | arad = get_value('probe ','arad ') |
|---|
| 186 | dorad = get_value('probe ','radiate ') .ne. 0 |
|---|
| 187 | dodamp = get_option('damp ') .ne. 0 |
|---|
| 188 | dorand = get_option('quantum ') .ne. 0 |
|---|
| 189 | call dcopy(orbit0,orbit,6) |
|---|
| 190 | |
|---|
| 191 | doupdate = get_option('update ') .ne. 0 |
|---|
| 192 | |
|---|
| 193 | !--- loop over turns |
|---|
| 194 | nobs = 0 |
|---|
| 195 | do turn = 1, turns |
|---|
| 196 | |
|---|
| 197 | if (doupdate) call trupdate(turn) |
|---|
| 198 | |
|---|
| 199 | j = restart_sequ() |
|---|
| 200 | nlm = 0 |
|---|
| 201 | sum=zero |
|---|
| 202 | !--- loop over nodes |
|---|
| 203 | 10 continue |
|---|
| 204 | bbd_pos=j |
|---|
| 205 | if (turn .eq. 1) then |
|---|
| 206 | code = node_value('mad8_type ') |
|---|
| 207 | if(code.eq.39) code=15 |
|---|
| 208 | if(code.eq.38) code=24 |
|---|
| 209 | el = node_value('l ') |
|---|
| 210 | code_buf(nlm+1) = code |
|---|
| 211 | l_buf(nlm+1) = el |
|---|
| 212 | !hbu get current node name |
|---|
| 213 | call element_name(el_name,len(el_name)) |
|---|
| 214 | if (.not.(is_drift() .or. is_thin() .or. is_quad() .or. is_matrix())) then |
|---|
| 215 | print*," " |
|---|
| 216 | print*,"code: ",code," el: ",el," THICK ELEMENT FOUND" |
|---|
| 217 | sum = node_value('name ') |
|---|
| 218 | print*," " |
|---|
| 219 | print*,"Track dies nicely" |
|---|
| 220 | print*,"Thick lenses will get nowhere" |
|---|
| 221 | print*,"MAKETHIN will save you" |
|---|
| 222 | print*," " |
|---|
| 223 | print*," " |
|---|
| 224 | ! Better to use stop. If use return, irrelevant tracking output |
|---|
| 225 | ! is printed |
|---|
| 226 | stop |
|---|
| 227 | ! call aafail('TRRUN', & |
|---|
| 228 | ! '----element with length found : CONVERT STRUCTURE WITH '// & |
|---|
| 229 | ! 'MAKETHIN') |
|---|
| 230 | endif |
|---|
| 231 | else |
|---|
| 232 | el = l_buf(nlm+1) |
|---|
| 233 | code = code_buf(nlm+1) |
|---|
| 234 | endif |
|---|
| 235 | if (switch .eq. 1) then |
|---|
| 236 | nobs = node_value('obs_point ') |
|---|
| 237 | endif |
|---|
| 238 | !-------- Misalignment at beginning of element (from twissfs.f) |
|---|
| 239 | if (code .ne. 1) then |
|---|
| 240 | call dzero(al_errors, align_max) |
|---|
| 241 | n_align = node_al_errors(al_errors) |
|---|
| 242 | if (n_align .ne. 0) then |
|---|
| 243 | do i = 1, jmax |
|---|
| 244 | call dcopy(z(1,i),zz,6) |
|---|
| 245 | call tmali1(zz,al_errors, betas, gammas,z(1,i), re) |
|---|
| 246 | enddo |
|---|
| 247 | endif |
|---|
| 248 | endif |
|---|
| 249 | !-------- Track through element // suppress dxt 13.12.04 |
|---|
| 250 | call ttmap(code,el,z,jmax,dxt,dyt,sum,turn,part_id, & |
|---|
| 251 | last_turn,last_pos, last_orbit,aperflag,maxaper,al_errors,onepass) |
|---|
| 252 | !-------- Misalignment at end of element (from twissfs.f) |
|---|
| 253 | if (code .ne. 1) then |
|---|
| 254 | if (n_align .ne. 0) then |
|---|
| 255 | do i = 1, jmax |
|---|
| 256 | call dcopy(z(1,i),zz,6) |
|---|
| 257 | call tmali2(el,zz, al_errors, betas, gammas,z(1,i), re) |
|---|
| 258 | enddo |
|---|
| 259 | endif |
|---|
| 260 | endif |
|---|
| 261 | nlm = nlm+1 |
|---|
| 262 | if (nobs .gt. 0) then |
|---|
| 263 | call dzero(obs_orb,6) |
|---|
| 264 | call get_node_vector('obs_orbit ', lobs, obs_orb) |
|---|
| 265 | if (lobs .lt. 6) & |
|---|
| 266 | call aafail('TRACK', 'obs. point orbit not found') |
|---|
| 267 | if (onetable) then |
|---|
| 268 | !hbu |
|---|
| 269 | spos=sum |
|---|
| 270 | !hbu get current node name |
|---|
| 271 | call element_name(el_name,len(el_name)) |
|---|
| 272 | !hbu |
|---|
| 273 | call tt_putone(jmax, turn, tot_segm, segment, part_id, & |
|---|
| 274 | z, obs_orb,spos,nlm,el_name) |
|---|
| 275 | else |
|---|
| 276 | if (mod(turn, ffile) .eq. 0) then |
|---|
| 277 | do i = 1, jmax |
|---|
| 278 | !hbu add spos |
|---|
| 279 | call tt_puttab(part_id(i), turn, nobs, z(1,i), obs_orb, & |
|---|
| 280 | spos) |
|---|
| 281 | enddo |
|---|
| 282 | endif |
|---|
| 283 | endif |
|---|
| 284 | endif |
|---|
| 285 | if (advance_node().ne.0) then |
|---|
| 286 | j=j+1 |
|---|
| 287 | go to 10 |
|---|
| 288 | endif |
|---|
| 289 | !--- end of loop over nodes |
|---|
| 290 | if (switch .eq. 1) then |
|---|
| 291 | if (mod(turn, ffile) .eq. 0) then |
|---|
| 292 | if (turn .eq. turns) last_out = .true. |
|---|
| 293 | if (onetable) then |
|---|
| 294 | !hbu |
|---|
| 295 | spos=sum |
|---|
| 296 | !hbu spos added |
|---|
| 297 | call tt_putone(jmax, turn, tot_segm, segment, part_id, & |
|---|
| 298 | z, orbit0,spos,nlm,el_name) |
|---|
| 299 | else |
|---|
| 300 | do i = 1, jmax |
|---|
| 301 | !hbu |
|---|
| 302 | call tt_puttab(part_id(i), turn, 1, z(1,i), orbit0,spos) |
|---|
| 303 | enddo |
|---|
| 304 | endif |
|---|
| 305 | endif |
|---|
| 306 | else |
|---|
| 307 | do i = 1, jmax |
|---|
| 308 | do j = 1, 6 |
|---|
| 309 | coords(j,turn,i) = z(j,i) - orbit0(j) |
|---|
| 310 | enddo |
|---|
| 311 | enddo |
|---|
| 312 | endif |
|---|
| 313 | if (jmax .eq. 0 .or. (switch .gt. 1 .and. jmax .lt. j_tot)) & |
|---|
| 314 | goto 20 |
|---|
| 315 | if (switch .eq. 2 .and. info) then |
|---|
| 316 | if (mod(turn,100) .eq. 0) print *, 'turn :', turn |
|---|
| 317 | endif |
|---|
| 318 | enddo |
|---|
| 319 | !--- end of loop over turns |
|---|
| 320 | 20 continue |
|---|
| 321 | if (switch .gt. 1 .and. jmax .lt. j_tot) then |
|---|
| 322 | e_flag = 1 |
|---|
| 323 | return |
|---|
| 324 | endif |
|---|
| 325 | do i = 1, jmax |
|---|
| 326 | last_turn(part_id(i)) = min(turns, turn) |
|---|
| 327 | last_pos(part_id(i)) = sum |
|---|
| 328 | do j = 1, 6 |
|---|
| 329 | last_orbit(j,part_id(i)) = z(j,i) |
|---|
| 330 | enddo |
|---|
| 331 | enddo |
|---|
| 332 | turn = min(turn, turns) |
|---|
| 333 | !--- enter last turn in tables if not done already |
|---|
| 334 | if (.not. last_out) then |
|---|
| 335 | if (switch .eq. 1) then |
|---|
| 336 | if (onetable) then |
|---|
| 337 | !hbu |
|---|
| 338 | spos=sum |
|---|
| 339 | !hbu get current node name |
|---|
| 340 | call element_name(el_name,len(el_name)) |
|---|
| 341 | !hbu spos added |
|---|
| 342 | call tt_putone(jmax, turn, tot_segm, segment, part_id, & |
|---|
| 343 | z, orbit0,spos,nlm,el_name) |
|---|
| 344 | else |
|---|
| 345 | do i = 1, jmax |
|---|
| 346 | !hbu |
|---|
| 347 | call tt_puttab(part_id(i), turn, 1, z(1,i), orbit0,spos) |
|---|
| 348 | enddo |
|---|
| 349 | endif |
|---|
| 350 | endif |
|---|
| 351 | endif |
|---|
| 352 | !--- enter last turn in summary table |
|---|
| 353 | do i = 1,j_tot |
|---|
| 354 | tmp_d = i |
|---|
| 355 | call double_to_table_curr('tracksumm ', 'number ', tmp_d) |
|---|
| 356 | tmp_d = last_turn(i) |
|---|
| 357 | call double_to_table_curr('tracksumm ', 'turn ', tmp_d) |
|---|
| 358 | do j = 1, 6 |
|---|
| 359 | tmp_d = last_orbit(j,i) - orbit0(j) |
|---|
| 360 | call double_to_table_curr('tracksumm ', vec_names(j), tmp_d) |
|---|
| 361 | enddo |
|---|
| 362 | !hbu |
|---|
| 363 | spos=last_pos(i) |
|---|
| 364 | !hbu |
|---|
| 365 | call double_to_table_curr('tracksumm ',vec_names(7),spos) |
|---|
| 366 | call augment_count('tracksumm ') |
|---|
| 367 | enddo |
|---|
| 368 | end subroutine trrun |
|---|
| 369 | |
|---|
| 370 | |
|---|
| 371 | subroutine ttmap(code,el,track,ktrack,dxt,dyt,sum,turn,part_id, & |
|---|
| 372 | last_turn,last_pos, last_orbit,aperflag,maxaper,al_errors, & |
|---|
| 373 | onepass) |
|---|
| 374 | |
|---|
| 375 | use twtrrfi |
|---|
| 376 | use twiss0fi |
|---|
| 377 | use name_lenfi |
|---|
| 378 | implicit none |
|---|
| 379 | |
|---|
| 380 | !----------------------------------------------------------------------* |
|---|
| 381 | ! Purpose: * |
|---|
| 382 | ! - Test apertures * |
|---|
| 383 | ! - Track through thin lenses ONLY. * |
|---|
| 384 | ! * |
|---|
| 385 | ! Input/output: * |
|---|
| 386 | ! SUM (double) Accumulated length. * |
|---|
| 387 | ! TRACK(6,*)(double) Track coordinates: (X, PX, Y, PY, T, PT). * |
|---|
| 388 | ! NUMBER(*) (integer) Number of current track. * |
|---|
| 389 | ! KTRACK (integer) number of surviving tracks. * |
|---|
| 390 | !----------------------------------------------------------------------* |
|---|
| 391 | logical aperflag,fmap,onepass |
|---|
| 392 | integer turn,code,ktrack,part_id(*),last_turn(*),nn,jtrk, & |
|---|
| 393 | get_option |
|---|
| 394 | double precision apx,apy,apr,el,sum,node_value,track(6,*), & |
|---|
| 395 | last_pos(*),last_orbit(6,*),parvec(26),get_value,ct,tmp, & |
|---|
| 396 | aperture(maxnaper),one,maxaper(6), zero,al_errors(align_max),st, & |
|---|
| 397 | theta,ek(6),re(6,6),te(6,6,6),craporb(6),dxt(*),dyt(*),offset(2), & |
|---|
| 398 | offx,offy |
|---|
| 399 | character(name_len) aptype |
|---|
| 400 | parameter(zero = 0.d0, one=1d0) |
|---|
| 401 | |
|---|
| 402 | fmap=.false. |
|---|
| 403 | call dzero(ek,6) |
|---|
| 404 | call dzero(craporb,6) |
|---|
| 405 | call dzero(re,36) |
|---|
| 406 | call dzero(te,216) |
|---|
| 407 | |
|---|
| 408 | !---- Drift space |
|---|
| 409 | if(code.eq.1) then |
|---|
| 410 | call ttdrf(el,track,ktrack) |
|---|
| 411 | go to 502 |
|---|
| 412 | endif |
|---|
| 413 | |
|---|
| 414 | !---- Rotate trajectory before entry |
|---|
| 415 | theta = node_value('tilt ') |
|---|
| 416 | if (theta .ne. zero) then |
|---|
| 417 | st = sin(theta) |
|---|
| 418 | ct = cos(theta) |
|---|
| 419 | do jtrk = 1,ktrack |
|---|
| 420 | tmp = track(1,jtrk) |
|---|
| 421 | track(1,jtrk) = ct * tmp + st * track(3,jtrk) |
|---|
| 422 | track(3,jtrk) = ct * track(3,jtrk) - st * tmp |
|---|
| 423 | tmp = track(2,jtrk) |
|---|
| 424 | track(2,jtrk) = ct * tmp + st * track(4,jtrk) |
|---|
| 425 | track(4,jtrk) = ct * track(4,jtrk) - st * tmp |
|---|
| 426 | enddo |
|---|
| 427 | endif |
|---|
| 428 | |
|---|
| 429 | !---- Translation of particles // AK 23.02.2006 |
|---|
| 430 | !---- useful for combination of different transfer lines or rings |
|---|
| 431 | if(code.eq.36) then |
|---|
| 432 | if (onepass) then |
|---|
| 433 | call tttrans(track,ktrack) |
|---|
| 434 | go to 500 |
|---|
| 435 | endif |
|---|
| 436 | endif |
|---|
| 437 | |
|---|
| 438 | !---- Beam-beam, standard 4D, no aperture |
|---|
| 439 | if(code.eq.22) then |
|---|
| 440 | parvec(5)=get_value('probe ', 'arad ') |
|---|
| 441 | parvec(6)=node_value('charge ') * get_value('probe ', 'npart ') |
|---|
| 442 | parvec(7)=get_value('probe ','gamma ') |
|---|
| 443 | call ttbb(track, ktrack, parvec) |
|---|
| 444 | go to 500 |
|---|
| 445 | endif |
|---|
| 446 | |
|---|
| 447 | !---- Special colllimator aperture check taken out AK 20071211 |
|---|
| 448 | !!---- Collimator with elliptic aperture. |
|---|
| 449 | ! if(code.eq.20) then |
|---|
| 450 | ! apx = node_value('xsize ') |
|---|
| 451 | ! apy = node_value('ysize ') |
|---|
| 452 | ! if(apx.eq.zero) then |
|---|
| 453 | ! apx=maxaper(1) |
|---|
| 454 | ! endif |
|---|
| 455 | ! if(apy.eq.zero) then |
|---|
| 456 | ! apy=maxaper(3) |
|---|
| 457 | ! endif |
|---|
| 458 | ! call trcoll(1, apx, apy, turn, sum, part_id, last_turn, & |
|---|
| 459 | ! last_pos, last_orbit, track, ktrack,al_errors) |
|---|
| 460 | ! go to 500 |
|---|
| 461 | ! endif |
|---|
| 462 | !!---- Collimator with rectangular aperture. |
|---|
| 463 | ! if(code.eq.21) then |
|---|
| 464 | ! apx = node_value('xsize ') |
|---|
| 465 | ! apy = node_value('ysize ') |
|---|
| 466 | ! if(apx.eq.zero) then |
|---|
| 467 | ! apx=maxaper(1) |
|---|
| 468 | ! endif |
|---|
| 469 | ! if(apy.eq.zero) then |
|---|
| 470 | ! apy=maxaper(3) |
|---|
| 471 | ! endif |
|---|
| 472 | ! call trcoll(2, apx, apy, turn, sum, part_id, last_turn, & |
|---|
| 473 | ! last_pos, last_orbit, track, ktrack,al_errors) |
|---|
| 474 | ! go to 500 |
|---|
| 475 | ! endif |
|---|
| 476 | |
|---|
| 477 | |
|---|
| 478 | !---- Test aperture. ALL ELEMENTS BUT DRIFTS |
|---|
| 479 | aperflag = get_option('aperture ') .ne. 0 |
|---|
| 480 | if(aperflag) then |
|---|
| 481 | nn=name_len |
|---|
| 482 | call node_string('apertype ',aptype,nn) |
|---|
| 483 | call dzero(aperture,maxnaper) |
|---|
| 484 | call get_node_vector('aperture ',nn,aperture) |
|---|
| 485 | call dzero(offset,2) |
|---|
| 486 | call get_node_vector('aper_offset ',nn,offset) |
|---|
| 487 | |
|---|
| 488 | offx = offset(1) |
|---|
| 489 | offy = offset(2) |
|---|
| 490 | ! print *, " TYPE ",aptype & |
|---|
| 491 | ! "values x y lhc",aperture(1),aperture(2),aperture(3) |
|---|
| 492 | !------------ ellipse case ---------------------------------- |
|---|
| 493 | if(aptype.eq.'ellipse') then |
|---|
| 494 | apx = aperture(1) |
|---|
| 495 | apy = aperture(2) |
|---|
| 496 | call trcoll(1, apx, apy, turn, sum, part_id, last_turn, & |
|---|
| 497 | last_pos, last_orbit, track, ktrack,al_errors,offx,offy) |
|---|
| 498 | !------------ circle case ---------------------------------- |
|---|
| 499 | else if(aptype.eq.'circle') then |
|---|
| 500 | apx = aperture(1) |
|---|
| 501 | ! print *,"radius of circle in element",apx |
|---|
| 502 | if(apx.eq.zero) then |
|---|
| 503 | apx = maxaper(1) |
|---|
| 504 | ! print *,"radius of circle by default",apx |
|---|
| 505 | endif |
|---|
| 506 | apy = apx |
|---|
| 507 | ! print *,"circle, radius= ",apx |
|---|
| 508 | call trcoll(1, apx, apy, turn, sum, part_id, last_turn, & |
|---|
| 509 | last_pos, last_orbit, track,ktrack,al_errors,offx,offy) |
|---|
| 510 | !------------ rectangle case ---------------------------------- |
|---|
| 511 | else if(aptype.eq.'rectangle') then |
|---|
| 512 | apx = aperture(1) |
|---|
| 513 | apy = aperture(2) |
|---|
| 514 | call trcoll(2, apx, apy, turn, sum, part_id, last_turn, & |
|---|
| 515 | last_pos, last_orbit, track,ktrack,al_errors,offx,offy) |
|---|
| 516 | !-------------Racetrack type , added by Yipeng SUN 21-10-2008--- |
|---|
| 517 | else if(aptype.eq.'racetrack') then |
|---|
| 518 | apx = aperture(1) |
|---|
| 519 | apy = aperture(2) |
|---|
| 520 | apr = aperture(3) |
|---|
| 521 | call trcoll1(4, apx, apy, turn, sum, part_id, last_turn, & |
|---|
| 522 | last_pos,last_orbit,track,ktrack,al_errors,apr,offx,offy) |
|---|
| 523 | !------------ LHC screen case ---------------------------------- |
|---|
| 524 | else if(aptype.eq.'lhcscreen') then |
|---|
| 525 | ! print *, "LHC screen start, Xrect= ", |
|---|
| 526 | ! aperture(1)," Yrect= ",aperture(2)," Rcirc= ",aperture(3) |
|---|
| 527 | apx = aperture(3) |
|---|
| 528 | apy = aperture(3) |
|---|
| 529 | !JMJ! Making essential changes in AV's absence, 16/7/2003 |
|---|
| 530 | !JMJ! this tests whether the particle is outside the circumscribing |
|---|
| 531 | !JMJ! circle. |
|---|
| 532 | call trcoll(1, apx, apy, turn, sum, part_id, last_turn, & |
|---|
| 533 | last_pos, last_orbit, track,ktrack,al_errors,offx,offy) |
|---|
| 534 | !JMJ! |
|---|
| 535 | !JMJ! This tests whether particles are outside the space bounded by |
|---|
| 536 | !JMJ! two or four lines intersecting the circle. |
|---|
| 537 | !JMJ! The previous version checked that it was outside a rectangle, |
|---|
| 538 | !JMJ! one of whose dimensions was zero (in the LHC) so particles |
|---|
| 539 | !JMJ! were ALWAYS lost on the beam screen. |
|---|
| 540 | !JMJ! The new version of the test works on the understanding that |
|---|
| 541 | !JMJ! values of aperture(1) or aperture(2) greater than aperture(3) |
|---|
| 542 | !JMJ! have no meaning and that specifying a zero value is equivalent. |
|---|
| 543 | !JMJ! The most general aperture described by the present |
|---|
| 544 | !JMJ! implementation of LHCSCREEN is a rectangle with rounded |
|---|
| 545 | !JMJ! corners. |
|---|
| 546 | !JMJ! N.B. apx and apy already have the value aperture(3); the "if"s |
|---|
| 547 | !JMJ! ensure that they don't get set to zero. |
|---|
| 548 | if(aperture(1).gt.0.) apx = aperture(1) |
|---|
| 549 | if(aperture(2).gt.0.) apy = aperture(2) |
|---|
| 550 | call trcoll(2, apx, apy, turn, sum, part_id, last_turn, & |
|---|
| 551 | last_pos, last_orbit, track,ktrack,al_errors,offx,offy) |
|---|
| 552 | ! print *, "LHC screen end" |
|---|
| 553 | !------------ marguerite case ---------------------------------- |
|---|
| 554 | else if(aptype.eq.'marguerite') then |
|---|
| 555 | apx = aperture(1) |
|---|
| 556 | apy = aperture(2) |
|---|
| 557 | call trcoll(3, apx, apy, turn, sum, part_id, last_turn, & |
|---|
| 558 | last_pos, last_orbit, track,ktrack,al_errors,offx,offy) |
|---|
| 559 | !------------ rectellipse case ---------------------------------- |
|---|
| 560 | else if(aptype.eq.'rectellipse') then |
|---|
| 561 | !***** test ellipse |
|---|
| 562 | apx = aperture(3) |
|---|
| 563 | apy = aperture(4) |
|---|
| 564 | call trcoll(1, apx, apy, turn, sum, part_id, last_turn, & |
|---|
| 565 | last_pos, last_orbit, track,ktrack,al_errors,offx,offy) |
|---|
| 566 | !***** test rectangle |
|---|
| 567 | apx = aperture(1) |
|---|
| 568 | apy = aperture(2) |
|---|
| 569 | call trcoll(2, apx, apy, turn, sum, part_id, last_turn, & |
|---|
| 570 | last_pos, last_orbit, track,ktrack,al_errors,offx,offy) |
|---|
| 571 | ! print*, " test apertures" |
|---|
| 572 | ! print*, " apx=",apx, " apy=",apy," apxell=",apxell, |
|---|
| 573 | ! " apyell=",apyell |
|---|
| 574 | ! call trcoll(3, apx, apy, turn, sum, part_id, last_turn, & |
|---|
| 575 | ! last_pos, last_orbit, track,ktrack,al_errors) |
|---|
| 576 | endif |
|---|
| 577 | ! else |
|---|
| 578 | !!---- Check on collimator xsize/ysize even if user did not use 'aperture' |
|---|
| 579 | !!---- option in track command. This simulates roughly the old MAD-X |
|---|
| 580 | !!---- behaviour. |
|---|
| 581 | !!---- Collimator with elliptic aperture. |
|---|
| 582 | ! if(code.eq.20) then |
|---|
| 583 | ! apx = node_value('xsize ') |
|---|
| 584 | ! apy = node_value('ysize ') |
|---|
| 585 | ! if(apx.eq.zero) then |
|---|
| 586 | ! apx=maxaper(1) |
|---|
| 587 | ! endif |
|---|
| 588 | ! if(apy.eq.zero) then |
|---|
| 589 | ! apy=maxaper(3) |
|---|
| 590 | ! endif |
|---|
| 591 | ! call trcoll(1, apx, apy, turn, sum, part_id, last_turn, & |
|---|
| 592 | ! last_pos, last_orbit, track, ktrack,al_errors) |
|---|
| 593 | ! go to 500 |
|---|
| 594 | ! endif |
|---|
| 595 | !!---- Collimator with rectangular aperture. |
|---|
| 596 | ! if(code.eq.21) then |
|---|
| 597 | ! apx = node_value('xsize ') |
|---|
| 598 | ! apy = node_value('ysize ') |
|---|
| 599 | ! if(apx.eq.zero) then |
|---|
| 600 | ! apx=maxaper(1) |
|---|
| 601 | ! endif |
|---|
| 602 | ! if(apy.eq.zero) then |
|---|
| 603 | ! apy=maxaper(3) |
|---|
| 604 | ! endif |
|---|
| 605 | ! call trcoll(2, apx, apy, turn, sum, part_id, last_turn, & |
|---|
| 606 | ! last_pos, last_orbit, track, ktrack,al_errors) |
|---|
| 607 | ! go to 500 |
|---|
| 608 | ! endif |
|---|
| 609 | endif |
|---|
| 610 | !---- END OF IF(APERFLAG) ---------------- |
|---|
| 611 | |
|---|
| 612 | |
|---|
| 613 | !-- switch on element type BUT DRIFT, COLLIMATORS, BEAM_BEAM / 13.03.03 |
|---|
| 614 | !-- 500 has been specified at the relevant places in go to below |
|---|
| 615 | !-- code =1 for drift, treated above, go to 500 directly |
|---|
| 616 | ! print *," CODE ",code |
|---|
| 617 | go to ( 500, 20, 30, 40, 50, 60, 70, 80, 90, 100, & |
|---|
| 618 | 110, 120, 130, 140, 150, 160, 170, 180, 190, 500, & |
|---|
| 619 | 500, 500, 230, 240, 250, 260, 270, 280, 290, 300, 310, 320, & |
|---|
| 620 | ! 330, 500, 350, 360, 370,500,500,400,410,500, 500, 500, 500), code |
|---|
| 621 | ! Use this line to enable non-linear thin lens |
|---|
| 622 | ! 330, 500, 350, 360, 370,500,500,400,410,420,500, 500, 500, 500), code |
|---|
| 623 | ! Enable non-linear thin lens and RF-Multipole |
|---|
| 624 | 330, 500, 350, 360, 370,500,500,400,410,420,430, 500, 500, 500), code |
|---|
| 625 | ! |
|---|
| 626 | !---- Make sure that nothing is execute if element is not known |
|---|
| 627 | go to 500 |
|---|
| 628 | ! |
|---|
| 629 | !---- Bending magnet. OBSOLETE, to be kept for go to |
|---|
| 630 | 20 continue ! RBEND |
|---|
| 631 | 30 continue ! SBEND |
|---|
| 632 | go to 500 |
|---|
| 633 | !---- Arbitrary matrix. OBSOLETE, to be kept for go to |
|---|
| 634 | 40 continue |
|---|
| 635 | call tmarb(.false.,.false.,craporb,fmap,ek,re,te) |
|---|
| 636 | call tttrak(ek,re,track,ktrack) |
|---|
| 637 | go to 500 |
|---|
| 638 | !---- Quadrupole. OBSOLETE, to be kept for go to |
|---|
| 639 | !---- AL: not so fast, if here it's a thick quadrupole |
|---|
| 640 | 50 continue |
|---|
| 641 | call tttquad(track,ktrack) |
|---|
| 642 | go to 500 |
|---|
| 643 | !---- Sextupole. OBSOLETE, to be kept for go to |
|---|
| 644 | 60 continue |
|---|
| 645 | go to 500 |
|---|
| 646 | !---- Octupole. OBSOLETE, to be kept for go to |
|---|
| 647 | 70 continue |
|---|
| 648 | go to 500 |
|---|
| 649 | !---- Monitors, beam instrument., MONITOR, HMONITOR, VMONITOR |
|---|
| 650 | 170 continue |
|---|
| 651 | 180 continue |
|---|
| 652 | 190 continue |
|---|
| 653 | go to 500 |
|---|
| 654 | !---- Multipole |
|---|
| 655 | 80 continue |
|---|
| 656 | call ttmult(track,ktrack,dxt,dyt,turn) |
|---|
| 657 | go to 500 |
|---|
| 658 | !---- Solenoid. |
|---|
| 659 | 90 continue |
|---|
| 660 | call trsol(track, ktrack) |
|---|
| 661 | go to 500 |
|---|
| 662 | !---- RF cavity. |
|---|
| 663 | 100 continue |
|---|
| 664 | call ttrf(track,ktrack) |
|---|
| 665 | go to 500 |
|---|
| 666 | !---- Electrostatic separator. |
|---|
| 667 | 110 continue |
|---|
| 668 | call ttsep(el, track, ktrack) |
|---|
| 669 | go to 500 |
|---|
| 670 | !---- Rotation around s-axis. |
|---|
| 671 | 120 continue |
|---|
| 672 | call ttsrot(track, ktrack) |
|---|
| 673 | go to 500 |
|---|
| 674 | !---- Rotation around y-axis. |
|---|
| 675 | 130 continue |
|---|
| 676 | ! call ttyrot(track, ktrack) |
|---|
| 677 | go to 500 |
|---|
| 678 | !---- Correctors. |
|---|
| 679 | 140 continue |
|---|
| 680 | 150 continue |
|---|
| 681 | 160 continue |
|---|
| 682 | call ttcorr(el, track, ktrack, turn) |
|---|
| 683 | go to 500 |
|---|
| 684 | !---- ECollimator, RCollimator, BeamBeam, Lump. ?? |
|---|
| 685 | 230 continue |
|---|
| 686 | go to 500 |
|---|
| 687 | !---- Instrument |
|---|
| 688 | 240 continue |
|---|
| 689 | go to 500 |
|---|
| 690 | !---- Marker. |
|---|
| 691 | 250 continue |
|---|
| 692 | go to 500 |
|---|
| 693 | !---- General bend (dipole, quadrupole, and skew quadrupole). |
|---|
| 694 | 260 continue |
|---|
| 695 | go to 500 |
|---|
| 696 | !---- LCAV cavity. |
|---|
| 697 | 270 continue |
|---|
| 698 | ! call ttlcav(el, track, ktrack) |
|---|
| 699 | go to 500 |
|---|
| 700 | !---- Reserved. (PROFILE,WIRE,SLMONITOR,BLMONITOR,IMONITOR) |
|---|
| 701 | 280 continue |
|---|
| 702 | 290 continue |
|---|
| 703 | 300 continue |
|---|
| 704 | 310 continue |
|---|
| 705 | 320 continue |
|---|
| 706 | go to 500 |
|---|
| 707 | !---- Dipole edge |
|---|
| 708 | 330 continue |
|---|
| 709 | call ttdpdg(track,ktrack) |
|---|
| 710 | go to 500 |
|---|
| 711 | !---- Changeref ??? |
|---|
| 712 | 350 continue |
|---|
| 713 | go to 500 |
|---|
| 714 | !---- Translation |
|---|
| 715 | 360 continue |
|---|
| 716 | go to 500 |
|---|
| 717 | !---- Crab cavity. |
|---|
| 718 | 370 continue |
|---|
| 719 | call ttcrabrf(track,ktrack,turn) |
|---|
| 720 | go to 500 |
|---|
| 721 | !---- h ac dipole. |
|---|
| 722 | 400 continue |
|---|
| 723 | call tthacdip(track,ktrack,turn) |
|---|
| 724 | go to 500 |
|---|
| 725 | !---- v ac dipole. |
|---|
| 726 | 410 continue |
|---|
| 727 | call ttvacdip(track,ktrack,turn) |
|---|
| 728 | go to 500 |
|---|
| 729 | !---- nonlinear elliptical lens |
|---|
| 730 | 420 continue |
|---|
| 731 | call ttnllens(track,ktrack) |
|---|
| 732 | go to 500 |
|---|
| 733 | !---- rf multipoles |
|---|
| 734 | 430 continue |
|---|
| 735 | call ttrfmult(track,ktrack,turn) |
|---|
| 736 | go to 500 |
|---|
| 737 | |
|---|
| 738 | !---- Solenoid. |
|---|
| 739 | |
|---|
| 740 | |
|---|
| 741 | 500 continue |
|---|
| 742 | |
|---|
| 743 | !---- Rotate trajectory at exit |
|---|
| 744 | if (theta .ne. zero) then |
|---|
| 745 | do jtrk = 1,ktrack |
|---|
| 746 | tmp = track(1,jtrk) |
|---|
| 747 | track(1,jtrk) = ct * tmp - st * track(3,jtrk) |
|---|
| 748 | track(3,jtrk) = ct * track(3,jtrk) + st * tmp |
|---|
| 749 | tmp = track(2,jtrk) |
|---|
| 750 | track(2,jtrk) = ct * tmp - st * track(4,jtrk) |
|---|
| 751 | track(4,jtrk) = ct * track(4,jtrk) + st * tmp |
|---|
| 752 | enddo |
|---|
| 753 | endif |
|---|
| 754 | |
|---|
| 755 | !---- Accumulate length. |
|---|
| 756 | 502 sum = sum + el |
|---|
| 757 | return |
|---|
| 758 | end subroutine ttmap |
|---|
| 759 | |
|---|
| 760 | subroutine ttmult(track, ktrack,dxt,dyt,turn) |
|---|
| 761 | |
|---|
| 762 | use twtrrfi |
|---|
| 763 | use name_lenfi |
|---|
| 764 | use trackfi |
|---|
| 765 | implicit none |
|---|
| 766 | |
|---|
| 767 | !----------------------------------------------------------------------* |
|---|
| 768 | ! Purpose: * |
|---|
| 769 | ! Track particle through a general thin multipole. * |
|---|
| 770 | ! Input/output: * |
|---|
| 771 | ! TRACK(6,*)(double) Track coordinates: (X, PX, Y, PY, T, PT). * |
|---|
| 772 | ! KTRACK (integer) Number of surviving tracks. * |
|---|
| 773 | ! dxt (double) local buffer * |
|---|
| 774 | ! dyt (double) local buffer * |
|---|
| 775 | !----------------------------------------------------------------------* |
|---|
| 776 | logical first |
|---|
| 777 | integer iord,jtrk,nd,nord,ktrack,j,n_ferr,nn,ns,node_fd_errors, & |
|---|
| 778 | get_option,turn,noisemax,nn1,in |
|---|
| 779 | double precision const,curv,dbi,dbr,dipi,dipr,dx,dy,elrad, & |
|---|
| 780 | pt,px,py,rfac,rpt1,rpt2,rpx1,rpx2,rpy1,rpy2, & |
|---|
| 781 | f_errors(0:maxferr),field(2,0:maxmul),vals(2,0:maxmul), & |
|---|
| 782 | ordinv(maxmul),track(6,*),dxt(*),dyt(*),normal(0:maxmul), & |
|---|
| 783 | skew(0:maxmul),bvk,node_value,zero,one,two,three,half,ttt, & |
|---|
| 784 | npeak(100), nlag(100), ntune(100), temp,noise |
|---|
| 785 | |
|---|
| 786 | parameter(zero=0d0,one=1d0,two=2d0,three=3d0,half=5d-1) |
|---|
| 787 | character(name_len) aptype |
|---|
| 788 | |
|---|
| 789 | save first,ordinv |
|---|
| 790 | data first / .true. / |
|---|
| 791 | |
|---|
| 792 | |
|---|
| 793 | !---- Precompute reciprocals of orders. |
|---|
| 794 | if (first) then |
|---|
| 795 | do iord = 1, maxmul |
|---|
| 796 | ordinv(iord) = one / float(iord) |
|---|
| 797 | enddo |
|---|
| 798 | first = .false. |
|---|
| 799 | endif |
|---|
| 800 | call dzero(f_errors, maxferr+1) |
|---|
| 801 | n_ferr = node_fd_errors(f_errors) |
|---|
| 802 | bvk = node_value('other_bv ') |
|---|
| 803 | !---- Multipole length for radiation. |
|---|
| 804 | elrad = node_value('lrad ') |
|---|
| 805 | noise = node_value('noise ') |
|---|
| 806 | !---- Multipole components. |
|---|
| 807 | call dzero(normal,maxmul+1) |
|---|
| 808 | call dzero(skew,maxmul+1) |
|---|
| 809 | call get_node_vector('knl ',nn,normal) |
|---|
| 810 | call get_node_vector('ksl ',ns,skew) |
|---|
| 811 | nd = 2 * max(nn, ns) |
|---|
| 812 | call dzero(vals,2*(maxmul+1)) |
|---|
| 813 | |
|---|
| 814 | if(noise .eq. 1) then |
|---|
| 815 | nn1=name_len |
|---|
| 816 | noisemax = node_value('noisemax ') |
|---|
| 817 | call dzero(npeak,noisemax) |
|---|
| 818 | call dzero(ntune,noisemax) |
|---|
| 819 | call dzero(nlag,noisemax) |
|---|
| 820 | call get_node_vector('npeak ',nn1,npeak) |
|---|
| 821 | call get_node_vector('ntune ',nn1,ntune) |
|---|
| 822 | call get_node_vector('nlag ',nn1,nlag) |
|---|
| 823 | |
|---|
| 824 | temp = 0 |
|---|
| 825 | do in = 1, noisemax |
|---|
| 826 | temp = temp + npeak(in) * sin(nlag(in) + ntune(in) * turn) |
|---|
| 827 | enddo |
|---|
| 828 | |
|---|
| 829 | |
|---|
| 830 | ! temp = npeak * sin(nlag + ntune * turn) |
|---|
| 831 | do iord = 0, nn |
|---|
| 832 | vals(1,iord) = normal(iord) * (1+temp) |
|---|
| 833 | enddo |
|---|
| 834 | do iord = 0, ns |
|---|
| 835 | vals(2,iord) = skew(iord) * (1+temp) |
|---|
| 836 | enddo |
|---|
| 837 | else |
|---|
| 838 | do iord = 0, nn |
|---|
| 839 | vals(1,iord) = normal(iord) |
|---|
| 840 | enddo |
|---|
| 841 | do iord = 0, ns |
|---|
| 842 | vals(2,iord) = skew(iord) |
|---|
| 843 | enddo |
|---|
| 844 | endif |
|---|
| 845 | |
|---|
| 846 | ! do iord = 0, nn |
|---|
| 847 | ! vals(1,iord) = normal(iord) |
|---|
| 848 | ! enddo |
|---|
| 849 | ! do iord = 0, ns |
|---|
| 850 | ! vals(2,iord) = skew(iord) |
|---|
| 851 | ! enddo |
|---|
| 852 | !---- Field error vals. |
|---|
| 853 | call dzero(field,2*(maxmul+1)) |
|---|
| 854 | if (n_ferr .gt. 0) then |
|---|
| 855 | call dcopy(f_errors,field,n_ferr) |
|---|
| 856 | endif |
|---|
| 857 | !-----added FrankS, 10-12-2008 |
|---|
| 858 | nd = 2 * max(nn, ns, n_ferr/2-1) |
|---|
| 859 | !---- Dipole error. |
|---|
| 860 | ! dbr = bvk * field(1,0) / (one + deltas) |
|---|
| 861 | ! dbi = bvk * field(2,0) / (one + deltas) |
|---|
| 862 | dbr = bvk * field(1,0) |
|---|
| 863 | dbi = bvk * field(2,0) |
|---|
| 864 | !---- Nominal dipole strength. |
|---|
| 865 | ! dipr = bvk * vals(1,0) / (one + deltas) |
|---|
| 866 | ! dipi = bvk * vals(2,0) / (one + deltas) |
|---|
| 867 | dipr = bvk * vals(1,0) |
|---|
| 868 | dipi = bvk * vals(2,0) |
|---|
| 869 | !---- Other components and errors. |
|---|
| 870 | nord = 0 |
|---|
| 871 | do iord = 1, nd/2 |
|---|
| 872 | do j = 1, 2 |
|---|
| 873 | ! field(j,iord) = bvk * (vals(j,iord) + field(j,iord)) & |
|---|
| 874 | ! / (one + deltas) |
|---|
| 875 | field(j,iord) = bvk * (vals(j,iord) + field(j,iord)) |
|---|
| 876 | if (field(j,iord) .ne. zero) nord = iord |
|---|
| 877 | enddo |
|---|
| 878 | enddo |
|---|
| 879 | !---- Pure dipole: only quadrupole kicks according to lrad. |
|---|
| 880 | if (nord .eq. 0) then |
|---|
| 881 | do jtrk = 1,ktrack |
|---|
| 882 | dxt(jtrk) = zero |
|---|
| 883 | dyt(jtrk) = zero |
|---|
| 884 | enddo |
|---|
| 885 | !----------- introduction of dipole focusing |
|---|
| 886 | if(elrad.gt.zero.and.get_option('thin_foc ').eq.1) then |
|---|
| 887 | do jtrk = 1,ktrack |
|---|
| 888 | dxt(jtrk) = dipr*dipr*track(1,jtrk)/elrad |
|---|
| 889 | dyt(jtrk) = dipi*dipi*track(3,jtrk)/elrad |
|---|
| 890 | enddo |
|---|
| 891 | endif |
|---|
| 892 | !---- Accumulate multipole kick from highest multipole to quadrupole. |
|---|
| 893 | else |
|---|
| 894 | do jtrk = 1,ktrack |
|---|
| 895 | dxt(jtrk) = & |
|---|
| 896 | field(1,nord)*track(1,jtrk) - field(2,nord)*track(3,jtrk) |
|---|
| 897 | dyt(jtrk) = & |
|---|
| 898 | field(1,nord)*track(3,jtrk) + field(2,nord)*track(1,jtrk) |
|---|
| 899 | enddo |
|---|
| 900 | |
|---|
| 901 | do iord = nord - 1, 1, -1 |
|---|
| 902 | do jtrk = 1,ktrack |
|---|
| 903 | dx = dxt(jtrk)*ordinv(iord+1) + field(1,iord) |
|---|
| 904 | dy = dyt(jtrk)*ordinv(iord+1) + field(2,iord) |
|---|
| 905 | dxt(jtrk) = dx*track(1,jtrk) - dy*track(3,jtrk) |
|---|
| 906 | dyt(jtrk) = dx*track(3,jtrk) + dy*track(1,jtrk) |
|---|
| 907 | enddo |
|---|
| 908 | enddo |
|---|
| 909 | ! do jtrk = 1,ktrack |
|---|
| 910 | ! dxt(jtrk) = dxt(jtrk) / (one + deltas) |
|---|
| 911 | ! dyt(jtrk) = dyt(jtrk) / (one + deltas) |
|---|
| 912 | ! enddo |
|---|
| 913 | if(elrad.gt.zero.and.get_option('thin_foc ').eq.1) then |
|---|
| 914 | do jtrk = 1,ktrack |
|---|
| 915 | dxt(jtrk) = dxt(jtrk) + dipr*dipr*track(1,jtrk)/elrad |
|---|
| 916 | dyt(jtrk) = dyt(jtrk) + dipi*dipi*track(3,jtrk)/elrad |
|---|
| 917 | enddo |
|---|
| 918 | endif |
|---|
| 919 | endif |
|---|
| 920 | |
|---|
| 921 | !---- Radiation loss at entrance. |
|---|
| 922 | if (dorad .and. elrad .ne. 0) then |
|---|
| 923 | const = arad * gammas**3 / three |
|---|
| 924 | |
|---|
| 925 | !---- Full damping. |
|---|
| 926 | if (dodamp) then |
|---|
| 927 | do jtrk = 1,ktrack |
|---|
| 928 | curv = sqrt((dipr + dxt(jtrk))**2 + & |
|---|
| 929 | (dipi + dyt(jtrk))**2) / elrad |
|---|
| 930 | |
|---|
| 931 | if (dorand) then |
|---|
| 932 | call trphot(elrad,curv,rfac,deltas) |
|---|
| 933 | else |
|---|
| 934 | rfac = const * curv**2 * elrad |
|---|
| 935 | endif |
|---|
| 936 | |
|---|
| 937 | px = track(2,jtrk) |
|---|
| 938 | py = track(4,jtrk) |
|---|
| 939 | pt = track(6,jtrk) |
|---|
| 940 | track(2,jtrk) = px - rfac * (one + pt) * px |
|---|
| 941 | track(4,jtrk) = py - rfac * (one + pt) * py |
|---|
| 942 | track(6,jtrk) = pt - rfac * (one + pt) ** 2 |
|---|
| 943 | enddo |
|---|
| 944 | |
|---|
| 945 | !---- Energy loss like for closed orbit. |
|---|
| 946 | else |
|---|
| 947 | |
|---|
| 948 | !---- Store energy loss on closed orbit. |
|---|
| 949 | rfac = const * ((dipr + dxt(1))**2 + (dipi + dyt(1))**2) |
|---|
| 950 | rpx1 = rfac * (one + track(6,1)) * track(2,1) |
|---|
| 951 | rpy1 = rfac * (one + track(6,1)) * track(4,1) |
|---|
| 952 | rpt1 = rfac * (one + track(6,1)) ** 2 |
|---|
| 953 | |
|---|
| 954 | do jtrk = 1,ktrack |
|---|
| 955 | track(2,jtrk) = track(2,jtrk) - rpx1 |
|---|
| 956 | track(4,jtrk) = track(4,jtrk) - rpy1 |
|---|
| 957 | track(6,jtrk) = track(6,jtrk) - rpt1 |
|---|
| 958 | enddo |
|---|
| 959 | |
|---|
| 960 | endif |
|---|
| 961 | endif |
|---|
| 962 | |
|---|
| 963 | !---- Apply multipole effect including dipole. |
|---|
| 964 | do jtrk = 1,ktrack |
|---|
| 965 | ! Added for correct Ripken implementation of formulae |
|---|
| 966 | ttt = sqrt(one+two*track(6,jtrk)*bet0i+track(6,jtrk)**2) |
|---|
| 967 | ! track(2,jtrk) = track(2,jtrk) - & |
|---|
| 968 | ! (dbr + dxt(jtrk) - dipr * (deltas + beti*track(6,jtrk))) |
|---|
| 969 | ! track(4,jtrk) = track(4,jtrk) + & |
|---|
| 970 | ! (dbi + dyt(jtrk) - dipi * (deltas + beti*track(6,jtrk))) |
|---|
| 971 | ! track(5,jtrk) = track(5,jtrk) & |
|---|
| 972 | ! - (dipr*track(1,jtrk) - dipi*track(3,jtrk)) * beti |
|---|
| 973 | track(2,jtrk) = track(2,jtrk) - & |
|---|
| 974 | (dbr + dxt(jtrk) - dipr * (ttt - one)) |
|---|
| 975 | track(4,jtrk) = track(4,jtrk) + & |
|---|
| 976 | (dbi + dyt(jtrk) - dipi * (ttt - one)) |
|---|
| 977 | track(5,jtrk) = track(5,jtrk) - & |
|---|
| 978 | (dipr*track(1,jtrk) - dipi*track(3,jtrk)) * & |
|---|
| 979 | ((one + bet0*track(6,jtrk))/ttt)*bet0i |
|---|
| 980 | enddo |
|---|
| 981 | |
|---|
| 982 | !---- Radiation loss at exit. |
|---|
| 983 | if (dorad .and. elrad .ne. 0) then |
|---|
| 984 | |
|---|
| 985 | !---- Full damping. |
|---|
| 986 | if (dodamp) then |
|---|
| 987 | do jtrk = 1,ktrack |
|---|
| 988 | curv = sqrt((dipr + dxt(jtrk))**2 + & |
|---|
| 989 | (dipi + dyt(jtrk))**2) / elrad |
|---|
| 990 | |
|---|
| 991 | if (dorand) then |
|---|
| 992 | call trphot(elrad,curv,rfac,deltas) |
|---|
| 993 | else |
|---|
| 994 | rfac = const * curv**2 * elrad |
|---|
| 995 | endif |
|---|
| 996 | |
|---|
| 997 | px = track(2,jtrk) |
|---|
| 998 | py = track(4,jtrk) |
|---|
| 999 | pt = track(6,jtrk) |
|---|
| 1000 | track(2,jtrk) = px - rfac * (one + pt) * px |
|---|
| 1001 | track(4,jtrk) = py - rfac * (one + pt) * py |
|---|
| 1002 | track(6,jtrk) = pt - rfac * (one + pt) ** 2 |
|---|
| 1003 | enddo |
|---|
| 1004 | |
|---|
| 1005 | !---- Energy loss like for closed orbit. |
|---|
| 1006 | else |
|---|
| 1007 | |
|---|
| 1008 | !---- Store energy loss on closed orbit. |
|---|
| 1009 | rfac = const * ((dipr + dxt(1))**2 + (dipi + dyt(1))**2) |
|---|
| 1010 | rpx2 = rfac * (one + track(6,1)) * track(2,1) |
|---|
| 1011 | rpy2 = rfac * (one + track(6,1)) * track(4,1) |
|---|
| 1012 | rpt2 = rfac * (one + track(6,1)) ** 2 |
|---|
| 1013 | |
|---|
| 1014 | do jtrk = 1,ktrack |
|---|
| 1015 | track(2,jtrk) = track(2,jtrk) - rpx2 |
|---|
| 1016 | track(4,jtrk) = track(4,jtrk) - rpy2 |
|---|
| 1017 | track(6,jtrk) = track(6,jtrk) - rpt2 |
|---|
| 1018 | enddo |
|---|
| 1019 | endif |
|---|
| 1020 | endif |
|---|
| 1021 | |
|---|
| 1022 | end subroutine ttmult |
|---|
| 1023 | |
|---|
| 1024 | subroutine ttsrot(track,ktrack) |
|---|
| 1025 | |
|---|
| 1026 | use trackfi |
|---|
| 1027 | implicit none |
|---|
| 1028 | |
|---|
| 1029 | !----------------------------------------------------------------------* |
|---|
| 1030 | ! Purpose: * |
|---|
| 1031 | ! Coordinate change due to rotation about s axis * |
|---|
| 1032 | ! Input/output: * |
|---|
| 1033 | ! TRACK(6,*)(double) Track coordinates: (X, PX, Y, PY, T, PT). * |
|---|
| 1034 | ! KTRACK (integer) number of surviving tracks. tmsrot * |
|---|
| 1035 | !----------------------------------------------------------------------* |
|---|
| 1036 | integer itrack,ktrack,j |
|---|
| 1037 | double precision track(6,*),psi,node_value,ct,st,trb(4) |
|---|
| 1038 | psi = node_value('angle ') |
|---|
| 1039 | ct = cos(psi) |
|---|
| 1040 | st = -sin(psi) |
|---|
| 1041 | do itrack = 1, ktrack |
|---|
| 1042 | do j = 1,4 |
|---|
| 1043 | trb(j)=track(j,itrack) |
|---|
| 1044 | enddo |
|---|
| 1045 | track(1,itrack) = trb(1)*ct-trb(3)*st |
|---|
| 1046 | track(2,itrack) = trb(2)*ct-trb(4)*st |
|---|
| 1047 | track(3,itrack) = trb(1)*st+trb(3)*ct |
|---|
| 1048 | track(4,itrack) = trb(2)*st+trb(4)*ct |
|---|
| 1049 | enddo |
|---|
| 1050 | end subroutine ttsrot |
|---|
| 1051 | |
|---|
| 1052 | |
|---|
| 1053 | subroutine ttyrot(track,ktrack) |
|---|
| 1054 | |
|---|
| 1055 | use trackfi |
|---|
| 1056 | implicit none |
|---|
| 1057 | |
|---|
| 1058 | !----------------------------------------------------------------------* |
|---|
| 1059 | ! Purpose: * |
|---|
| 1060 | ! Coordinate change due to rotation about y axis * |
|---|
| 1061 | ! Input/output: * |
|---|
| 1062 | ! TRACK(6,*)(double) Track coordinates: (X, PX, Y, PY, T, PT). * |
|---|
| 1063 | ! KTRACK (integer) number of surviving tracks. tmsrot * |
|---|
| 1064 | !----------------------------------------------------------------------* |
|---|
| 1065 | integer itrack,ktrack,j |
|---|
| 1066 | double precision track(6,*),theta,node_value,ct,st,trb(6) |
|---|
| 1067 | theta = node_value('angle ') |
|---|
| 1068 | ct = cos(theta) |
|---|
| 1069 | st = -sin(theta) |
|---|
| 1070 | do itrack = 1, ktrack |
|---|
| 1071 | do j = 1,6 |
|---|
| 1072 | trb(j)=track(j,itrack) |
|---|
| 1073 | enddo |
|---|
| 1074 | track(1,itrack) = trb(1)*ct-trb(5)*st |
|---|
| 1075 | track(2,itrack) = trb(2)*ct-trb(6)*st |
|---|
| 1076 | track(5,itrack) = trb(1)*st+trb(5)*ct |
|---|
| 1077 | track(6,itrack) = trb(2)*st+trb(6)*ct |
|---|
| 1078 | enddo |
|---|
| 1079 | end subroutine ttyrot |
|---|
| 1080 | |
|---|
| 1081 | subroutine ttdrf(el,track,ktrack) |
|---|
| 1082 | |
|---|
| 1083 | use trackfi |
|---|
| 1084 | implicit none |
|---|
| 1085 | |
|---|
| 1086 | !----------------------------------------------------------------------* |
|---|
| 1087 | ! Purpose: * |
|---|
| 1088 | ! Track a set of particle through a drift space. * |
|---|
| 1089 | ! Input/output: * |
|---|
| 1090 | ! TRACK(6,*)(double) Track coordinates: (X, PX, Y, PY, T, PT). * |
|---|
| 1091 | ! KTRACK (integer) number of surviving tracks. * |
|---|
| 1092 | ! Output: * |
|---|
| 1093 | ! EL (double) Length of quadrupole. * |
|---|
| 1094 | !----------------------------------------------------------------------* |
|---|
| 1095 | integer itrack,ktrack |
|---|
| 1096 | double precision el,pt,px,py,track(6,*),ttt,one,two |
|---|
| 1097 | parameter(one=1d0,two=2d0) |
|---|
| 1098 | |
|---|
| 1099 | ! picked from trturn in madx.ss |
|---|
| 1100 | do itrack = 1, ktrack |
|---|
| 1101 | px = track(2,itrack) |
|---|
| 1102 | py = track(4,itrack) |
|---|
| 1103 | pt = track(6,itrack) |
|---|
| 1104 | ttt = el/sqrt(one+two*pt*bet0i+pt**2 - px**2 - py**2) |
|---|
| 1105 | track(1,itrack) = track(1,itrack) + ttt*px |
|---|
| 1106 | track(3,itrack) = track(3,itrack) + ttt*py |
|---|
| 1107 | ! track(5,itrack) = track(5,itrack) & |
|---|
| 1108 | ! + el*(beti + pt * dtbyds) - (beti+pt)*ttt |
|---|
| 1109 | !---- AK 20060413 |
|---|
| 1110 | !---- Ripken DESY-95-189 p.36 |
|---|
| 1111 | track(5,itrack) = track(5,itrack) & |
|---|
| 1112 | + bet0i*(el - (one + bet0*pt) * ttt) |
|---|
| 1113 | enddo |
|---|
| 1114 | end subroutine ttdrf |
|---|
| 1115 | subroutine ttrf(track,ktrack) |
|---|
| 1116 | |
|---|
| 1117 | implicit none |
|---|
| 1118 | |
|---|
| 1119 | !----------------------------------------------------------------------* |
|---|
| 1120 | ! Purpose: * |
|---|
| 1121 | ! Track a set of trajectories through a thin cavity (zero length). * |
|---|
| 1122 | ! The cavity is sandwiched between two drift spaces of half length. * |
|---|
| 1123 | ! Input/output: * |
|---|
| 1124 | ! TRACK(6,*)(double) Track coordinates: (X, PX, Y, PY, T, PT). * |
|---|
| 1125 | ! KTRACK (integer) number of surviving tracks. * |
|---|
| 1126 | ! Output: * |
|---|
| 1127 | ! EL (double) Length of quadrupole. * |
|---|
| 1128 | !----------------------------------------------------------------------* |
|---|
| 1129 | ! Modified: 06-JAN-1999, T. Raubenheimer (SLAC) * |
|---|
| 1130 | ! Modified to allow wakefield tracking however no modification to * |
|---|
| 1131 | ! the logic and the nominal energy is not updated -- see routine * |
|---|
| 1132 | ! TTLCAV to change the nominal energy * |
|---|
| 1133 | !----------------------------------------------------------------------* |
|---|
| 1134 | integer itrack,ktrack |
|---|
| 1135 | double precision omega,phirf,pt,rff,rfl,rfv,track(6,*),clight, & |
|---|
| 1136 | twopi,vrf,pc0,get_variable,node_value,get_value,one,two,half,bvk, & |
|---|
| 1137 | ten3m,ten6p |
|---|
| 1138 | ! double precision px,py,ttt,beti,el1 |
|---|
| 1139 | parameter(one=1d0,two=2d0,half=5d-1,ten3m=1d-3,ten6p=1d6) |
|---|
| 1140 | |
|---|
| 1141 | !---- Initialize |
|---|
| 1142 | clight=get_variable('clight ') |
|---|
| 1143 | twopi=get_variable('twopi ') |
|---|
| 1144 | |
|---|
| 1145 | !---- BV flag |
|---|
| 1146 | bvk = node_value('other_bv ') |
|---|
| 1147 | |
|---|
| 1148 | !---- Fetch data. |
|---|
| 1149 | ! el = node_value('l ') |
|---|
| 1150 | ! el1 = node_value('l ') |
|---|
| 1151 | rfv = bvk * node_value('volt ') |
|---|
| 1152 | rff = node_value('freq ') |
|---|
| 1153 | rfl = node_value('lag ') |
|---|
| 1154 | ! deltap = get_value('probe ','deltap ') |
|---|
| 1155 | ! deltas = get_variable('track_deltap ') |
|---|
| 1156 | ! pc = get_value('probe ','pc ') |
|---|
| 1157 | pc0 = get_value('beam ','pc ') |
|---|
| 1158 | ! betas = get_value('probe ','beta ') |
|---|
| 1159 | ! gammas= get_value('probe ','gamma ') |
|---|
| 1160 | ! dtbyds = get_value('probe ','dtbyds ') |
|---|
| 1161 | |
|---|
| 1162 | ! print *,"RF cav. volt=",rfv, " freq.",rff |
|---|
| 1163 | |
|---|
| 1164 | !*---- Get the longitudinal wakefield filename (parameter #17). |
|---|
| 1165 | ! if (iq(lcelm+melen+15*mcsiz-2) .eq. 61) then |
|---|
| 1166 | ! lstr = iq(lcelm+melen+15*mcsiz) |
|---|
| 1167 | ! call uhtoc(iq(lq(lcelm-17)+1), mcwrd, lfile, 80) |
|---|
| 1168 | ! else |
|---|
| 1169 | ! lfile = ' ' |
|---|
| 1170 | ! endif |
|---|
| 1171 | |
|---|
| 1172 | !*---- Get the transverse wakefield filename (parameter #18). |
|---|
| 1173 | ! if (iq(lcelm+melen+16*mcsiz-2) .eq. 61) then |
|---|
| 1174 | ! lstr = iq(lcelm+melen+16*mcsiz) |
|---|
| 1175 | ! call uhtoc(iq(lq(lcelm-18)+1), mcwrd, tfile, 80) |
|---|
| 1176 | ! else |
|---|
| 1177 | ! tfile = ' ' |
|---|
| 1178 | ! endif |
|---|
| 1179 | |
|---|
| 1180 | !*---- If there are wakefields split the cavity. |
|---|
| 1181 | ! if (lfile .ne. ' ' .or. tfile .ne. ' ') then |
|---|
| 1182 | ! el1 = el / two |
|---|
| 1183 | ! rfv = rfv / two |
|---|
| 1184 | ! lwake = .true. |
|---|
| 1185 | ! else |
|---|
| 1186 | ! el1 = el |
|---|
| 1187 | ! lwake = .false. |
|---|
| 1188 | ! endif |
|---|
| 1189 | !---- Set up. |
|---|
| 1190 | omega = rff * (ten6p * twopi / clight) |
|---|
| 1191 | ! vrf = rfv * ten3m / (pc * (one + deltas)) |
|---|
| 1192 | vrf = rfv * ten3m / pc0 |
|---|
| 1193 | phirf = rfl * twopi |
|---|
| 1194 | ! dl = el * half |
|---|
| 1195 | ! bi2gi2 = one / (betas * gammas) ** 2 |
|---|
| 1196 | !---- Loop for all particles. |
|---|
| 1197 | ! do itrack = 1, ktrack |
|---|
| 1198 | ! pt = track(6,itrack) |
|---|
| 1199 | ! pt = pt + vrf * sin(phirf - omega * track(5,itrack)) |
|---|
| 1200 | ! track(6,itrack) = pt |
|---|
| 1201 | !! print *," pt ",pt, " track(5,itrack)",track(5,itrack) |
|---|
| 1202 | ! enddo |
|---|
| 1203 | do itrack = 1, ktrack |
|---|
| 1204 | pt = track(6,itrack) |
|---|
| 1205 | pt = pt + vrf * sin(phirf - omega * track(5,itrack)) |
|---|
| 1206 | track(6,itrack) = pt |
|---|
| 1207 | ! print *," pt ",pt, " track(5,itrack)",track(5,itrack) |
|---|
| 1208 | enddo |
|---|
| 1209 | |
|---|
| 1210 | !*---- If there were wakefields, track the wakes and then the 2nd half |
|---|
| 1211 | !* of the cavity. |
|---|
| 1212 | ! if (lwake) then |
|---|
| 1213 | ! call ttwake(two*el1, nbin, binmax, lfile, tfile, ener1, track, |
|---|
| 1214 | ! + ktrack) |
|---|
| 1215 | ! |
|---|
| 1216 | !*---- Track 2nd half of cavity -- loop for all particles. |
|---|
| 1217 | ! do 20 itrack = 1, ktrack |
|---|
| 1218 | ! |
|---|
| 1219 | !*---- Drift to centre. |
|---|
| 1220 | ! px = track(2,itrack) |
|---|
| 1221 | ! py = track(4,itrack) |
|---|
| 1222 | ! pt = track(6,itrack) |
|---|
| 1223 | ! ttt = one/sqrt(one+two*pt*beti+pt**2 - px**2 - py**2) |
|---|
| 1224 | ! track(1,itrack) = track(1,itrack) + dl*ttt*px |
|---|
| 1225 | ! track(3,itrack) = track(3,itrack) + dl*ttt*py |
|---|
| 1226 | ! track(5,itrack) = track(5,itrack) |
|---|
| 1227 | ! + + dl*(beti - (beti+pt)*ttt) + dl*pt*dtbyds |
|---|
| 1228 | ! |
|---|
| 1229 | !*---- Acceleration. |
|---|
| 1230 | ! pt = pt + vrf * sin(phirf - omega * track(5,itrack)) |
|---|
| 1231 | ! track(6,itrack) = pt |
|---|
| 1232 | ! |
|---|
| 1233 | !*---- Drift to end. |
|---|
| 1234 | ! ttt = one/sqrt(one+two*pt*beti+pt**2 - px**2 - py**2) |
|---|
| 1235 | ! track(1,itrack) = track(1,itrack) + dl*ttt*px |
|---|
| 1236 | ! track(3,itrack) = track(3,itrack) + dl*ttt*py |
|---|
| 1237 | ! track(5,itrack) = track(5,itrack) |
|---|
| 1238 | ! + + dl*(beti - (beti+pt)*ttt) + dl*pt*dtbyds |
|---|
| 1239 | ! 20 continue |
|---|
| 1240 | ! endif |
|---|
| 1241 | end subroutine ttrf |
|---|
| 1242 | subroutine ttcrabrf(track,ktrack,turn) |
|---|
| 1243 | |
|---|
| 1244 | implicit none |
|---|
| 1245 | |
|---|
| 1246 | !----------------------------------------------------------------------* |
|---|
| 1247 | ! Purpose: * |
|---|
| 1248 | ! Track a set of trajectories through a thin crab cavity (zero l.) * |
|---|
| 1249 | ! The crab cavity is sandwiched between 2 drift spaces half length. * |
|---|
| 1250 | ! Input/output: * |
|---|
| 1251 | ! TRACK(6,*)(double) Track coordinates: (X, PX, Y, PY, T, PT). * |
|---|
| 1252 | ! KTRACK (integer) number of surviving tracks. * |
|---|
| 1253 | ! Output: * |
|---|
| 1254 | ! EL (double) Length of quadrupole. * |
|---|
| 1255 | !----------------------------------------------------------------------* |
|---|
| 1256 | ! Added: R. Calaga/F. Zimmermann (10/06, 11/07, 09/10) * |
|---|
| 1257 | !----------------------------------------------------------------------* |
|---|
| 1258 | integer itrack,ktrack,turn,t1,t2,t3,t4,p1,p2 |
|---|
| 1259 | double precision omega,phirf,pt,rff,rfl,rfv,eph,track(6,*),clight,twopi,vrf,pc0, & |
|---|
| 1260 | get_variable,node_value,get_value,one,two,half,ten3m,ten6p,px |
|---|
| 1261 | ! double precision px,py,ttt,beti,el1 |
|---|
| 1262 | parameter(one=1d0,two=2d0,half=5d-1,ten3m=1d-3,ten6p=1d6) |
|---|
| 1263 | |
|---|
| 1264 | !---- Initialize |
|---|
| 1265 | clight=get_variable('clight ') |
|---|
| 1266 | twopi=get_variable('twopi ') |
|---|
| 1267 | |
|---|
| 1268 | !---- Fetch data. |
|---|
| 1269 | rfv = node_value('volt ') |
|---|
| 1270 | rff = node_value('freq ') |
|---|
| 1271 | rfl = node_value('lag ') |
|---|
| 1272 | pc0 = get_value('beam ','pc ') |
|---|
| 1273 | |
|---|
| 1274 | !--- specify voltage ramp (relative) |
|---|
| 1275 | t1 = node_value('rv1 ') |
|---|
| 1276 | t2 = node_value('rv2 ') + t1 |
|---|
| 1277 | t3 = node_value('rv3 ') + t2 |
|---|
| 1278 | t4 = node_value('rv4 ') + t3 |
|---|
| 1279 | |
|---|
| 1280 | !--- specify phase change in two steps |
|---|
| 1281 | p1 = node_value('rph1 ') |
|---|
| 1282 | p2 = node_value('rph2 ') + p1 |
|---|
| 1283 | eph = node_value('lagf ') |
|---|
| 1284 | |
|---|
| 1285 | !---- Set up voltage ramp. |
|---|
| 1286 | omega = rff * (ten6p * twopi / clight) |
|---|
| 1287 | |
|---|
| 1288 | if (turn .lt. t1) then |
|---|
| 1289 | vrf = 0.0 |
|---|
| 1290 | else if (turn .ge. t1 .and. turn .lt. t2) then |
|---|
| 1291 | vrf = (turn-t1)*rfv*ten3m/pc0/(t2-t1) |
|---|
| 1292 | else if (turn .ge. t2 .and. turn .lt. t3) then |
|---|
| 1293 | vrf = rfv*ten3m/pc0 |
|---|
| 1294 | else if (turn .ge. t3 .and. turn .lt. t4) then |
|---|
| 1295 | vrf = (t4-turn)*rfv*ten3m/pc0/(t4-t3) |
|---|
| 1296 | else |
|---|
| 1297 | vrf = 0.0 |
|---|
| 1298 | endif |
|---|
| 1299 | |
|---|
| 1300 | !---- Set up phase ramp. |
|---|
| 1301 | if (turn .lt. p1) then |
|---|
| 1302 | phirf = rfl*twopi |
|---|
| 1303 | else if (turn .ge. p1 .and. turn .lt. p2) then |
|---|
| 1304 | phirf = (turn-p1)*eph*twopi/(p2-p1) |
|---|
| 1305 | else |
|---|
| 1306 | phirf= eph*twopi |
|---|
| 1307 | endif |
|---|
| 1308 | |
|---|
| 1309 | ! print*," turn: ",turn, " phase: ", phirf*360/twopi |
|---|
| 1310 | |
|---|
| 1311 | |
|---|
| 1312 | do itrack = 1, ktrack |
|---|
| 1313 | px = track(2,itrack) & |
|---|
| 1314 | + vrf * sin(phirf - omega * track(5,itrack)) |
|---|
| 1315 | pt = track(6,itrack) & |
|---|
| 1316 | - omega* vrf * track(1,itrack) * & |
|---|
| 1317 | cos(phirf - omega * track(5,itrack)) |
|---|
| 1318 | |
|---|
| 1319 | !---- track(2,jtrk) = track(2,jtrk) |
|---|
| 1320 | ! pt = track(6,itrack) |
|---|
| 1321 | ! pt = pt + vrf * sin(phirf - omega * track(5,itrack)) |
|---|
| 1322 | |
|---|
| 1323 | track(2,itrack) = px |
|---|
| 1324 | track(6,itrack) = pt |
|---|
| 1325 | enddo |
|---|
| 1326 | end subroutine ttcrabrf |
|---|
| 1327 | subroutine ttsep(el,track,ktrack) |
|---|
| 1328 | |
|---|
| 1329 | use twtrrfi |
|---|
| 1330 | use trackfi |
|---|
| 1331 | implicit none |
|---|
| 1332 | |
|---|
| 1333 | !----------------------------------------------------------------------* |
|---|
| 1334 | ! Purpose: * |
|---|
| 1335 | ! ttsep tracks particle through an electrostatic separator with an * |
|---|
| 1336 | ! homogeneous, purely perpendicularly acting electric field. The * |
|---|
| 1337 | ! element is given as a thick element and then sliced by makethin. * |
|---|
| 1338 | ! Input: * |
|---|
| 1339 | ! EL (double) Length of the element * |
|---|
| 1340 | ! KTRACK: (integer) Number of surviving tracks * |
|---|
| 1341 | ! Input/Output: * |
|---|
| 1342 | ! TRACK(6,*) (double) Track coordinantes: (X, PX, Y, PY, T, PT) * |
|---|
| 1343 | !----------------------------------------------------------------------* |
|---|
| 1344 | integer ktrack,itrack |
|---|
| 1345 | double precision track(6,*),node_value,get_value |
|---|
| 1346 | double precision el,ex,ey,tilt,charge,mass,pc,beta,deltap,kick0 |
|---|
| 1347 | double precision cos_tilt,sin_tilt,efieldx,efieldy,pt |
|---|
| 1348 | double precision one,ten3m |
|---|
| 1349 | parameter(one=1.0d0,ten3m=1.0d-3) |
|---|
| 1350 | ! |
|---|
| 1351 | ex=node_value('ex_l ') |
|---|
| 1352 | ey=node_value('ey_l ') |
|---|
| 1353 | tilt=node_value('tilt ') |
|---|
| 1354 | cos_tilt=cos(tilt) |
|---|
| 1355 | sin_tilt=sin(tilt) |
|---|
| 1356 | |
|---|
| 1357 | charge=get_value('probe ','charge ') |
|---|
| 1358 | mass =get_value('probe ','mass ') |
|---|
| 1359 | pc =get_value('probe ','pc ') |
|---|
| 1360 | beta =get_value('probe ','beta ') |
|---|
| 1361 | efieldx=ex*cos_tilt + ey*sin_tilt |
|---|
| 1362 | efieldy=-ex*sin_tilt + ey*cos_tilt |
|---|
| 1363 | do itrack = 1, ktrack |
|---|
| 1364 | pt=track(6,itrack) |
|---|
| 1365 | deltap=sqrt(one-one/beta/beta+(pt+one/beta)**2) - one |
|---|
| 1366 | kick0=charge*ten3m/pc/(one+deltap)/beta |
|---|
| 1367 | track(2,itrack) = track(2,itrack) + kick0*efieldx |
|---|
| 1368 | track(4,itrack) = track(4,itrack) + kick0*efieldy |
|---|
| 1369 | end do |
|---|
| 1370 | |
|---|
| 1371 | end subroutine ttsep |
|---|
| 1372 | subroutine ttcorr(el,track,ktrack,turn) |
|---|
| 1373 | |
|---|
| 1374 | use twtrrfi |
|---|
| 1375 | use trackfi |
|---|
| 1376 | implicit none |
|---|
| 1377 | |
|---|
| 1378 | !----------------------------------------------------------------------* |
|---|
| 1379 | ! Purpose: * |
|---|
| 1380 | ! Track particle through an orbit corrector. * |
|---|
| 1381 | ! The corrector is sandwiched between two half-length drift spaces. * |
|---|
| 1382 | ! Input/output: * |
|---|
| 1383 | ! TRACK(6,*)(double) Track coordinates: (X, PX, Y, PY, T, PT). * |
|---|
| 1384 | ! KTRACK (integer) number of surviving tracks. * |
|---|
| 1385 | ! Output: * |
|---|
| 1386 | ! EL (double) Length of quadrupole. * |
|---|
| 1387 | !----------------------------------------------------------------------* |
|---|
| 1388 | ! logical dorad,dodamp,dorand |
|---|
| 1389 | integer itrack,ktrack,n_ferr,node_fd_errors,code,bvk,i,get_option |
|---|
| 1390 | integer sinkick,turn |
|---|
| 1391 | double precision bi2gi2,bil2,curv,dpx,dpy,el,pt,px,py,rfac,rpt, & |
|---|
| 1392 | rpx,rpy,track(6,*),xkick,ykick, & |
|---|
| 1393 | div,f_errors(0:maxferr),field(2),get_variable,get_value, & |
|---|
| 1394 | node_value,zero,one,two,three,half,twopi |
|---|
| 1395 | double precision dpxx,dpyy |
|---|
| 1396 | double precision temp,sinpeak,sintune,sinphase |
|---|
| 1397 | parameter(zero=0d0,one=1d0,two=2d0,three=3d0,half=5d-1) |
|---|
| 1398 | |
|---|
| 1399 | !---- Initialize. |
|---|
| 1400 | twopi=get_variable('twopi ') |
|---|
| 1401 | bvk = node_value('other_bv ') |
|---|
| 1402 | deltas = get_variable('track_deltap ') |
|---|
| 1403 | arad = get_value('probe ','arad ') |
|---|
| 1404 | betas = get_value('probe ','beta ') |
|---|
| 1405 | gammas = get_value('probe ','gamma ') |
|---|
| 1406 | dtbyds = get_value('probe ','dtbyds ') |
|---|
| 1407 | code = node_value('mad8_type ') |
|---|
| 1408 | if(code.eq.39) code=15 |
|---|
| 1409 | if(code.eq.38) code=24 |
|---|
| 1410 | call dzero(f_errors, maxferr+1) |
|---|
| 1411 | n_ferr = node_fd_errors(f_errors) |
|---|
| 1412 | dorad = get_value('probe ','radiate ') .ne. 0 |
|---|
| 1413 | dodamp = get_option('damp ') .ne. 0 |
|---|
| 1414 | dorand = get_option('quantum ') .ne. 0 |
|---|
| 1415 | if (el .eq. zero) then |
|---|
| 1416 | div = one |
|---|
| 1417 | else |
|---|
| 1418 | div = el |
|---|
| 1419 | endif |
|---|
| 1420 | |
|---|
| 1421 | do i = 1, 2 |
|---|
| 1422 | field(i) = zero |
|---|
| 1423 | enddo |
|---|
| 1424 | if (n_ferr .gt. 0) call dcopy(f_errors, field, min(2, n_ferr)) |
|---|
| 1425 | if (code.eq.14) then |
|---|
| 1426 | xkick=bvk*(node_value('kick ')+node_value('chkick ')+ & |
|---|
| 1427 | field(1)/div) |
|---|
| 1428 | ykick=zero |
|---|
| 1429 | else if (code.eq.15) then |
|---|
| 1430 | xkick=bvk*(node_value('hkick ')+node_value('chkick ')+ & |
|---|
| 1431 | field(1)/div) |
|---|
| 1432 | ykick=bvk*(node_value('vkick ')+node_value('cvkick ')+ & |
|---|
| 1433 | field(2)/div) |
|---|
| 1434 | else if(code.eq.16) then |
|---|
| 1435 | xkick=zero |
|---|
| 1436 | ykick=bvk*(node_value('kick ')+node_value('cvkick ')+ & |
|---|
| 1437 | field(2)/div) |
|---|
| 1438 | else |
|---|
| 1439 | xkick=zero |
|---|
| 1440 | ykick=zero |
|---|
| 1441 | endif |
|---|
| 1442 | |
|---|
| 1443 | !---- Sinusoidal kick |
|---|
| 1444 | sinkick = node_value('sinkick ') |
|---|
| 1445 | if(sinkick.eq.1) then |
|---|
| 1446 | sinpeak = node_value('sinpeak ') |
|---|
| 1447 | sintune = node_value('sintune ') |
|---|
| 1448 | sinphase = node_value('sinphase ') |
|---|
| 1449 | temp = sinpeak * sin(sinphase + twopi * sintune * turn) |
|---|
| 1450 | if (code.eq.14) then |
|---|
| 1451 | xkick=xkick+temp |
|---|
| 1452 | else if (code.eq.15) then |
|---|
| 1453 | xkick=xkick+temp |
|---|
| 1454 | ykick=ykick+temp |
|---|
| 1455 | else if(code.eq.16) then |
|---|
| 1456 | ykick=ykick+temp |
|---|
| 1457 | endif |
|---|
| 1458 | endif |
|---|
| 1459 | |
|---|
| 1460 | !---- Sum up total kicks. |
|---|
| 1461 | dpx = xkick / (one + deltas) |
|---|
| 1462 | dpy = ykick / (one + deltas) |
|---|
| 1463 | !---- Thin lens 6D version |
|---|
| 1464 | dpxx = xkick |
|---|
| 1465 | dpyy = ykick |
|---|
| 1466 | |
|---|
| 1467 | bil2 = el / (two * betas) |
|---|
| 1468 | bi2gi2 = one / (betas * gammas) ** 2 |
|---|
| 1469 | |
|---|
| 1470 | !---- Half radiation effects at entrance. |
|---|
| 1471 | if (dorad .and. el .ne. 0) then |
|---|
| 1472 | if (dodamp .and. dorand) then |
|---|
| 1473 | curv = sqrt(dpx**2 + dpy**2) / el |
|---|
| 1474 | else |
|---|
| 1475 | rfac = arad * gammas**3 * (dpx**2 + dpy**2) / (three * el) |
|---|
| 1476 | endif |
|---|
| 1477 | |
|---|
| 1478 | !---- Full damping. |
|---|
| 1479 | if (dodamp) then |
|---|
| 1480 | do itrack = 1, ktrack |
|---|
| 1481 | if (dorand) call trphot(el,curv,rfac,deltas) |
|---|
| 1482 | px = track(2,itrack) |
|---|
| 1483 | py = track(4,itrack) |
|---|
| 1484 | pt = track(6,itrack) |
|---|
| 1485 | track(2,itrack) = px - rfac * (one + pt) * px |
|---|
| 1486 | track(4,itrack) = py - rfac * (one + pt) * py |
|---|
| 1487 | track(6,itrack) = pt - rfac * (one + pt) ** 2 |
|---|
| 1488 | enddo |
|---|
| 1489 | |
|---|
| 1490 | !---- Energy loss as for closed orbit. |
|---|
| 1491 | else |
|---|
| 1492 | rpx = rfac * (one + track(6,1)) * track(2,1) |
|---|
| 1493 | rpy = rfac * (one + track(6,1)) * track(4,1) |
|---|
| 1494 | rpt = rfac * (one + track(6,1)) ** 2 |
|---|
| 1495 | |
|---|
| 1496 | do itrack = 1, ktrack |
|---|
| 1497 | track(2,itrack) = track(2,itrack) - rpx |
|---|
| 1498 | track(4,itrack) = track(4,itrack) - rpy |
|---|
| 1499 | track(6,itrack) = track(6,itrack) - rpt |
|---|
| 1500 | enddo |
|---|
| 1501 | endif |
|---|
| 1502 | endif |
|---|
| 1503 | |
|---|
| 1504 | !---- Thick lens code taken out! AK 21.04.2006 |
|---|
| 1505 | !!---- Half kick at entrance. |
|---|
| 1506 | ! do itrack = 1, ktrack |
|---|
| 1507 | ! px = track(2,itrack) + half * dpx |
|---|
| 1508 | ! py = track(4,itrack) + half * dpy |
|---|
| 1509 | ! pt = track(6,itrack) |
|---|
| 1510 | ! |
|---|
| 1511 | !!---- Drift through corrector. |
|---|
| 1512 | ! d = (one - pt / betas) * el |
|---|
| 1513 | ! track(1,itrack) = track(1,itrack) + px * d |
|---|
| 1514 | ! track(3,itrack) = track(3,itrack) + py * d |
|---|
| 1515 | ! track(5,itrack) = track(5,itrack) + d * bi2gi2 * pt - & |
|---|
| 1516 | ! bil2 * (px**2 + py**2 + bi2gi2*pt**2) + el*pt*dtbyds |
|---|
| 1517 | ! |
|---|
| 1518 | !!---- Half kick at exit. |
|---|
| 1519 | ! track(2,itrack) = px + half * dpx |
|---|
| 1520 | ! track(4,itrack) = py + half * dpy |
|---|
| 1521 | ! track(6,itrack) = pt |
|---|
| 1522 | ! enddo |
|---|
| 1523 | |
|---|
| 1524 | !---- Kick at dipole corrector magnet |
|---|
| 1525 | ! including PT-dependence |
|---|
| 1526 | do itrack = 1, ktrack |
|---|
| 1527 | px = track(2,itrack) |
|---|
| 1528 | py = track(4,itrack) |
|---|
| 1529 | pt = track(6,itrack) |
|---|
| 1530 | ! ttt = sqrt(one+two*pt*bet0i+pt**2 - px**2 - py**2) |
|---|
| 1531 | ! ddd = sqrt(one+two*pt*bet0i+pt**2) |
|---|
| 1532 | ! track(2,itrack) = px + dpxx*ttt/ddd |
|---|
| 1533 | ! track(4,itrack) = py + dpyy*ttt/ddd |
|---|
| 1534 | ! ttt = sqrt(one+two*pt*bet0i+pt**2 - dpxx**2 - dpyy**2) |
|---|
| 1535 | ! ddd = sqrt(one+two*pt*bet0i+pt**2) |
|---|
| 1536 | ! track(2,itrack) = px + dpxx*ttt/ddd |
|---|
| 1537 | ! track(4,itrack) = py + dpyy*ttt/ddd |
|---|
| 1538 | |
|---|
| 1539 | ! ttt = sqrt(one+two*pt*bet0i+pt**2 - px**2 - py**2) |
|---|
| 1540 | ! ddd = sqrt(one+two*pt*bet0i+pt**2) |
|---|
| 1541 | ! |
|---|
| 1542 | ! xp = px/ttt + dpxx/ddd ! Apply kick in standard coordinates! |
|---|
| 1543 | ! yp = py/ttt + dpyy/ddd |
|---|
| 1544 | ! |
|---|
| 1545 | ! ttt = sqrt(one + xp**2 + yp**2) |
|---|
| 1546 | ! |
|---|
| 1547 | ! px = xp*ddd/ttt ! Transform back to canonical coordinates |
|---|
| 1548 | ! py = yp*ddd/ttt |
|---|
| 1549 | ! |
|---|
| 1550 | ! track(2,itrack) = px |
|---|
| 1551 | ! track(4,itrack) = py |
|---|
| 1552 | |
|---|
| 1553 | track(2,itrack) = px + dpxx |
|---|
| 1554 | track(4,itrack) = py + dpyy |
|---|
| 1555 | |
|---|
| 1556 | ! Add time of flight effects (stolen from Ripken-Dipole) |
|---|
| 1557 | ! track(5,itrack) = track(5,itrack) - & |
|---|
| 1558 | ! (dpxx*track(1,itrack) - dpyy*track(3,itrack)) * & |
|---|
| 1559 | ! ((one + bet0*track(6,itrack))/ddd)*bet0i |
|---|
| 1560 | |
|---|
| 1561 | enddo |
|---|
| 1562 | |
|---|
| 1563 | !---- Half radiation effects at exit. |
|---|
| 1564 | ! If not random, use same RFAC as at entrance. |
|---|
| 1565 | if (dorad .and. el .ne. 0) then |
|---|
| 1566 | |
|---|
| 1567 | !---- Full damping. |
|---|
| 1568 | if (dodamp) then |
|---|
| 1569 | do itrack = 1, ktrack |
|---|
| 1570 | if (dorand) call trphot(el,curv,rfac,deltas) |
|---|
| 1571 | px = track(2,itrack) |
|---|
| 1572 | py = track(4,itrack) |
|---|
| 1573 | pt = track(6,itrack) |
|---|
| 1574 | track(2,itrack) = px - rfac * (one + pt) * px |
|---|
| 1575 | track(4,itrack) = py - rfac * (one + pt) * py |
|---|
| 1576 | track(6,itrack) = pt - rfac * (one + pt) ** 2 |
|---|
| 1577 | enddo |
|---|
| 1578 | |
|---|
| 1579 | !---- Energy loss as for closed orbit. |
|---|
| 1580 | else |
|---|
| 1581 | rpx = rfac * (one + track(6,1)) * track(2,1) |
|---|
| 1582 | rpy = rfac * (one + track(6,1)) * track(4,1) |
|---|
| 1583 | rpt = rfac * (one + track(6,1)) ** 2 |
|---|
| 1584 | |
|---|
| 1585 | do itrack = 1, ktrack |
|---|
| 1586 | track(2,itrack) = track(2,itrack) - rpx |
|---|
| 1587 | track(4,itrack) = track(4,itrack) - rpy |
|---|
| 1588 | track(6,itrack) = track(6,itrack) - rpt |
|---|
| 1589 | enddo |
|---|
| 1590 | endif |
|---|
| 1591 | endif |
|---|
| 1592 | |
|---|
| 1593 | end subroutine ttcorr |
|---|
| 1594 | subroutine tt_putone(npart,turn,tot_segm,segment,part_id,z,orbit0,& |
|---|
| 1595 | spos,ielem,el_name) |
|---|
| 1596 | |
|---|
| 1597 | use name_lenfi |
|---|
| 1598 | implicit none |
|---|
| 1599 | |
|---|
| 1600 | !hbu added spos, ielem, el_name |
|---|
| 1601 | !----------------------------------------------------------------------* |
|---|
| 1602 | !--- purpose: enter all particle coordinates in one table * |
|---|
| 1603 | ! input: * |
|---|
| 1604 | ! npart (int) number of particles * |
|---|
| 1605 | ! turn (int) turn number * |
|---|
| 1606 | ! tot_segm (int) total (target) number of entries * |
|---|
| 1607 | ! segment(int) current segment count * |
|---|
| 1608 | ! part_id (int array) particle identifiers * |
|---|
| 1609 | ! z (double (6,*)) particle orbits * |
|---|
| 1610 | ! orbit0 (double array) reference orbit * |
|---|
| 1611 | !----------------------------------------------------------------------* |
|---|
| 1612 | integer i,j,npart,turn,tot_segm,segment,part_id(*),length |
|---|
| 1613 | double precision z(6,*),orbit0(6),tmp,tt,ss |
|---|
| 1614 | !hbu was *36 allow longer info |
|---|
| 1615 | character(120) table,comment |
|---|
| 1616 | !hbu |
|---|
| 1617 | integer ielem |
|---|
| 1618 | !hbu name of element |
|---|
| 1619 | character(name_len) el_name |
|---|
| 1620 | !hbu |
|---|
| 1621 | double precision spos |
|---|
| 1622 | !hbu |
|---|
| 1623 | character(4) vec_names(7) |
|---|
| 1624 | !hbu |
|---|
| 1625 | data vec_names / 'x', 'px', 'y', 'py', 't', 'pt','s' / |
|---|
| 1626 | data table / 'trackone' / |
|---|
| 1627 | |
|---|
| 1628 | !hbu |
|---|
| 1629 | length = len(comment) |
|---|
| 1630 | segment = segment + 1 |
|---|
| 1631 | !hbu |
|---|
| 1632 | write(comment, '(''#segment'',4i8,1X,A)') & |
|---|
| 1633 | segment,tot_segm,npart,ielem,el_name |
|---|
| 1634 | call comment_to_table_curr(table, comment, length) |
|---|
| 1635 | tt = turn |
|---|
| 1636 | do i = 1, npart |
|---|
| 1637 | call double_to_table_curr(table, 'turn ', tt) |
|---|
| 1638 | ss = part_id(i) |
|---|
| 1639 | call double_to_table_curr(table, 'number ', ss) |
|---|
| 1640 | do j = 1, 6 |
|---|
| 1641 | tmp = z(j,i) - orbit0(j) |
|---|
| 1642 | call double_to_table_curr(table, vec_names(j), tmp) |
|---|
| 1643 | enddo |
|---|
| 1644 | !hbu spos |
|---|
| 1645 | call double_to_table_curr(table,vec_names(7),spos) |
|---|
| 1646 | call augment_count(table) |
|---|
| 1647 | enddo |
|---|
| 1648 | end subroutine tt_putone |
|---|
| 1649 | |
|---|
| 1650 | |
|---|
| 1651 | |
|---|
| 1652 | subroutine tt_puttab(npart,turn,nobs,orbit,orbit0,spos) |
|---|
| 1653 | |
|---|
| 1654 | implicit none |
|---|
| 1655 | |
|---|
| 1656 | !hbu added spos |
|---|
| 1657 | !----------------------------------------------------------------------* |
|---|
| 1658 | !--- purpose: enter particle coordinates in table * |
|---|
| 1659 | ! input: * |
|---|
| 1660 | ! npart (int) particle number * |
|---|
| 1661 | ! turn (int) turn number * |
|---|
| 1662 | ! nobs (int) observation point number * |
|---|
| 1663 | ! orbit (double array) particle orbit * |
|---|
| 1664 | ! orbit0 (double array) reference orbit * |
|---|
| 1665 | !----------------------------------------------------------------------* |
|---|
| 1666 | integer npart,turn,j,nobs |
|---|
| 1667 | double precision orbit(6),orbit0(6),tmp,tt,tn |
|---|
| 1668 | double precision energy, get_value |
|---|
| 1669 | character(36) table |
|---|
| 1670 | !hbu |
|---|
| 1671 | double precision spos |
|---|
| 1672 | !hbu |
|---|
| 1673 | character(4) vec_names(7) |
|---|
| 1674 | !hbu |
|---|
| 1675 | data vec_names / 'x', 'px', 'y', 'py', 't', 'pt', 's' / |
|---|
| 1676 | data table / 'track.obs$$$$.p$$$$' / |
|---|
| 1677 | |
|---|
| 1678 | tt = turn |
|---|
| 1679 | tn = npart |
|---|
| 1680 | |
|---|
| 1681 | write(table(10:13), '(i4.4)') nobs |
|---|
| 1682 | write(table(16:19), '(i4.4)') npart |
|---|
| 1683 | |
|---|
| 1684 | energy = get_value('probe ','energy ') |
|---|
| 1685 | |
|---|
| 1686 | call double_to_table_curr(table, 'turn ', tt) ! the number of the cur |
|---|
| 1687 | call double_to_table_curr(table, 'number ', tn) ! the number of the cur |
|---|
| 1688 | call double_to_table_curr(table, 'e ', energy) |
|---|
| 1689 | |
|---|
| 1690 | do j = 1, 6 |
|---|
| 1691 | tmp = orbit(j) - orbit0(j) |
|---|
| 1692 | call double_to_table_curr(table, vec_names(j), tmp) |
|---|
| 1693 | enddo |
|---|
| 1694 | !hbu spos |
|---|
| 1695 | call double_to_table_curr(table,vec_names(7),spos) |
|---|
| 1696 | call augment_count(table) |
|---|
| 1697 | end subroutine tt_puttab |
|---|
| 1698 | |
|---|
| 1699 | subroutine trinicmd(switch,orbit0,eigen,jend,z,turns,coords) |
|---|
| 1700 | |
|---|
| 1701 | use bbfi |
|---|
| 1702 | use trackfi |
|---|
| 1703 | implicit none |
|---|
| 1704 | |
|---|
| 1705 | !----------------------------------------------------------------------* |
|---|
| 1706 | ! Purpose: * |
|---|
| 1707 | ! Define initial conditions for all particles to be tracked * |
|---|
| 1708 | ! input: * |
|---|
| 1709 | ! switch (int) 1: run, 2: dynap fastune, 3: dynap aperture * |
|---|
| 1710 | ! orbit0(6) - closed orbit * |
|---|
| 1711 | ! x, px, y, py, t, deltap, fx, phix, fy, phiy, ft, phit * |
|---|
| 1712 | ! - raw coordinates from start list * |
|---|
| 1713 | ! eigen - Eigenvectors * |
|---|
| 1714 | ! output: * |
|---|
| 1715 | ! jend - number of particles to track * |
|---|
| 1716 | ! z(6,jend) - Transformed cartesian coordinates incl. c.o. * |
|---|
| 1717 | ! coords dp(6,0:turns,npart) (only switch > 1) particle coords. * |
|---|
| 1718 | !----------------------------------------------------------------------* |
|---|
| 1719 | logical zgiv,zngiv |
|---|
| 1720 | integer j,jend,k,kp,kq,next_start,itype(23),switch,turns,jt |
|---|
| 1721 | double precision phi,track(12),zstart(12),twopi,z(6,*),zn(6), & |
|---|
| 1722 | ex,ey,et,orbit0(6),eigen(6,6),x,px,y,py,t,deltae,fx,phix,fy,phiy, & |
|---|
| 1723 | ft,phit,get_value,get_variable,zero,deltax,coords(6,0:turns,*) |
|---|
| 1724 | double precision deltap,one |
|---|
| 1725 | parameter(zero=0d0,one=1d0) |
|---|
| 1726 | character(120) msg(2) |
|---|
| 1727 | |
|---|
| 1728 | deltap = get_variable('track_deltap ') |
|---|
| 1729 | !---- Initialise orbit, emittances and eigenvectors etc. |
|---|
| 1730 | j = 0 |
|---|
| 1731 | twopi = get_variable('twopi ') |
|---|
| 1732 | ex = get_value('probe ','ex ') |
|---|
| 1733 | ey = get_value('probe ','ey ') |
|---|
| 1734 | et = get_value('probe ','et ') |
|---|
| 1735 | bet0 = get_value('beam ','beta ') |
|---|
| 1736 | bet0i = one / bet0 |
|---|
| 1737 | !-----get x add-on for lyaponuv calculation from dynap table |
|---|
| 1738 | deltax = get_value('dynap ', 'lyapunov ') |
|---|
| 1739 | 1 continue |
|---|
| 1740 | jt = next_start(x,px,y,py,t,deltae,fx,phix,fy,phiy,ft,phit) |
|---|
| 1741 | if (switch.gt.1) then |
|---|
| 1742 | j = 2*jt-1 |
|---|
| 1743 | else |
|---|
| 1744 | j = jt |
|---|
| 1745 | endif |
|---|
| 1746 | if (j .ne. -1 .and. j.ne.0) then |
|---|
| 1747 | if (switch .lt. 3 .or. j .eq. 1) then |
|---|
| 1748 | jend = j |
|---|
| 1749 | if (switch.gt.1) jend=jend+1 |
|---|
| 1750 | !---- Get start coordinates |
|---|
| 1751 | track(1) = x |
|---|
| 1752 | track(2) = px |
|---|
| 1753 | track(3) = y |
|---|
| 1754 | track(4) = py |
|---|
| 1755 | track(5) = t |
|---|
| 1756 | ! track(6) = deltae |
|---|
| 1757 | !---- Here we absorb deltap into the PT variable |
|---|
| 1758 | track(6) = deltae + sqrt((one+deltap)**2+bet0i**2-one)-bet0i |
|---|
| 1759 | track(7) = fx |
|---|
| 1760 | track(8) = phix |
|---|
| 1761 | track(9) = fy |
|---|
| 1762 | track(10) = phiy |
|---|
| 1763 | track(11) = ft |
|---|
| 1764 | track(12) = phit |
|---|
| 1765 | do k = 1,12 |
|---|
| 1766 | if(abs(track(k)).ne.zero) then |
|---|
| 1767 | itype(k) = 1 |
|---|
| 1768 | else |
|---|
| 1769 | itype(k) = 0 |
|---|
| 1770 | endif |
|---|
| 1771 | enddo |
|---|
| 1772 | !---- Normalized coordinates. |
|---|
| 1773 | do kq = 1, 5, 2 |
|---|
| 1774 | kp = kq + 1 |
|---|
| 1775 | phi = twopi * track(kq+7) |
|---|
| 1776 | zn(kq) = track(kq+6) * cos(phi) |
|---|
| 1777 | zn(kp) = - track(kq+6) * sin(phi) |
|---|
| 1778 | enddo |
|---|
| 1779 | !---- Transform to unnormalized coordinates and refer to closed orbit. |
|---|
| 1780 | zgiv = .false. |
|---|
| 1781 | zngiv = .false. |
|---|
| 1782 | do k = 1, 6 |
|---|
| 1783 | if (itype(k) .ne. 0) zgiv = .true. |
|---|
| 1784 | if (itype(k+6) .ne. 0) zngiv = .true. |
|---|
| 1785 | zstart(k) = track(k) & |
|---|
| 1786 | + sqrt(ex) * (eigen(k,1) * zn(1) + eigen(k,2) * zn(2)) & |
|---|
| 1787 | + sqrt(ey) * (eigen(k,3) * zn(3) + eigen(k,4) * zn(4)) & |
|---|
| 1788 | + sqrt(et) * (eigen(k,5) * zn(5) + eigen(k,6) * zn(6)) |
|---|
| 1789 | enddo |
|---|
| 1790 | if (switch .gt. 1) then |
|---|
| 1791 | !--- keep initial coordinates for dynap |
|---|
| 1792 | do k = 1, 6 |
|---|
| 1793 | coords(k,0,j) = zstart(k) |
|---|
| 1794 | enddo |
|---|
| 1795 | endif |
|---|
| 1796 | !---- Warn user about possible data conflict. |
|---|
| 1797 | if (zgiv .and. zngiv) then |
|---|
| 1798 | msg(1) = 'Absolute and normalized coordinates given,' |
|---|
| 1799 | msg(2) = 'Superposition used.' |
|---|
| 1800 | call aawarn('START-1: ', msg(1)) |
|---|
| 1801 | call aawarn('START-2: ', msg(2)) |
|---|
| 1802 | endif |
|---|
| 1803 | do k = 1, 6 |
|---|
| 1804 | z(k,j) = orbit0(k) + zstart(k) |
|---|
| 1805 | enddo |
|---|
| 1806 | if (switch.gt.1) then |
|---|
| 1807 | do k = 1, 6 |
|---|
| 1808 | z(k,j+1) = z(k,j) |
|---|
| 1809 | coords(k,0,j+1) = z(k,j+1) |
|---|
| 1810 | enddo |
|---|
| 1811 | z(1,j+1) = z(1,j) + deltax |
|---|
| 1812 | coords(1,0,j+1) = z(1,j+1) |
|---|
| 1813 | endif |
|---|
| 1814 | endif |
|---|
| 1815 | goto 1 |
|---|
| 1816 | endif |
|---|
| 1817 | ! if (switch .eq. 3) then |
|---|
| 1818 | !--- create second particle with x add-on |
|---|
| 1819 | ! deltax = get_value('dynap ', 'lyapunov ') |
|---|
| 1820 | ! jend = 2 |
|---|
| 1821 | ! z(1,jend) = z(1,1) + deltax |
|---|
| 1822 | ! do k = 2, 6 |
|---|
| 1823 | ! z(k,jend) = z(k,1) |
|---|
| 1824 | ! enddo |
|---|
| 1825 | ! endif |
|---|
| 1826 | end subroutine trinicmd |
|---|
| 1827 | subroutine trbegn(rt,eigen) |
|---|
| 1828 | |
|---|
| 1829 | implicit none |
|---|
| 1830 | |
|---|
| 1831 | !----------------------------------------------------------------------* |
|---|
| 1832 | ! Purpose: * |
|---|
| 1833 | ! Initialize tracking mode; TRACK command execution. * |
|---|
| 1834 | !----------------------------------------------------------------------* |
|---|
| 1835 | logical m66sta,onepass |
|---|
| 1836 | integer get_option |
|---|
| 1837 | double precision rt(6,6),reval(6),aival(6),eigen(6,6),zero |
|---|
| 1838 | parameter(zero=0d0) |
|---|
| 1839 | |
|---|
| 1840 | !---- Initialize |
|---|
| 1841 | onepass = get_option('onepass ') .ne. 0 |
|---|
| 1842 | !---- One-pass system: Do not normalize. |
|---|
| 1843 | if (onepass) then |
|---|
| 1844 | call m66one(eigen) |
|---|
| 1845 | else |
|---|
| 1846 | !---- Find eigenvectors. |
|---|
| 1847 | if (m66sta(rt)) then |
|---|
| 1848 | call laseig(rt,reval,aival,eigen) |
|---|
| 1849 | else |
|---|
| 1850 | call ladeig(rt,reval,aival,eigen) |
|---|
| 1851 | endif |
|---|
| 1852 | endif |
|---|
| 1853 | end subroutine trbegn |
|---|
| 1854 | subroutine ttdpdg(track, ktrack) |
|---|
| 1855 | |
|---|
| 1856 | implicit none |
|---|
| 1857 | |
|---|
| 1858 | !----------------------------------------------------------------------* |
|---|
| 1859 | ! Purpose: computes the effect of the dipole edge * |
|---|
| 1860 | ! Input/Output: ktrack , number of surviving trajectories * |
|---|
| 1861 | ! track , trajectory coordinates * |
|---|
| 1862 | !----------------------------------------------------------------------* |
|---|
| 1863 | integer ktrack,itrack |
|---|
| 1864 | double precision fint,e1,h,hgap,corr,rw(6,6),tw(6,6,6),track(6,*),& |
|---|
| 1865 | node_value,zero,one |
|---|
| 1866 | parameter(zero=0d0,one=1d0) |
|---|
| 1867 | |
|---|
| 1868 | call m66one(rw) |
|---|
| 1869 | ! call dzero(tw, 216) |
|---|
| 1870 | |
|---|
| 1871 | e1 = node_value('e1 ') |
|---|
| 1872 | h = node_value('h ') |
|---|
| 1873 | hgap = node_value('hgap ') |
|---|
| 1874 | fint = node_value('fint ') |
|---|
| 1875 | corr = (h + h) * hgap * fint |
|---|
| 1876 | ! print*,"------------------------------------------ " |
|---|
| 1877 | !---- Fringe fields effects computed from the TWISS routine tmfrng |
|---|
| 1878 | ! tmfrng returns the matrix elements rw(used) and tw(unused) |
|---|
| 1879 | ! No radiation effects as it is a pure thin lens with no lrad |
|---|
| 1880 | call tmfrng(.false.,h,zero,e1,zero,zero,corr,rw,tw) |
|---|
| 1881 | do itrack = 1, ktrack |
|---|
| 1882 | track(2,itrack) = track(2,itrack) + rw(2,1)*track(1,itrack) |
|---|
| 1883 | track(4,itrack) = track(4,itrack) + rw(4,3)*track(3,itrack) |
|---|
| 1884 | enddo |
|---|
| 1885 | return |
|---|
| 1886 | end subroutine ttdpdg |
|---|
| 1887 | |
|---|
| 1888 | subroutine trsol(track,ktrack) |
|---|
| 1889 | |
|---|
| 1890 | implicit none |
|---|
| 1891 | |
|---|
| 1892 | !----------------------------------------------------------------------* |
|---|
| 1893 | ! Purpose: * |
|---|
| 1894 | ! Track a set of particles through a thin solenoid. * |
|---|
| 1895 | ! Input/output: * |
|---|
| 1896 | ! TRACK(6,*)(double) Track coordinates: (X, PX, Y, PY, T, PT). * |
|---|
| 1897 | ! KTRACK (integer) number of surviving tracks. * |
|---|
| 1898 | ! Output: * |
|---|
| 1899 | ! EL (double) Length of solenoid. * |
|---|
| 1900 | !----------------------------------------------------------------------* |
|---|
| 1901 | integer itrack,ktrack |
|---|
| 1902 | double precision beta |
|---|
| 1903 | double precision get_value, node_value |
|---|
| 1904 | double precision track(6,*),one,two |
|---|
| 1905 | double precision sk,skl,sks,sksl,cosTh,sinTh,Q,R,Z |
|---|
| 1906 | double precision xf,yf,pxf,pyf,sigf,psigf,bvk |
|---|
| 1907 | double precision onedp,fpsig,fppsig |
|---|
| 1908 | parameter(one=1d0,two=2d0) |
|---|
| 1909 | ! |
|---|
| 1910 | !---- Initialize. |
|---|
| 1911 | ! dtbyds = get_value('probe ','dtbyds ') |
|---|
| 1912 | ! gamma = get_value('probe ','gamma ') |
|---|
| 1913 | beta = get_value('probe ','beta ') |
|---|
| 1914 | ! deltap = get_value('probe ','deltap ') |
|---|
| 1915 | ! |
|---|
| 1916 | !---- Get solenoid parameters |
|---|
| 1917 | ! elrad = node_value('lrad ') |
|---|
| 1918 | sksl = node_value('ksi ') |
|---|
| 1919 | sks = node_value('ks ') |
|---|
| 1920 | |
|---|
| 1921 | !---- BV flag |
|---|
| 1922 | bvk = node_value('other_bv ') |
|---|
| 1923 | sks = sks * bvk |
|---|
| 1924 | sksl = sksl * bvk |
|---|
| 1925 | |
|---|
| 1926 | ! |
|---|
| 1927 | !---- Set up strengths |
|---|
| 1928 | ! sk = sks / two / (one + deltap) |
|---|
| 1929 | sk = sks / two |
|---|
| 1930 | skl = sksl / two |
|---|
| 1931 | ! |
|---|
| 1932 | !---- Loop over particles |
|---|
| 1933 | ! |
|---|
| 1934 | do itrack = 1, ktrack |
|---|
| 1935 | ! |
|---|
| 1936 | ! Ripken formulae p.28 (3.35 and 3.36) |
|---|
| 1937 | xf = track(1,itrack) |
|---|
| 1938 | yf = track(3,itrack) |
|---|
| 1939 | psigf= track(6,itrack)/beta |
|---|
| 1940 | ! |
|---|
| 1941 | ! We do not use a constant deltap!!!!! WE use full 6D formulae! |
|---|
| 1942 | onedp = sqrt(one+2*psigf+(beta**2)*(psigf**2)) |
|---|
| 1943 | fpsig = onedp - one |
|---|
| 1944 | fppsig = (one+(beta**2)*psigf)/onedp |
|---|
| 1945 | ! |
|---|
| 1946 | ! Set up C,S, Q,R,Z |
|---|
| 1947 | cosTh = cos(skl/onedp) |
|---|
| 1948 | sinTh = sin(skl/onedp) |
|---|
| 1949 | Q = -skl*sk/onedp |
|---|
| 1950 | R = fppsig/(onedp**2)*skl*sk |
|---|
| 1951 | Z = fppsig/(onedp**2)*skl |
|---|
| 1952 | ! |
|---|
| 1953 | ! |
|---|
| 1954 | pxf = track(2,itrack) + xf*Q |
|---|
| 1955 | pyf = track(4,itrack) + yf*Q |
|---|
| 1956 | sigf = track(5,itrack)*beta - 0.5D0*(xf**2 + yf**2)*R |
|---|
| 1957 | |
|---|
| 1958 | ! Ripken formulae p.29 (3.37) |
|---|
| 1959 | ! |
|---|
| 1960 | track(1,itrack) = xf*cosTh + yf*sinTh |
|---|
| 1961 | track(2,itrack) = pxf*cosTh + pyf*sinTh |
|---|
| 1962 | track(3,itrack) = -xf*sinTh + yf*cosTh |
|---|
| 1963 | track(4,itrack) = -pxf*sinTh + pyf*cosTh |
|---|
| 1964 | track(5,itrack) = (sigf + (xf*pyf - yf*pxf)*Z)/beta |
|---|
| 1965 | ! track(6,itrack) = psigf*beta |
|---|
| 1966 | |
|---|
| 1967 | enddo |
|---|
| 1968 | end subroutine trsol |
|---|
| 1969 | |
|---|
| 1970 | subroutine tttrans(track,ktrack) |
|---|
| 1971 | |
|---|
| 1972 | implicit none |
|---|
| 1973 | |
|---|
| 1974 | !----------------------------------------------------------------------* |
|---|
| 1975 | ! Purpose: * |
|---|
| 1976 | ! Translation of a set of particles. * |
|---|
| 1977 | ! Input/output: * |
|---|
| 1978 | ! TRACK(6,*)(double) Track coordinates: (X, PX, Y, PY, T, PT). * |
|---|
| 1979 | !----------------------------------------------------------------------* |
|---|
| 1980 | integer itrack,ktrack |
|---|
| 1981 | double precision node_value |
|---|
| 1982 | double precision track(6,*) |
|---|
| 1983 | double precision t_x,t_px,t_y,t_py,t_sig,t_psig |
|---|
| 1984 | ! |
|---|
| 1985 | !---- Initialize. |
|---|
| 1986 | ! |
|---|
| 1987 | !---- Get translation parameters |
|---|
| 1988 | t_x = node_value('x ') |
|---|
| 1989 | t_px = node_value('px ') |
|---|
| 1990 | t_y = node_value('y ') |
|---|
| 1991 | t_py = node_value('py ') |
|---|
| 1992 | t_sig = node_value('t ') |
|---|
| 1993 | t_psig = node_value('pt ') |
|---|
| 1994 | ! |
|---|
| 1995 | !---- Loop over particles |
|---|
| 1996 | ! |
|---|
| 1997 | do itrack = 1, ktrack |
|---|
| 1998 | ! |
|---|
| 1999 | ! Add vector to particle coordinates |
|---|
| 2000 | track(1,itrack) = track(1,itrack) + t_x |
|---|
| 2001 | track(2,itrack) = track(2,itrack) + t_px |
|---|
| 2002 | track(3,itrack) = track(3,itrack) + t_y |
|---|
| 2003 | track(4,itrack) = track(4,itrack) + t_py |
|---|
| 2004 | track(5,itrack) = track(5,itrack) + t_sig |
|---|
| 2005 | track(6,itrack) = track(6,itrack) + t_psig |
|---|
| 2006 | ! |
|---|
| 2007 | enddo |
|---|
| 2008 | end subroutine tttrans |
|---|
| 2009 | |
|---|
| 2010 | subroutine tttrak(ek,re,track,ktrack) |
|---|
| 2011 | |
|---|
| 2012 | implicit none |
|---|
| 2013 | |
|---|
| 2014 | !----------------------------------------------------------------------* |
|---|
| 2015 | ! Purpose: * |
|---|
| 2016 | ! Track a set of particle with a given TRANSPORT map. * |
|---|
| 2017 | ! Input: * |
|---|
| 2018 | ! EK(6) (real) Kick. * |
|---|
| 2019 | ! RE(6,6) (real) First-order terms. * |
|---|
| 2020 | ! Input/output: * |
|---|
| 2021 | ! TRACK(6,*)(real) Track coordinates: (X, PX, Y, PY, T, PT). * |
|---|
| 2022 | ! KTRACK (integer) number of surviving tracks. * |
|---|
| 2023 | !----------------------------------------------------------------------* |
|---|
| 2024 | integer i,itrack,ktrack |
|---|
| 2025 | double precision ek(6),re(36),temp(6),track(6,*) |
|---|
| 2026 | ! double precision se(36) |
|---|
| 2027 | |
|---|
| 2028 | do itrack = 1, ktrack |
|---|
| 2029 | ! do i = 1, 36 |
|---|
| 2030 | ! se(i) = re(i) & |
|---|
| 2031 | ! + te(i,1)*track(1,itrack) + te(i,2) * track(2,itrack) & |
|---|
| 2032 | ! + te(i,3)*track(3,itrack) + te(i,4) * track(4,itrack) & |
|---|
| 2033 | ! + te(i,5)*track(5,itrack) + te(i,6) * track(6,itrack) |
|---|
| 2034 | ! enddo |
|---|
| 2035 | do i = 1, 6 |
|---|
| 2036 | temp(i) = ek(i) & |
|---|
| 2037 | + re(i) * track(1,itrack) + re(i+ 6) * track(2,itrack) & |
|---|
| 2038 | + re(i+12) * track(3,itrack) + re(i+18) * track(4,itrack) & |
|---|
| 2039 | + re(i+24) * track(5,itrack) + re(i+30) * track(6,itrack) |
|---|
| 2040 | enddo |
|---|
| 2041 | do i = 1, 6 |
|---|
| 2042 | track(i,itrack) = temp(i) |
|---|
| 2043 | enddo |
|---|
| 2044 | enddo |
|---|
| 2045 | end subroutine tttrak |
|---|
| 2046 | |
|---|
| 2047 | subroutine trupdate(turn) |
|---|
| 2048 | |
|---|
| 2049 | implicit none |
|---|
| 2050 | |
|---|
| 2051 | !----------------------------------------------------------------------* |
|---|
| 2052 | ! Purpose: * |
|---|
| 2053 | ! Update global turnnumber-VARIABLE, * |
|---|
| 2054 | ! which can be used to make other variables turn-dependent. * |
|---|
| 2055 | ! Input/output: * |
|---|
| 2056 | ! turn (integer) Current turn number. * |
|---|
| 2057 | !----------------------------------------------------------------------* |
|---|
| 2058 | integer turn |
|---|
| 2059 | character(50) cmd |
|---|
| 2060 | |
|---|
| 2061 | !---- call pro_input('TR$TURN := turn;') |
|---|
| 2062 | write(cmd, '(''tr$turni := '',i8,'' ; '')') turn |
|---|
| 2063 | call pro_input(cmd) |
|---|
| 2064 | write(cmd, '(''exec, tr$macro($tr$turni) ; '')') |
|---|
| 2065 | call pro_input(cmd) |
|---|
| 2066 | end subroutine trupdate |
|---|
| 2067 | |
|---|
| 2068 | subroutine trclor(orbit0) |
|---|
| 2069 | |
|---|
| 2070 | use twiss0fi |
|---|
| 2071 | use name_lenfi |
|---|
| 2072 | use trackfi |
|---|
| 2073 | implicit none |
|---|
| 2074 | |
|---|
| 2075 | !----------------------------------------------------------------------* |
|---|
| 2076 | ! Purpose: * |
|---|
| 2077 | ! - Find 6D closed orbit * |
|---|
| 2078 | ! - Use Newton approximation * |
|---|
| 2079 | ! * |
|---|
| 2080 | ! Input/output: * |
|---|
| 2081 | ! ORBIT0(6) (double) Closed Orbit Coordinates @ sequence start * |
|---|
| 2082 | ! (X, PX, Y, PY, T, PT) * |
|---|
| 2083 | !----------------------------------------------------------------------* |
|---|
| 2084 | double precision zero, one |
|---|
| 2085 | parameter(zero=0d0,one=1d0) |
|---|
| 2086 | double precision orbit0(6),z(6,7),zz(6),z0(6,7),z00(6,7),a(6,7),deltap,ddd(6) |
|---|
| 2087 | integer itra, itmax, j, bbd_pos, j_tot |
|---|
| 2088 | integer code |
|---|
| 2089 | double precision el,dxt(200),dyt(200) |
|---|
| 2090 | integer restart_sequ |
|---|
| 2091 | integer advance_node, get_option, node_al_errors |
|---|
| 2092 | double precision node_value, get_value, get_variable |
|---|
| 2093 | integer n_align |
|---|
| 2094 | double precision al_errors(align_max) |
|---|
| 2095 | double precision re(6,6) |
|---|
| 2096 | |
|---|
| 2097 | logical aperflag, onepass |
|---|
| 2098 | integer pmax, turn, turnmax, part_id(1),last_turn(1) |
|---|
| 2099 | double precision sum, orbit(6) |
|---|
| 2100 | double precision last_pos(6),last_orbit(6,1),maxaper(6) |
|---|
| 2101 | |
|---|
| 2102 | integer nint, ndble, nchar, int_arr(1), char_l |
|---|
| 2103 | character(12) char_a |
|---|
| 2104 | |
|---|
| 2105 | integer i,k, irank |
|---|
| 2106 | |
|---|
| 2107 | double precision cotol, err |
|---|
| 2108 | |
|---|
| 2109 | logical is_drift, is_thin, is_quad, is_matrix |
|---|
| 2110 | |
|---|
| 2111 | print *," " |
|---|
| 2112 | ! print *," AK special version 2007/12/13" |
|---|
| 2113 | ! print *," =============================" |
|---|
| 2114 | |
|---|
| 2115 | ! print *," Modified by Yipeng SUN 19/11/2008 " |
|---|
| 2116 | ! print *," =======================" |
|---|
| 2117 | print *," Full 6D closed orbit search." |
|---|
| 2118 | print*," Initial value of 6-D closed orbit from Twiss: " |
|---|
| 2119 | print*," orbit0 ",orbit0 |
|---|
| 2120 | |
|---|
| 2121 | !---- Initialize needed variables |
|---|
| 2122 | turn = 1 |
|---|
| 2123 | turnmax = 1 |
|---|
| 2124 | itmax = 10 |
|---|
| 2125 | pmax = 7 |
|---|
| 2126 | sum = zero |
|---|
| 2127 | aperflag = .false. |
|---|
| 2128 | onepass = .true. |
|---|
| 2129 | |
|---|
| 2130 | do k=1,7 |
|---|
| 2131 | call dcopy(orbit0,z(1,k),6) |
|---|
| 2132 | enddo |
|---|
| 2133 | |
|---|
| 2134 | ! ddd(1) = orbit0(1)/100000 |
|---|
| 2135 | ! ddd(2) = orbit0(2)/100000 |
|---|
| 2136 | ! ddd(3) = orbit0(3)/100000 |
|---|
| 2137 | ! ddd(4) = orbit0(4)/100000 |
|---|
| 2138 | ! ddd(5) = orbit0(5)/100000 |
|---|
| 2139 | ! ddd(6) = orbit0(6)/100000 |
|---|
| 2140 | |
|---|
| 2141 | ! for on momentum pt =0 |
|---|
| 2142 | ! ddd(1) = 1e-15 |
|---|
| 2143 | ! ddd(2) = 1e-15 |
|---|
| 2144 | ! ddd(3) = 1e-15 |
|---|
| 2145 | ! ddd(4) = 1e-15 |
|---|
| 2146 | ! ddd(5) = 1e-15 |
|---|
| 2147 | ! ddd(6) = 1e-15 |
|---|
| 2148 | !-------------------------------------- |
|---|
| 2149 | ddd(1) = 1e-15 |
|---|
| 2150 | ddd(2) = 1e-15 |
|---|
| 2151 | ddd(3) = 1e-15 |
|---|
| 2152 | ddd(4) = 1e-15 |
|---|
| 2153 | ddd(5) = 1e-15 |
|---|
| 2154 | ddd(6) = 1e-15 |
|---|
| 2155 | |
|---|
| 2156 | ! do k=1,6 |
|---|
| 2157 | ! z(k,k+1) = z(k,k+1) + ddd(k) |
|---|
| 2158 | ! enddo |
|---|
| 2159 | |
|---|
| 2160 | do k=1,7 |
|---|
| 2161 | call dcopy(z(1,k),z0(1,k),6) |
|---|
| 2162 | call dcopy(z(1,k),z00(1,k),6) |
|---|
| 2163 | enddo |
|---|
| 2164 | |
|---|
| 2165 | !--- jmax may be reduced by particle loss - keep number in j_tot |
|---|
| 2166 | j_tot = pmax |
|---|
| 2167 | !--- get vector of six coordinate maxapers (both RUN and DYNAP) |
|---|
| 2168 | call comm_para('maxaper ',nint,ndble,nchar,int_arr,maxaper, & |
|---|
| 2169 | char_a, char_l) |
|---|
| 2170 | !--- set particle id |
|---|
| 2171 | |
|---|
| 2172 | |
|---|
| 2173 | |
|---|
| 2174 | |
|---|
| 2175 | cotol = get_variable('twiss_tol ') |
|---|
| 2176 | |
|---|
| 2177 | !---- Initialize kinematics and orbit |
|---|
| 2178 | bet0 = get_value('beam ','beta ') |
|---|
| 2179 | betas = get_value('probe ','beta ') |
|---|
| 2180 | gammas= get_value('probe ','gamma ') |
|---|
| 2181 | bet0i = one / bet0 |
|---|
| 2182 | beti = one / betas |
|---|
| 2183 | dtbyds = get_value('probe ','dtbyds ') |
|---|
| 2184 | deltas = get_variable('track_deltap ') |
|---|
| 2185 | deltap = get_value('probe ','deltap ') |
|---|
| 2186 | arad = get_value('probe ','arad ') |
|---|
| 2187 | dorad = get_value('probe ','radiate ') .ne. 0 |
|---|
| 2188 | dodamp = get_option('damp ') .ne. 0 |
|---|
| 2189 | dorand = get_option('quantum ') .ne. 0 |
|---|
| 2190 | call dcopy(orbit0,orbit,6) |
|---|
| 2191 | |
|---|
| 2192 | |
|---|
| 2193 | !---- Iteration for closed orbit. |
|---|
| 2194 | do itra = 1, itmax |
|---|
| 2195 | ! print*," before tracking " |
|---|
| 2196 | ! print*," itra: ",itra, " z: ", z |
|---|
| 2197 | j = restart_sequ() |
|---|
| 2198 | |
|---|
| 2199 | !---- loop over nodes |
|---|
| 2200 | 10 continue |
|---|
| 2201 | |
|---|
| 2202 | |
|---|
| 2203 | bbd_pos = j |
|---|
| 2204 | code = node_value('mad8_type ') |
|---|
| 2205 | if(code.eq.39) code=15 |
|---|
| 2206 | if(code.eq.38) code=24 |
|---|
| 2207 | el = node_value('l ') |
|---|
| 2208 | if (itra .eq. 1) then |
|---|
| 2209 | if (.not.(is_drift() .or. is_thin() .or. is_quad() .or. is_matrix())) then |
|---|
| 2210 | print*," " |
|---|
| 2211 | print*,"code: ",code," el: ",el," THICK ELEMENT FOUND" |
|---|
| 2212 | print*," " |
|---|
| 2213 | print*,"Track dies nicely" |
|---|
| 2214 | print*,"Thick lenses will get nowhere" |
|---|
| 2215 | print*,"MAKETHIN will save you" |
|---|
| 2216 | print*," " |
|---|
| 2217 | print*," " |
|---|
| 2218 | stop |
|---|
| 2219 | endif |
|---|
| 2220 | endif |
|---|
| 2221 | |
|---|
| 2222 | |
|---|
| 2223 | |
|---|
| 2224 | !-------- Misalignment at beginning of element (from twissfs.f) |
|---|
| 2225 | if (code .ne. 1) then |
|---|
| 2226 | call dzero(al_errors, align_max) |
|---|
| 2227 | n_align = node_al_errors(al_errors) |
|---|
| 2228 | if (n_align .ne. 0) then |
|---|
| 2229 | do i = 1, pmax |
|---|
| 2230 | call dcopy(z(1,i),zz,6) |
|---|
| 2231 | call tmali1(zz,al_errors, betas, gammas,z(1,i), re) |
|---|
| 2232 | enddo |
|---|
| 2233 | endif |
|---|
| 2234 | endif |
|---|
| 2235 | |
|---|
| 2236 | !-------- Track through element |
|---|
| 2237 | call ttmap(code,el,z,pmax,dxt,dyt,sum,turn,part_id, & |
|---|
| 2238 | last_turn,last_pos, last_orbit,aperflag,maxaper,al_errors,onepass) |
|---|
| 2239 | |
|---|
| 2240 | !-------- Misalignment at end of element (from twissfs.f) |
|---|
| 2241 | if (code .ne. 1) then |
|---|
| 2242 | if (n_align .ne. 0) then |
|---|
| 2243 | do i = 1, pmax |
|---|
| 2244 | call dcopy(z(1,i),zz,6) |
|---|
| 2245 | call tmali2(el,zz, al_errors, betas, gammas,z(1,i), re) |
|---|
| 2246 | enddo |
|---|
| 2247 | endif |
|---|
| 2248 | endif |
|---|
| 2249 | |
|---|
| 2250 | |
|---|
| 2251 | if (advance_node().ne.0) then |
|---|
| 2252 | j=j+1 |
|---|
| 2253 | goto 10 |
|---|
| 2254 | endif |
|---|
| 2255 | !---- end of loop over nodes |
|---|
| 2256 | |
|---|
| 2257 | |
|---|
| 2258 | !---- construct one-turn map |
|---|
| 2259 | do k=1,6 |
|---|
| 2260 | do i=1,6 |
|---|
| 2261 | a(i,k) = (z(i,k+1) - z(i,1))/ddd(i) |
|---|
| 2262 | enddo |
|---|
| 2263 | enddo |
|---|
| 2264 | !---- Solve for dynamic case. |
|---|
| 2265 | err = zero |
|---|
| 2266 | do i= 1,6 |
|---|
| 2267 | a(i,i) = a(i,i) - one |
|---|
| 2268 | a(i,7) = z(i,1) - z0(i,1) |
|---|
| 2269 | err = max(abs(a(i,7)), err) |
|---|
| 2270 | enddo |
|---|
| 2271 | |
|---|
| 2272 | call solver(a,6,1,irank) |
|---|
| 2273 | if (irank.lt.6) go to 100 |
|---|
| 2274 | do i = 1, 6 |
|---|
| 2275 | z0(i,1) = z0(i,1) - a(i,7) |
|---|
| 2276 | z00(i,1) = z00(i,1) - a(i,7) |
|---|
| 2277 | enddo |
|---|
| 2278 | !---- Solve for dynamic case. |
|---|
| 2279 | ! do i = 1, 6 |
|---|
| 2280 | ! print*," a(i,7) ",a(i,7) |
|---|
| 2281 | ! enddo |
|---|
| 2282 | |
|---|
| 2283 | ! if (err.lt.cotol) then |
|---|
| 2284 | |
|---|
| 2285 | ! print *, ' ' |
|---|
| 2286 | ! print '(''iteration: '',i3,'' error: '',1p,e14.6,'' deltap: '' & |
|---|
| 2287 | ! , & |
|---|
| 2288 | ! 1p,e14.6)',itra,err,deltap |
|---|
| 2289 | ! print '(''orbit: '', 1p,6e14.6)', z0 |
|---|
| 2290 | |
|---|
| 2291 | ! goto 110 |
|---|
| 2292 | |
|---|
| 2293 | ! endif |
|---|
| 2294 | |
|---|
| 2295 | |
|---|
| 2296 | |
|---|
| 2297 | do k=2,7 |
|---|
| 2298 | call dcopy(z00(1,1),z0(1,k),6) |
|---|
| 2299 | enddo |
|---|
| 2300 | |
|---|
| 2301 | do k=1,6 |
|---|
| 2302 | z0(k,k+1) = z0(k,k+1) + ddd(k) |
|---|
| 2303 | enddo |
|---|
| 2304 | |
|---|
| 2305 | do k=1,7 |
|---|
| 2306 | call dcopy(z0(1,k),z(1,k),6) |
|---|
| 2307 | enddo |
|---|
| 2308 | |
|---|
| 2309 | |
|---|
| 2310 | !---- end of Iteration |
|---|
| 2311 | enddo |
|---|
| 2312 | |
|---|
| 2313 | call dcopy(z0(1,1),orbit0,6) |
|---|
| 2314 | ! print *,"madX::trrun.F::trclor()" |
|---|
| 2315 | ! print *," Iteration maximum reached. NO convergence." |
|---|
| 2316 | goto 110 |
|---|
| 2317 | |
|---|
| 2318 | 100 continue |
|---|
| 2319 | print *," Singular matrix occurred during closed orbit search." |
|---|
| 2320 | |
|---|
| 2321 | 110 continue |
|---|
| 2322 | print *, ' ' |
|---|
| 2323 | print *, '6D closed orbit found by subroutine trclor ' |
|---|
| 2324 | print '(''iteration: '',i3,'' error: '',1p,e14.6,'' deltap: '',1p,e14.6)',itra,err,deltap |
|---|
| 2325 | print '(''orbit: '', 1p,6e14.6)', orbit0 |
|---|
| 2326 | |
|---|
| 2327 | end subroutine trclor |
|---|
| 2328 | |
|---|
| 2329 | subroutine ttnllens(track,ktrack) |
|---|
| 2330 | |
|---|
| 2331 | implicit none |
|---|
| 2332 | |
|---|
| 2333 | !----------------------------------------------------------------------* |
|---|
| 2334 | ! Purpose: * |
|---|
| 2335 | ! Track a set of trajectories through a thin nonlinear lens * |
|---|
| 2336 | ! Input/output: * |
|---|
| 2337 | ! TRACK(6,*)(double) Track coordinates: (X, PX, Y, PY, T, PT). * |
|---|
| 2338 | ! KTRACK (integer) number of surviving tracks. * |
|---|
| 2339 | !----------------------------------------------------------------------* |
|---|
| 2340 | ! Added by Alexander Valishev 22 Oct 2010 * |
|---|
| 2341 | !----------------------------------------------------------------------* |
|---|
| 2342 | integer itrack,ktrack |
|---|
| 2343 | double precision track(6,*),node_value,get_variable, & |
|---|
| 2344 | pi,one,two,knll,cnll, & |
|---|
| 2345 | dd, u, v, dUu, dUv, dux, duy, dvx, dvy, x, y |
|---|
| 2346 | parameter(one=1d0,two=2d0) |
|---|
| 2347 | |
|---|
| 2348 | cnll=node_value('cnll ') |
|---|
| 2349 | knll=node_value('knll ')/cnll |
|---|
| 2350 | pi=get_variable('pi ') |
|---|
| 2351 | |
|---|
| 2352 | do itrack = 1, ktrack |
|---|
| 2353 | |
|---|
| 2354 | x = track(1,itrack)/cnll |
|---|
| 2355 | y = track(3,itrack)/cnll |
|---|
| 2356 | |
|---|
| 2357 | u=0.5*sqrt((x-1)**2+y**2)+0.5*sqrt((x+1)**2+y**2) |
|---|
| 2358 | v=0.5*sqrt((x+1)**2+y**2)-0.5*sqrt((x-1)**2+y**2) |
|---|
| 2359 | if (u.eq.1d0) then |
|---|
| 2360 | dd=0 |
|---|
| 2361 | else |
|---|
| 2362 | dd=u**2*log(u+sqrt(u*u-1))/sqrt(u**2-1) |
|---|
| 2363 | endif |
|---|
| 2364 | dUu=(u+log(u+sqrt(u*u-1))*sqrt(u**2-1)+dd)/(u**2-v**2) & |
|---|
| 2365 | -2*u*(u*log(u+sqrt(u*u-1))*sqrt(u**2-1) & |
|---|
| 2366 | +v*(acos(v)-0.5*pi)*sqrt(1-v**2)) /(u**2-v**2)**2 |
|---|
| 2367 | dUv=2*v*(u*log(u+sqrt(u*u-1))*sqrt(u**2-1) & |
|---|
| 2368 | +v*(acos(v)-0.5*pi)*sqrt(1-v**2)) /(u**2-v**2)**2 & |
|---|
| 2369 | -(v-(acos(v)-0.5*pi)*sqrt(1-v**2)+v**2*(acos(v)-0.5*pi)/sqrt(1-v**2))& |
|---|
| 2370 | /(u**2-v**2) |
|---|
| 2371 | dux=0.5*(x-1)/sqrt((x-1)**2+y**2) +0.5*(x+1)/sqrt((x+1)**2+y**2) |
|---|
| 2372 | duy=0.5*y/sqrt((x-1)**2+y**2) +0.5*y/sqrt((x+1)**2+y**2) |
|---|
| 2373 | dvx=0.5*(x+1)/sqrt((x+1)**2+y**2) -0.5*(x-1)/sqrt((x-1)**2+y**2) |
|---|
| 2374 | dvy=0.5*y/sqrt((x+1)**2+y**2) -0.5*y/sqrt((x-1)**2+y**2) |
|---|
| 2375 | |
|---|
| 2376 | track(2,itrack)=track(2,itrack)+knll*(dUu*dux+dUv*dvx) |
|---|
| 2377 | track(4,itrack)=track(4,itrack)+knll*(dUu*duy+dUv*dvy) |
|---|
| 2378 | |
|---|
| 2379 | enddo |
|---|
| 2380 | end subroutine ttnllens |
|---|
| 2381 | |
|---|
| 2382 | subroutine ttrfmult(track, ktrack, turn) |
|---|
| 2383 | |
|---|
| 2384 | use twtrrfi |
|---|
| 2385 | use trackfi |
|---|
| 2386 | implicit none |
|---|
| 2387 | |
|---|
| 2388 | !--------------------* |
|---|
| 2389 | ! Andrea Latina 2012 * |
|---|
| 2390 | !--------------------* |
|---|
| 2391 | !----------------------------------------------------------------------* |
|---|
| 2392 | ! Purpose: * |
|---|
| 2393 | ! Track particle through a general thin rf-multipole. * |
|---|
| 2394 | ! Input/output: * |
|---|
| 2395 | ! TRACK(6,*)(double) Track coordinates: (X, PX, Y, PY, T, PT). * |
|---|
| 2396 | ! KTRACK (integer) Number of surviving tracks. * |
|---|
| 2397 | !----------------------------------------------------------------------* |
|---|
| 2398 | |
|---|
| 2399 | double precision beta, angle, cangle, sangle, dtmp |
|---|
| 2400 | double precision node_value, get_value |
|---|
| 2401 | double precision bvk, deltap, elrad |
|---|
| 2402 | double precision f_errors(0:maxferr) |
|---|
| 2403 | double precision vals(2,0:maxmul) |
|---|
| 2404 | integer n_ferr, jtrk, iord, nord |
|---|
| 2405 | |
|---|
| 2406 | !--- AL: RF-multipole |
|---|
| 2407 | integer dummyi, nn, ns, node_fd_errors, turn, ktrack |
|---|
| 2408 | double precision normal(0:maxmul), skew(0:maxmul) |
|---|
| 2409 | double precision track(6,*), field(2,0:maxmul) |
|---|
| 2410 | double precision field_cos(2,0:maxmul) |
|---|
| 2411 | double precision field_sin(2,0:maxmul) |
|---|
| 2412 | double precision zero, one, two, three |
|---|
| 2413 | double precision pc, krf, rfac |
|---|
| 2414 | double precision pi, clight, ten3m |
|---|
| 2415 | double precision x, y, z, dpx, dpy, dpt |
|---|
| 2416 | double precision freq, volt, lag, harmon |
|---|
| 2417 | double precision pnl(0:maxmul), psl(0:maxmul) |
|---|
| 2418 | complex*16 ii, Cm2, Sm2, Cm1, Sm1, Cp0, Sp0, Cp1, Sp1 |
|---|
| 2419 | |
|---|
| 2420 | parameter ( pi = 3.14159265358979d0 ) |
|---|
| 2421 | parameter ( clight = 299792458d0 ) |
|---|
| 2422 | parameter ( zero=0d0, one=1d0, two=2d0, three=3d0, ten3m=1d-3) |
|---|
| 2423 | parameter ( ii=(0d0,1d0) ) |
|---|
| 2424 | |
|---|
| 2425 | !---- Zero the arrays |
|---|
| 2426 | call dzero(normal,maxmul+1) |
|---|
| 2427 | call dzero(skew,maxmul+1) |
|---|
| 2428 | call dzero(pnl,maxmul+1) |
|---|
| 2429 | call dzero(psl,maxmul+1) |
|---|
| 2430 | call dzero(f_errors,maxferr+1) |
|---|
| 2431 | call dzero(field,2*(maxmul+1)) |
|---|
| 2432 | |
|---|
| 2433 | !---- Read-in the parameters |
|---|
| 2434 | freq = node_value('freq '); |
|---|
| 2435 | volt = node_value('volt '); |
|---|
| 2436 | lag = node_value('lag '); |
|---|
| 2437 | harmon = node_value('harmon '); |
|---|
| 2438 | bvk = node_value('other_bv ') |
|---|
| 2439 | elrad = node_value('lrad ') |
|---|
| 2440 | deltap = get_value('probe ', 'deltap ') |
|---|
| 2441 | dorad = get_value('probe ','radiate ') .ne. zero |
|---|
| 2442 | arad = get_value('probe ','arad ') |
|---|
| 2443 | gammas = get_value('probe ','gamma ') |
|---|
| 2444 | beta = get_value('probe ','beta ') |
|---|
| 2445 | pc = get_value('probe ','pc ') |
|---|
| 2446 | n_ferr = node_fd_errors(f_errors); |
|---|
| 2447 | call get_node_vector('knl ', nn, normal) |
|---|
| 2448 | call get_node_vector('ksl ', ns, skew) |
|---|
| 2449 | call get_node_vector('pnl ', dummyi, pnl); |
|---|
| 2450 | call get_node_vector('psl ', dummyi, psl); |
|---|
| 2451 | |
|---|
| 2452 | rfac = zero |
|---|
| 2453 | |
|---|
| 2454 | !---- Set-up some parameters |
|---|
| 2455 | krf = 2*pi*freq*1d6/clight; |
|---|
| 2456 | |
|---|
| 2457 | if (n_ferr.gt.0) then |
|---|
| 2458 | call dcopy(f_errors,field,n_ferr) |
|---|
| 2459 | endif |
|---|
| 2460 | nord = max(nn, ns, n_ferr/2-1); |
|---|
| 2461 | |
|---|
| 2462 | !---- Prepare to calculate the kick and the matrix elements |
|---|
| 2463 | do jtrk = 1,ktrack |
|---|
| 2464 | x = track(1,jtrk); |
|---|
| 2465 | y = track(3,jtrk); |
|---|
| 2466 | z = track(5,jtrk); |
|---|
| 2467 | !---- Vector with strengths + field errors |
|---|
| 2468 | do iord = 0, nord; |
|---|
| 2469 | field_cos(1,iord) = bvk * (normal(iord) * cos(pnl(iord) * 2 * pi - krf * z) + field(1,iord)); |
|---|
| 2470 | field_sin(1,iord) = bvk * (normal(iord) * sin(pnl(iord) * 2 * pi - krf * z)); |
|---|
| 2471 | field_cos(2,iord) = bvk * (skew(iord) * cos(psl(iord) * 2 * pi - krf * z) + field(2,iord)); |
|---|
| 2472 | field_sin(2,iord) = bvk * (skew(iord) * sin(psl(iord) * 2 * pi - krf * z)); |
|---|
| 2473 | enddo |
|---|
| 2474 | Cm2 = 0d0; |
|---|
| 2475 | Sm2 = 0d0; |
|---|
| 2476 | Cm1 = 0d0; |
|---|
| 2477 | Sm1 = 0d0; |
|---|
| 2478 | Cp0 = 0d0; |
|---|
| 2479 | Sp0 = 0d0; |
|---|
| 2480 | Cp1 = 0d0; |
|---|
| 2481 | Sp1 = 0d0; |
|---|
| 2482 | do iord = nord, 0, -1 |
|---|
| 2483 | if (iord.ge.2) then |
|---|
| 2484 | Cm2 = Cm2 * (x+ii*y) / (iord-1) + field_cos(1,iord)+ii*field_cos(2,iord); |
|---|
| 2485 | Sm2 = Sm2 * (x+ii*y) / (iord-1) + field_sin(1,iord)+ii*field_sin(2,iord); |
|---|
| 2486 | endif |
|---|
| 2487 | if (iord.ge.1) then |
|---|
| 2488 | Cm1 = Cm1 * (x+ii*y) / (iord) + field_cos(1,iord)+ii*field_cos(2,iord); |
|---|
| 2489 | Sm1 = Sm1 * (x+ii*y) / (iord) + field_sin(1,iord)+ii*field_sin(2,iord); |
|---|
| 2490 | endif |
|---|
| 2491 | Cp0 = Cp0 * (x+ii*y) / (iord+1) + field_cos(1,iord)+ii*field_cos(2,iord); |
|---|
| 2492 | Sp0 = Sp0 * (x+ii*y) / (iord+1) + field_sin(1,iord)+ii*field_sin(2,iord); |
|---|
| 2493 | Cp1 = Cp1 * (x+ii*y) / (iord+2) + field_cos(1,iord)+ii*field_cos(2,iord); |
|---|
| 2494 | Sp1 = Sp1 * (x+ii*y) / (iord+2) + field_sin(1,iord)+ii*field_sin(2,iord); |
|---|
| 2495 | enddo |
|---|
| 2496 | Sp1 = Sp1 * (x+ii*y); |
|---|
| 2497 | Cp1 = Cp1 * (x+ii*y); |
|---|
| 2498 | |
|---|
| 2499 | !---- The kick |
|---|
| 2500 | dpx = -REAL(Cp0) / (one + track(6,jtrk)); |
|---|
| 2501 | dpy = AIMAG(Cp0) / (one + track(6,jtrk)); |
|---|
| 2502 | dpt = (volt * ten3m * sin(lag * 2 * pi - krf * z) / pc - krf * REAL(Sp1)) / (one + track(6,jtrk)); |
|---|
| 2503 | |
|---|
| 2504 | !---- Radiation effects at entrance. |
|---|
| 2505 | if (dorad .and. elrad .ne. zero) then |
|---|
| 2506 | rfac = arad * gammas**3 * (dpx**2+dpy**2) / (three*elrad) |
|---|
| 2507 | track(2,jtrk) = track(2,jtrk) - rfac * (one + track(6,jtrk)) * track(2,jtrk) |
|---|
| 2508 | track(4,jtrk) = track(4,jtrk) - rfac * (one + track(6,jtrk)) * track(4,jtrk) |
|---|
| 2509 | track(6,jtrk) = track(6,jtrk) - rfac * (one + track(6,jtrk)) ** 2 |
|---|
| 2510 | endif |
|---|
| 2511 | |
|---|
| 2512 | !---- Apply the kick |
|---|
| 2513 | track(2,jtrk) = track(2,jtrk) + dpx |
|---|
| 2514 | track(4,jtrk) = track(4,jtrk) + dpy |
|---|
| 2515 | track(6,jtrk) = track(6,jtrk) + dpt |
|---|
| 2516 | |
|---|
| 2517 | !---- Radiation effects at exit. |
|---|
| 2518 | if (dorad .and. elrad .ne. zero) then |
|---|
| 2519 | track(2,jtrk) = track(2,jtrk) - rfac * (one + track(6,jtrk)) * track(2,jtrk) |
|---|
| 2520 | track(4,jtrk) = track(4,jtrk) - rfac * (one + track(6,jtrk)) * track(4,jtrk) |
|---|
| 2521 | track(6,jtrk) = track(6,jtrk) - rfac * (one + track(6,jtrk)) ** 2 |
|---|
| 2522 | endif |
|---|
| 2523 | enddo |
|---|
| 2524 | |
|---|
| 2525 | end subroutine ttrfmult |
|---|
| 2526 | |
|---|
| 2527 | subroutine tttquad(track, ktrack) |
|---|
| 2528 | |
|---|
| 2529 | use twtrrfi |
|---|
| 2530 | use trackfi |
|---|
| 2531 | implicit none |
|---|
| 2532 | |
|---|
| 2533 | !--------------------* |
|---|
| 2534 | ! Andrea Latina 2012 * |
|---|
| 2535 | !--------------------* |
|---|
| 2536 | !----------------------------------------------------------------------* |
|---|
| 2537 | ! Purpose: * |
|---|
| 2538 | ! Track particle through a general thick quadrupole. * |
|---|
| 2539 | ! Input/output: * |
|---|
| 2540 | ! TRACK(6,*)(double) Track coordinates: (X, PX, Y, PY, T, PT). * |
|---|
| 2541 | ! KTRACK (integer) Number of surviving tracks. * |
|---|
| 2542 | !----------------------------------------------------------------------* |
|---|
| 2543 | |
|---|
| 2544 | double precision track(6,*) |
|---|
| 2545 | integer ktrack |
|---|
| 2546 | |
|---|
| 2547 | double precision node_value, get_value |
|---|
| 2548 | double precision k1, k1s, length, ksqrt, tmp |
|---|
| 2549 | double precision sx, cx, sy, cy, ct, st |
|---|
| 2550 | double precision x, px, y, py |
|---|
| 2551 | double precision bet0sqr |
|---|
| 2552 | logical skew, focusing |
|---|
| 2553 | integer jtrk |
|---|
| 2554 | |
|---|
| 2555 | double precision zero, one, two, sqrt2 |
|---|
| 2556 | parameter ( zero=0d0, one=1d0, two=2d0 ) |
|---|
| 2557 | parameter ( sqrt2=1.41421356237310d0 ) |
|---|
| 2558 | |
|---|
| 2559 | !---- Read-in the parameters |
|---|
| 2560 | bet0sqr = bet0*bet0; |
|---|
| 2561 | k1 = node_value('k1 '); |
|---|
| 2562 | k1s = node_value('k1s '); |
|---|
| 2563 | length = node_value('l '); |
|---|
| 2564 | |
|---|
| 2565 | if ((k1.ne.zero).and.(k1s.ne.zero)) then |
|---|
| 2566 | call aawarn('trrun: ',& |
|---|
| 2567 | 'a quadrupole cannot have both K1 and K1S different than zero'); |
|---|
| 2568 | return |
|---|
| 2569 | endif |
|---|
| 2570 | |
|---|
| 2571 | if (k1s.ne.zero) then |
|---|
| 2572 | skew = .true. |
|---|
| 2573 | focusing = k1s.gt.zero |
|---|
| 2574 | else |
|---|
| 2575 | skew = .false. |
|---|
| 2576 | focusing = k1.gt.zero |
|---|
| 2577 | endif |
|---|
| 2578 | |
|---|
| 2579 | !---- Prepare to calculate the kick and the matrix elements |
|---|
| 2580 | do jtrk = 1,ktrack |
|---|
| 2581 | !---- The particle position |
|---|
| 2582 | x = track(1,jtrk); |
|---|
| 2583 | px = track(2,jtrk); |
|---|
| 2584 | y = track(3,jtrk); |
|---|
| 2585 | py = track(4,jtrk); |
|---|
| 2586 | |
|---|
| 2587 | !!$ !---- Radiation effects at entrance. |
|---|
| 2588 | !!$ if (dorad .and. elrad .ne. zero) then |
|---|
| 2589 | !!$ rfac = arad * gammas**3 * (dpx**2+dpy**2) / (three*elrad) |
|---|
| 2590 | !!$ track(2,jtrk) = track(2,jtrk) - rfac * (one + track(6,jtrk)) * track(2,jtrk) |
|---|
| 2591 | !!$ track(4,jtrk) = track(4,jtrk) - rfac * (one + track(6,jtrk)) * track(4,jtrk) |
|---|
| 2592 | !!$ track(6,jtrk) = track(6,jtrk) - rfac * (one + track(6,jtrk)) ** 2 |
|---|
| 2593 | !!$ endif |
|---|
| 2594 | |
|---|
| 2595 | !---- If SKEW rotates by -45 degrees |
|---|
| 2596 | if (skew) then |
|---|
| 2597 | ct = sqrt2 / two |
|---|
| 2598 | st = -sqrt2 / two |
|---|
| 2599 | tmp = x |
|---|
| 2600 | x = ct * tmp + st * y |
|---|
| 2601 | y = ct * y - st * tmp |
|---|
| 2602 | tmp = px |
|---|
| 2603 | px = ct * tmp + st * py |
|---|
| 2604 | py = ct * py - st * tmp |
|---|
| 2605 | endif |
|---|
| 2606 | |
|---|
| 2607 | !---- Computes the kick |
|---|
| 2608 | if (skew) then |
|---|
| 2609 | ksqrt = sqrt(abs(k1s / (one + track(6, jtrk)))) |
|---|
| 2610 | else |
|---|
| 2611 | ksqrt = sqrt(abs(k1 / (one + track(6, jtrk)))) |
|---|
| 2612 | endif |
|---|
| 2613 | if (focusing) then |
|---|
| 2614 | cx = cos(ksqrt*length) |
|---|
| 2615 | sx = sin(ksqrt*length) |
|---|
| 2616 | cy = cosh(ksqrt*length) |
|---|
| 2617 | sy = sinh(ksqrt*length) |
|---|
| 2618 | tmp = x |
|---|
| 2619 | x = tmp * cx + px * sx / ksqrt |
|---|
| 2620 | px = -sx * ksqrt * tmp + px * cx |
|---|
| 2621 | tmp = y |
|---|
| 2622 | y = tmp * cy + py * sy / ksqrt |
|---|
| 2623 | py = sy * ksqrt * tmp + py * cy |
|---|
| 2624 | else |
|---|
| 2625 | cx = cosh(ksqrt*length) |
|---|
| 2626 | sx = sinh(ksqrt*length) |
|---|
| 2627 | cy = cos(ksqrt*length) |
|---|
| 2628 | sy = sin(ksqrt*length) |
|---|
| 2629 | tmp = x |
|---|
| 2630 | x = tmp * cx + px * sx / ksqrt |
|---|
| 2631 | px = sx * ksqrt * tmp + px * cx |
|---|
| 2632 | tmp = y |
|---|
| 2633 | y = tmp * cy + py * sy / ksqrt |
|---|
| 2634 | py = -sy * ksqrt * tmp + py * cy |
|---|
| 2635 | endif |
|---|
| 2636 | |
|---|
| 2637 | !---- If SKEW rotates by +45 degrees |
|---|
| 2638 | if (skew) then |
|---|
| 2639 | ct = sqrt2 / two |
|---|
| 2640 | st = sqrt2 / two |
|---|
| 2641 | tmp = x |
|---|
| 2642 | x = ct * tmp + st * y |
|---|
| 2643 | y = ct * y - st * tmp |
|---|
| 2644 | tmp = px |
|---|
| 2645 | px = ct * tmp + st * py |
|---|
| 2646 | py = ct * py - st * tmp |
|---|
| 2647 | endif |
|---|
| 2648 | |
|---|
| 2649 | !---- Applies the kick |
|---|
| 2650 | track(1,jtrk) = x |
|---|
| 2651 | track(2,jtrk) = px |
|---|
| 2652 | track(3,jtrk) = y |
|---|
| 2653 | track(4,jtrk) = py |
|---|
| 2654 | track(5,jtrk) = track(5,jtrk) + & |
|---|
| 2655 | (one - bet0sqr) * length * track(6,jtrk) / bet0sqr; |
|---|
| 2656 | |
|---|
| 2657 | !!$ !---- Radiation effects at exit. |
|---|
| 2658 | !!$ if (dorad .and. elrad .ne. zero) then |
|---|
| 2659 | !!$ track(2,jtrk) = track(2,jtrk) - rfac * (one + track(6,jtrk)) * track(2,jtrk) |
|---|
| 2660 | !!$ track(4,jtrk) = track(4,jtrk) - rfac * (one + track(6,jtrk)) * track(4,jtrk) |
|---|
| 2661 | !!$ track(6,jtrk) = track(6,jtrk) - rfac * (one + track(6,jtrk)) ** 2 |
|---|
| 2662 | !!$ endif |
|---|
| 2663 | |
|---|
| 2664 | enddo |
|---|
| 2665 | |
|---|
| 2666 | end subroutine tttquad |
|---|
![(please configure the [header_logo] section in trac.ini)](/PSPA/chrome/site/pspa.gif)
