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 |
---|