1 | !----------------------------------------------------------------------* |
---|
2 | subroutine dynap(eigen, coords, turns, npart, distvect, zn,dq, & |
---|
3 | onelog, turnnumber) |
---|
4 | use deltrafi |
---|
5 | implicit none |
---|
6 | |
---|
7 | |
---|
8 | |
---|
9 | integer turns,npart |
---|
10 | double precision eigen(6,6),coords(6,0:turns,*),get_value |
---|
11 | double precision zn(turns,6), distvect(turns),dq(2*turns), & |
---|
12 | onelog(turns), turnnumber(turns) |
---|
13 | |
---|
14 | write(*,*) ' entered dynap ' |
---|
15 | |
---|
16 | ! print *, 'eigenvectors' |
---|
17 | ! do i = 1, 6 |
---|
18 | ! print '(1p,6e12.4)', (eigen(i,j), j=1,6) |
---|
19 | ! enddo |
---|
20 | ! do k = 1, 2 |
---|
21 | ! print *, 'particle: ', k |
---|
22 | ! do i = 1, turns |
---|
23 | ! print '(1p,6e12.4)', (coords(j,i,k), j = 1, 6) |
---|
24 | ! enddo |
---|
25 | ! enddo |
---|
26 | |
---|
27 | fastune = get_value('dynap ','fastune ') .ne. 0 |
---|
28 | deltax = get_value('dynap ','lyapunov ') |
---|
29 | |
---|
30 | call trdynrun (eigen,coords,turns,npart,distvect,zn,onelog, & |
---|
31 | turnnumber,dq) |
---|
32 | ! call dynapfill() |
---|
33 | |
---|
34 | write(*,*) ' end dynap ' |
---|
35 | |
---|
36 | end subroutine dynap |
---|
37 | |
---|
38 | !********************************************************************** |
---|
39 | subroutine dynapfill() |
---|
40 | use dyntabfi |
---|
41 | use wmaxmin0fi |
---|
42 | implicit none |
---|
43 | |
---|
44 | |
---|
45 | !----------------------------------------------------------------------* |
---|
46 | ! Purpose: * |
---|
47 | ! writes the dynap output in the table "dynap" * |
---|
48 | ! Output: * |
---|
49 | ! dynapfrac (real) : fractional dynamic aperture * |
---|
50 | ! dktrturns (real) : number of turns after tracking * |
---|
51 | ! xend (real) : final x position w.r.t. closed orbit * |
---|
52 | ! pxend (real) : final x momentum w.r.t. closed orbit * |
---|
53 | ! yend (real) : final y position w.r.t. closed orbit * |
---|
54 | ! pyend (real) : final y momentum w.r.t. closed orbit * |
---|
55 | ! tend (real) : final longit. position w.r.t. closed orbit * |
---|
56 | ! ptend (real) : final longit. momentum w.r.t. closed orbit * |
---|
57 | ! wxmin (real) : minimum x betatron invariant during tracking * |
---|
58 | ! wxmax (real) : maximum x betatron invariant during tracking * |
---|
59 | ! wymin (real) : minimum y betatron invariant during tracking * |
---|
60 | ! wymax (real) : maximum y betatron invariant during tracking * |
---|
61 | ! wxymin (real) : minimum of (wx + wy) during tracking * |
---|
62 | ! wxymax (real) : maximum of (wx + wy) during tracking * |
---|
63 | ! smear (real) : 2.0 * (wxymax - wxymin) / (wxymax + wxymin) * |
---|
64 | ! yapunov (real) : interpolated Lyapunov exponent * |
---|
65 | !----------------------------------------------------------------------* |
---|
66 | |
---|
67 | call double_to_table_curr('dynap ', 'dynapfrac ', dynapfrac) |
---|
68 | call double_to_table_curr('dynap ', 'dktrturns ', dktrturns) |
---|
69 | call double_to_table_curr('dynap ', 'xend ', xend ) |
---|
70 | call double_to_table_curr('dynap ', 'pxend ', pxend ) |
---|
71 | call double_to_table_curr('dynap ', 'yend ', yend ) |
---|
72 | call double_to_table_curr('dynap ', 'pyend ',pyend ) |
---|
73 | call double_to_table_curr('dynap ', 'tend ', ptend ) |
---|
74 | call double_to_table_curr('dynap ', 'wxmin ', wxmin ) |
---|
75 | call double_to_table_curr('dynap ', 'wxmax ', wxmax ) |
---|
76 | call double_to_table_curr('dynap ', 'wymin ', wymin ) |
---|
77 | call double_to_table_curr('dynap ', 'wymax ', wymax ) |
---|
78 | call double_to_table_curr('dynap ', 'wxymin ', wxymin ) |
---|
79 | call double_to_table_curr('dynap ', 'wxymax ', wxymax ) |
---|
80 | call double_to_table_curr('dynap ', 'smear ', smear ) |
---|
81 | call double_to_table_curr('dynap ', 'yapunov ', yapunov ) |
---|
82 | |
---|
83 | write(*,*) ' yapunov = ', yapunov |
---|
84 | |
---|
85 | call augment_count('dynap ') |
---|
86 | end subroutine dynapfill |
---|
87 | !----------------- end of dynapfill subroutine -------------------------- |
---|
88 | !********************************************************************** |
---|
89 | subroutine dynaptunefill |
---|
90 | use tunesfi |
---|
91 | implicit none |
---|
92 | |
---|
93 | |
---|
94 | !----------------------------------------------------------------------* |
---|
95 | ! Purpose: * |
---|
96 | ! writes the dynap tunes in the table "dynaptune" * |
---|
97 | ! Output: * |
---|
98 | ! tunx (real) : x tune from FFT * |
---|
99 | ! tuny (real) : y tune from FFT * |
---|
100 | ! dtune (real) : tune error * |
---|
101 | !----------------------------------------------------------------------* |
---|
102 | double precision tx_tmp, ty_tmp |
---|
103 | tx_tmp = tunx |
---|
104 | if (tunx .gt. 0.5d0) tx_tmp = 1.d0 - tunx |
---|
105 | ty_tmp = tuny |
---|
106 | if (tuny .gt. 0.5d0) ty_tmp = 1.d0 - tuny |
---|
107 | call double_to_table_curr('dynaptune ', 'x ' , x0) |
---|
108 | call double_to_table_curr('dynaptune ', 'y ' , y0) |
---|
109 | call double_to_table_curr('dynaptune ', 'tunx ', tx_tmp) |
---|
110 | call double_to_table_curr('dynaptune ', 'tuny ', ty_tmp) |
---|
111 | call double_to_table_curr('dynaptune ', 'dtune ', dtune ) |
---|
112 | |
---|
113 | write(*,*) ' tunes = ', tunx, tuny, dtune |
---|
114 | |
---|
115 | call augment_count('dynaptune ') |
---|
116 | end subroutine dynaptunefill |
---|
117 | !----------------- end of dynaptunefill subroutine -------------------------- |
---|
118 | |
---|
119 | |
---|
120 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
---|
121 | subroutine trdynrun (eigen,coords,turns,npart,distvect,zn,onelog, & |
---|
122 | turnnumber,dq) |
---|
123 | |
---|
124 | use deltrafi |
---|
125 | use wmaxmin0fi |
---|
126 | use dyntabfi |
---|
127 | use tunesfi |
---|
128 | implicit none |
---|
129 | |
---|
130 | !----------------------------------------------------------------------* |
---|
131 | ! Purpose: * |
---|
132 | ! |
---|
133 | !----------------------------------------------------------------------* |
---|
134 | integer k,ix,iy,initt,turns,ktrturns,nturnhalf,i,j,npart |
---|
135 | double precision eigen(6,6),coords(6,0:turns,*),distvect(turns), & |
---|
136 | zn(turns,6),znt(6),track(6),fitlyap, & |
---|
137 | templyap,tunx1,tunx2,tuny1,tuny2,tuneabt,tuneabt2,zero,one,two, & |
---|
138 | onelog(turns), turnnumber(turns),dq(2*turns),dphi(3),dphi1,dphi2, & |
---|
139 | get_variable,pi,twopi |
---|
140 | parameter(zero=0d0,one=1d0,two=2d0) |
---|
141 | |
---|
142 | pi=get_variable('pi ') |
---|
143 | twopi=get_variable('twopi ') |
---|
144 | dktrturns = turns |
---|
145 | |
---|
146 | !---- Initialize max and min betatron invariants. |
---|
147 | wxmax = zero |
---|
148 | wymax = zero |
---|
149 | wxymax = zero |
---|
150 | wxmin = zero |
---|
151 | wymin = zero |
---|
152 | wxymin = zero |
---|
153 | |
---|
154 | !--- distance for second particle with x add-on |
---|
155 | ! deltax = get_value('dynap ', 'lyapunov ') |
---|
156 | ! jend = 2 |
---|
157 | ! z(1,jend) = z(1,1) + deltax |
---|
158 | ! do k = 2, 6 |
---|
159 | ! z(k,jend) = z(k,1) |
---|
160 | ! enddo |
---|
161 | |
---|
162 | !---- compute minimum and maximum amplitudes |
---|
163 | ! and normalized coordinates |
---|
164 | ! for all particles |
---|
165 | do k = 1, npart, 2 |
---|
166 | x0 = coords(1,0,k) |
---|
167 | y0 = coords(3,0,k) |
---|
168 | do i = 1, turns |
---|
169 | do j = 1, 6 |
---|
170 | track(j) = coords(j,i,k) |
---|
171 | ! write(*,*) track(j) |
---|
172 | end do |
---|
173 | call wmaxmin(track,eigen,znt) |
---|
174 | do j = 1, 6 |
---|
175 | zn(i,j)=znt(j) |
---|
176 | ! write(*,*) zn(i,j) |
---|
177 | end do |
---|
178 | end do |
---|
179 | |
---|
180 | !---- compute 'smear' |
---|
181 | smear = two * (wxymax - wxymin) / (wxymax + wxymin) |
---|
182 | |
---|
183 | ! write(*,*)' smear = ', smear |
---|
184 | |
---|
185 | !---- Fast tune calculation by interpolated FFT. |
---|
186 | if (fastune) then |
---|
187 | |
---|
188 | ktrturns = turns |
---|
189 | |
---|
190 | ix = 1 |
---|
191 | iy =3 |
---|
192 | initt = 0 |
---|
193 | if (ktrturns .le. 64) then |
---|
194 | tunx = tuneabt(zn, ix, initt, ktrturns, turns, dq) |
---|
195 | tuny = tuneabt(zn, iy, initt, ktrturns, turns, dq) |
---|
196 | else |
---|
197 | tunx = tuneabt2(zn, ix, initt, ktrturns, turns, dq) |
---|
198 | tuny = tuneabt2(zn, iy, initt, ktrturns, turns, dq) |
---|
199 | endif |
---|
200 | |
---|
201 | !---- Fast tune variation over half the number of turns. |
---|
202 | nturnhalf = ktrturns / 2 |
---|
203 | if (nturnhalf .le. 64) then |
---|
204 | initt = 0 |
---|
205 | tunx1 = tuneabt(zn, ix, initt, nturnhalf, turns, dq) |
---|
206 | tuny1 = tuneabt(zn, iy, initt, nturnhalf, turns, dq) |
---|
207 | initt = nturnhalf |
---|
208 | tunx2 = tuneabt(zn, ix, initt, nturnhalf, turns, dq) |
---|
209 | tuny2 = tuneabt(zn, iy, initt, nturnhalf, turns, dq) |
---|
210 | else |
---|
211 | initt = 0 |
---|
212 | tunx1 = tuneabt2(zn, ix, initt, nturnhalf, turns, dq) |
---|
213 | tuny1 = tuneabt2(zn, iy, initt, nturnhalf, turns, dq) |
---|
214 | initt = nturnhalf |
---|
215 | tunx2 = tuneabt2(zn, ix, initt, nturnhalf, turns, dq) |
---|
216 | tuny2 = tuneabt2(zn, iy, initt, nturnhalf, turns, dq) |
---|
217 | endif |
---|
218 | if (abs(tunx1-tunx).gt.0.4) tunx1 = 1-tunx1 |
---|
219 | if (abs(tunx2-tunx).gt.0.4) tunx2 = 1-tunx2 |
---|
220 | if (abs(tuny1-tuny).gt.0.4) tuny1 = 1-tuny1 |
---|
221 | if (abs(tuny2-tuny).gt.0.4) tuny2 = 1-tuny2 |
---|
222 | dtune = sqrt((tunx2 - tunx1)**2 + (tuny2 - tuny1)**2) |
---|
223 | write(*,*) tunx1,tunx2, ' ',tuny1,tuny2 |
---|
224 | call dynaptunefill |
---|
225 | endif |
---|
226 | |
---|
227 | !---- Lyapunov exponent calculation |
---|
228 | |
---|
229 | open(50,file="lyapunov.data") |
---|
230 | |
---|
231 | do i = 1, turns |
---|
232 | do j = 1, 6 |
---|
233 | track(j) = coords(j,i,k+1) |
---|
234 | ! write(*,*) track(j) |
---|
235 | end do |
---|
236 | call wmaxmin(track,eigen,znt) |
---|
237 | templyap = zero |
---|
238 | do j = 1, 3 |
---|
239 | if(znt(2*j).eq.zero.and.znt(2*j-1).eq.zero) then |
---|
240 | dphi1=zero |
---|
241 | else |
---|
242 | dphi1=atan2(znt(2*j),znt(2*j-1)) |
---|
243 | endif |
---|
244 | if(zn(i,2*j).eq.zero.and.zn(i,2*j-1).eq.zero) then |
---|
245 | dphi2=zero |
---|
246 | else |
---|
247 | dphi2=atan2(zn(i,2*j),zn(i,2*j-1)) |
---|
248 | endif |
---|
249 | dphi(j) = dphi1-dphi2 |
---|
250 | if (dphi(j).gt.pi) dphi(j)=dphi(j)-twopi |
---|
251 | if (dphi(j).lt.-pi) dphi(j)=dphi(j)+twopi |
---|
252 | templyap = templyap + (dphi(j))**2 |
---|
253 | end do |
---|
254 | distvect(i) = sqrt(templyap) |
---|
255 | write(50,*) i,distvect(i) |
---|
256 | ! write(*,*) distvect(i) |
---|
257 | end do |
---|
258 | ! write(*,*) ' templyap =', templyap |
---|
259 | yapunov = fitlyap(distvect,onelog,turnnumber,turns) |
---|
260 | call dynapfill() |
---|
261 | |
---|
262 | enddo |
---|
263 | end subroutine trdynrun |
---|
264 | |
---|
265 | !----------------- end of trdynrun subroutine -------------------------- |
---|
266 | |
---|
267 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
---|
268 | subroutine fft(data, nn, isign) |
---|
269 | |
---|
270 | implicit none |
---|
271 | !----------------------------------------------------------------------* |
---|
272 | ! Purpose: * |
---|
273 | ! Computes the FFT * |
---|
274 | ! Author: numerical receipes, pg. 395 * |
---|
275 | ! DATA: is a real array with the signal on input * |
---|
276 | ! with the Fourier transform on output. * |
---|
277 | ! N: is the number of data: must be a power of 2 * |
---|
278 | ! ISIGN=1: direct Fourier transform * |
---|
279 | ! ISIGN=-1: inverse Fourier transform * |
---|
280 | !----------------------------------------------------------------------* |
---|
281 | ! frankz: only removed a set of parameter statements - keep rest |
---|
282 | ! unchanged |
---|
283 | ! |
---|
284 | !---- Double precision version. |
---|
285 | integer i,isign,istep,j,m,mmax,n,nn |
---|
286 | double precision data(*),tempi,tempr,theta,wi,wpi,wpr,wr,wtemp, & |
---|
287 | get_variable,zero,half,one,two,twopi |
---|
288 | parameter(zero=0d0,half=0.5d0,one=1d0,two=2d0) |
---|
289 | |
---|
290 | !---- Initialize |
---|
291 | twopi=get_variable('twopi ') |
---|
292 | |
---|
293 | !---- Rearrange the data points. |
---|
294 | n = 2 * nn |
---|
295 | j = 1 |
---|
296 | do i = 1, n, 2 |
---|
297 | if(j .gt. i) then |
---|
298 | tempr = data(j) |
---|
299 | tempi = data(j+1) |
---|
300 | data(j) = data(i) |
---|
301 | data(j+1) = data(i+1) |
---|
302 | data(i) = tempr |
---|
303 | data(i+1) = tempi |
---|
304 | endif |
---|
305 | m = n / 2 |
---|
306 | 1 if (m .ge. 2 .and. j .gt. m) then |
---|
307 | j = j - m |
---|
308 | m = m / 2 |
---|
309 | go to 1 |
---|
310 | endif |
---|
311 | j = j + m |
---|
312 | enddo |
---|
313 | mmax = 2 |
---|
314 | 2 if (n .gt. mmax) then |
---|
315 | istep = 2 * mmax |
---|
316 | theta = twopi / (isign * mmax) |
---|
317 | wpr = - two * sin(half * theta)**2 |
---|
318 | wpi = sin(theta) |
---|
319 | wr = one |
---|
320 | wi = zero |
---|
321 | do m = 1, mmax, 2 |
---|
322 | do i = m, n, istep |
---|
323 | j = i + mmax |
---|
324 | tempr = wr * data(j) - wi * data(j+1) |
---|
325 | tempi = wr * data(j+1) + wi * data(j) |
---|
326 | data(j) = data(i) - tempr |
---|
327 | data(j+1) = data(i+1) - tempi |
---|
328 | data(i) = data(i) + tempr |
---|
329 | data(i+1) = data(i+1) + tempi |
---|
330 | enddo |
---|
331 | wtemp = wr |
---|
332 | wr = wr * wpr - wi * wpi + wr |
---|
333 | wi = wi * wpr + wtemp * wpi + wi |
---|
334 | enddo |
---|
335 | mmax = istep |
---|
336 | go to 2 |
---|
337 | endif |
---|
338 | |
---|
339 | end subroutine fft |
---|
340 | |
---|
341 | !----------------- end of fft subroutine -------------------------- |
---|
342 | |
---|
343 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
---|
344 | double precision function tuneabt(zn, ixy, initt, maxn, turns, dq) |
---|
345 | |
---|
346 | implicit none |
---|
347 | |
---|
348 | !----------------------------------------------------------------------* |
---|
349 | ! Purpose: * |
---|
350 | ! Computes the tune using formula (18) of CERN SL/95-84 (AP). * |
---|
351 | ! No filter. * |
---|
352 | ! (best suited for MAXN <= 64 TURNS) * |
---|
353 | ! zn(*,4) contains the 4 transverse coordinates for different turns * |
---|
354 | ! maxn is the number of turns * |
---|
355 | ! Authors: * |
---|
356 | ! R. Bartolini - CERN and Bologna University, * |
---|
357 | ! E. TODESCO - INFN and CERN. * |
---|
358 | !----------------------------------------------------------------------* |
---|
359 | integer mft,nft,ixy,initt,maxn,nftmax,npoint,mf,turns |
---|
360 | ! choose dimensions for now |
---|
361 | double precision zn(turns,6),dq(2*turns),ftmax,temp,cf1,cf2,cf3, & |
---|
362 | arg,assk,pi,get_variable,zero,one,two |
---|
363 | parameter(zero=0d0,one=1d0,two=2d0) |
---|
364 | |
---|
365 | !---- Initialize |
---|
366 | pi=get_variable('pi ') |
---|
367 | |
---|
368 | !---- Use first NPOINT points. |
---|
369 | mft = int(log(float(maxn)) / log(two)) |
---|
370 | npoint = 2**mft |
---|
371 | |
---|
372 | !---- Copy data to dq |
---|
373 | do mf = 1, npoint |
---|
374 | dq(2*mf-1) = zn(mf+initt,ixy) |
---|
375 | dq(2*mf) = zn(mf+initt,ixy+1) |
---|
376 | enddo |
---|
377 | call fft(dq, npoint, -1) |
---|
378 | |
---|
379 | !---- Search for maximum of Fourier spectrum. |
---|
380 | ftmax = zero |
---|
381 | nftmax = 0 |
---|
382 | do nft = 1, npoint |
---|
383 | temp = sqrt(dq(2*nft-1)**2 + dq(2*nft)**2) |
---|
384 | if (temp .gt. ftmax) then |
---|
385 | ftmax = temp |
---|
386 | nftmax = nft |
---|
387 | endif |
---|
388 | enddo |
---|
389 | |
---|
390 | !---- Improve estimate by interpolation. |
---|
391 | cf1 = sqrt(dq(2*nftmax-3)**2 + dq(2*nftmax-2)**2) |
---|
392 | cf2 = sqrt(dq(2*nftmax-1)**2 + dq(2*nftmax)**2) |
---|
393 | cf3 = sqrt(dq(2*nftmax+1)**2 + dq(2*nftmax+2)**2) |
---|
394 | if (cf3 .gt. cf1) then |
---|
395 | arg = sin(pi / npoint) / (cf2 / cf3 + cos(pi / npoint)) |
---|
396 | assk = float(nftmax) + npoint / pi * atan(arg) |
---|
397 | else |
---|
398 | arg = sin(pi / npoint) / (cf1 / cf2 + cos(pi / npoint)) |
---|
399 | assk = float(nftmax-1) + npoint / pi * atan(arg) |
---|
400 | endif |
---|
401 | tuneabt = one - (assk - one) / float(npoint) |
---|
402 | |
---|
403 | end function tuneabt |
---|
404 | !----------------- end of tuneabt function -------------------------- |
---|
405 | |
---|
406 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
---|
407 | double precision function tuneabt2(zn, ixy, initt, maxn, turns,dq) |
---|
408 | |
---|
409 | implicit none |
---|
410 | |
---|
411 | !----------------------------------------------------------------------* |
---|
412 | ! Purpose: * |
---|
413 | ! Computes the tune using the interpolated FFT with Hanning filter. * |
---|
414 | ! See CERN SL/95-84 formula (25). * |
---|
415 | ! (best suited for MAXN > 64 turns) * |
---|
416 | ! X, XP are the coordinates of the orbit, * |
---|
417 | ! MAXN is the length of the orbit. * |
---|
418 | ! Authors: * |
---|
419 | ! R. Bartolini - CERN and Bologna University, * |
---|
420 | ! E. TODESCO - INFN and CERN. * |
---|
421 | !----------------------------------------------------------------------* |
---|
422 | integer mft,nft,mf,ixy,initt,maxn,nftmax,npoint,nn,turns |
---|
423 | ! need dimensions? |
---|
424 | double precision zn(turns,6),dq(2*turns),ftmax,temp,step,cf1,cf2, & |
---|
425 | scra1,scra2,scra3,scra4,assk,co,si,p1,p2,pi,twopi,get_variable, & |
---|
426 | zero,one,two,cf3 |
---|
427 | parameter(zero=0d0,one=1d0,two=2d0) |
---|
428 | |
---|
429 | !---- Initialize |
---|
430 | pi=get_variable('pi ') |
---|
431 | twopi=get_variable('twopi ') |
---|
432 | |
---|
433 | !---- Use first NPOINT points. |
---|
434 | mft = int(log(float(maxn)) / log(two)) |
---|
435 | npoint = 2**mft |
---|
436 | |
---|
437 | !---- Copy data to local storage using Hanning filter. |
---|
438 | step = pi / npoint |
---|
439 | do mf = 1, npoint |
---|
440 | temp = sin(mf * step)**2 |
---|
441 | dq(2*mf-1) = temp * zn(mf+initt,ixy) |
---|
442 | dq(2*mf) = temp * zn(mf+initt,ixy+1) |
---|
443 | enddo |
---|
444 | call fft(dq, npoint, -1) |
---|
445 | |
---|
446 | !---- Search for maximum of Fourier spectrum. |
---|
447 | ftmax = zero |
---|
448 | nftmax = 0 |
---|
449 | do nft = 1, npoint |
---|
450 | temp = sqrt(dq(2*nft-1)**2 + dq(2*nft)**2) |
---|
451 | if (temp .gt. ftmax) then |
---|
452 | ftmax = temp |
---|
453 | nftmax = nft |
---|
454 | endif |
---|
455 | enddo |
---|
456 | cf1 = sqrt(dq(2*nftmax-3)**2 + dq(2*nftmax-2)**2) |
---|
457 | cf2 = sqrt(dq(2*nftmax-1)**2 + dq(2*nftmax)**2) |
---|
458 | cf3 = sqrt(dq(2*nftmax+1)**2 + dq(2*nftmax+2)**2) |
---|
459 | if (cf3 .gt. cf1) then |
---|
460 | p1 = cf2 |
---|
461 | p2 = cf3 |
---|
462 | nn = nftmax |
---|
463 | else |
---|
464 | p1 = cf1 |
---|
465 | p2 = cf2 |
---|
466 | nn = nftmax-1 |
---|
467 | endif |
---|
468 | |
---|
469 | !---- Interpolation. |
---|
470 | co = cos(twopi / npoint) |
---|
471 | si = sin(twopi / npoint) |
---|
472 | scra1 = co**2 * (p1 + p2)**2 - 2*p1*p2*(2*co**2 - co - one) |
---|
473 | scra2 = (p1 + p2*co) * (p1 - p2) |
---|
474 | scra3 = p1**2 + p2**2 + 2*p1*p2*co |
---|
475 | scra4 = (-scra2 + p2*sqrt(scra1)) / scra3 |
---|
476 | assk = nn + npoint / twopi * asin(si * scra4) |
---|
477 | tuneabt2 = one - (assk - one) / float(npoint) |
---|
478 | |
---|
479 | end function tuneabt2 |
---|
480 | !----------------- end of tuneabt2 function -------------------------- |
---|
481 | |
---|
482 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
---|
483 | double precision function fitlyap(distvect, onelog, turnnumber, & |
---|
484 | nturn) |
---|
485 | implicit none |
---|
486 | |
---|
487 | |
---|
488 | !----------------------------------------------------------------------* |
---|
489 | ! Purpose: |
---|
490 | ! Computes interpolated Lyapunov exponent. |
---|
491 | ! DISTVECT (normalized) distance between two companion particles |
---|
492 | ! NTURN is the number of turns. |
---|
493 | !----------------------------------------------------------------------* |
---|
494 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
---|
495 | ! needs: |
---|
496 | ! slopexy |
---|
497 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
---|
498 | integer mf,n1,n2,n3,npoint,nturn |
---|
499 | double precision onelog(*),turnnumber(*),fitlyap1,fitlyap2, & |
---|
500 | fitlyap3,deltalog1,deltalog2,deltalog3,distvect(*),slopexy,zero, & |
---|
501 | one,dlmax |
---|
502 | parameter(zero=0d0,one=1d0,dlmax=1d-5) |
---|
503 | |
---|
504 | do mf = 1, nturn |
---|
505 | ! dq(ilogd+mf) = log(distvect(mf)) |
---|
506 | ! dq(in+mf) = mf |
---|
507 | ! dq(ilogn+mf) = log(dq(in+mf)) |
---|
508 | if(distvect(mf).eq.zero) then |
---|
509 | onelog(mf) = zero |
---|
510 | else |
---|
511 | onelog(mf) = log(distvect(mf)) |
---|
512 | endif |
---|
513 | turnnumber(mf) = real(mf) |
---|
514 | ! turnlog(mf) = log(dble(turnnumber(mf))) |
---|
515 | enddo |
---|
516 | |
---|
517 | !---- Loglog fit over 3 subsequent periods of NPOINT = NTURN/4 TURNS |
---|
518 | ! starting at N1 = NPOINT + 1, i.e. at the second fourth |
---|
519 | npoint = int(nturn / 4) |
---|
520 | n1 = npoint + 1 |
---|
521 | n2 = n1 + npoint |
---|
522 | n3 = n2 + npoint |
---|
523 | |
---|
524 | !---- DELTALOG = deviation from 1 of loglog slope. |
---|
525 | ! |
---|
526 | deltalog1 = slopexy(turnnumber(n1), onelog(n1), npoint) - one |
---|
527 | deltalog2 = slopexy(turnnumber(n2), onelog(n2), npoint) - one |
---|
528 | deltalog3 = slopexy(turnnumber(n3), onelog(n3), npoint) - one |
---|
529 | |
---|
530 | if (deltalog1 .lt. dlmax .and. deltalog2 .lt. dlmax .and. & |
---|
531 | deltalog3 .lt. dlmax) then |
---|
532 | fitlyap = zero |
---|
533 | else |
---|
534 | fitlyap1 = slopexy(turnnumber(n1), onelog(n1), npoint) |
---|
535 | fitlyap2 = slopexy(turnnumber(n2), onelog(n2), npoint) |
---|
536 | fitlyap3 = slopexy(turnnumber(n3), onelog(n3), npoint) |
---|
537 | |
---|
538 | if (fitlyap1 .lt. fitlyap2) then |
---|
539 | fitlyap = fitlyap2 |
---|
540 | else |
---|
541 | fitlyap = fitlyap1 |
---|
542 | endif |
---|
543 | if (fitlyap .lt. fitlyap3) then |
---|
544 | fitlyap = fitlyap3 |
---|
545 | endif |
---|
546 | endif |
---|
547 | |
---|
548 | end function fitlyap |
---|
549 | !----------------- end of fitlyap function -------------------------- |
---|
550 | |
---|
551 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
---|
552 | subroutine wmaxmin(track,eigen,znt) |
---|
553 | |
---|
554 | use wmaxmin0fi |
---|
555 | implicit none |
---|
556 | |
---|
557 | !----------------------------------------------------------------------* |
---|
558 | ! Purpose: * |
---|
559 | ! Computes maximum and minimum betatron invariants during tracking. * |
---|
560 | ! Input: * |
---|
561 | ! TRACK(6,*)(real) Track coordinates: (X, PX, Y, PY, T, PT). * |
---|
562 | !----------------------------------------------------------------------* |
---|
563 | integer kp,kq |
---|
564 | double precision track(6),wx,wxy,wy,znt(6),eigen(6,6),zero |
---|
565 | parameter(zero=0d0) |
---|
566 | |
---|
567 | !---- Copy track coordinates. |
---|
568 | ! do 10 i = 1, 6 |
---|
569 | ! z(i) = track(i) |
---|
570 | ! 10 continue |
---|
571 | |
---|
572 | !---- Convert to normalized values. |
---|
573 | do kq = 1, 5, 2 |
---|
574 | kp = kq + 1 |
---|
575 | znt(kq) = eigen(2,kp) * track(1) - eigen(1,kp) * track(2) & |
---|
576 | + eigen(4,kp) * track(3) - eigen(3,kp) * track(4) & |
---|
577 | + eigen(6,kp) * track(5) - eigen(5,kp) * track(6) |
---|
578 | znt(kp) = eigen(1,kq) * track(2) - eigen(2,kq) * track(1) & |
---|
579 | + eigen(3,kq) * track(4) - eigen(4,kq) * track(3) & |
---|
580 | + eigen(5,kq) * track(6) - eigen(6,kq) * track(5) |
---|
581 | enddo |
---|
582 | |
---|
583 | !---- Convert to amplitudes (and phases: not computed). |
---|
584 | wx = znt(1)**2 + znt(2)**2 |
---|
585 | wy = znt(3)**2 + znt(4)**2 |
---|
586 | wxy = wx + wy |
---|
587 | |
---|
588 | !---- Compare to and redefine WMIN and WMAX in TRDYNAP. |
---|
589 | if (wx.gt.wxmax) then |
---|
590 | wxmax = wx |
---|
591 | else if (wx.lt.wxmin.or.wxmin.eq.zero) then |
---|
592 | wxmin = wx |
---|
593 | endif |
---|
594 | if (wy.gt.wymax) then |
---|
595 | wymax = wy |
---|
596 | else if (wy.lt.wymin.or.wymin.eq.zero) then |
---|
597 | wymin = wy |
---|
598 | endif |
---|
599 | if (wxy.gt.wxymax) then |
---|
600 | wxymax = wxy |
---|
601 | else if (wxy.lt.wxymin.or.wxymin.eq.zero) then |
---|
602 | wxymin = wxy |
---|
603 | endif |
---|
604 | |
---|
605 | end subroutine wmaxmin |
---|
606 | !----------------- end of wmaxmin subroutine -------------------------- |
---|
607 | |
---|
608 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
---|
609 | double precision function slopexy(vectorx, vectory, nturn) |
---|
610 | implicit none |
---|
611 | !----------------------------------------------------------------------* |
---|
612 | ! Purpose: * |
---|
613 | ! Computes slope a for linear fit of Y = a X + b. * |
---|
614 | ! VECTORX is the array of abscissas X, * |
---|
615 | ! VECTORX is the array of ordinates Y to interpolate, * |
---|
616 | ! NTURN is their dimension. * |
---|
617 | !----------------------------------------------------------------------* |
---|
618 | integer mf,nturn |
---|
619 | double precision vectorx(*),vectory(*),x2mean,xmean,xymean,ymean, & |
---|
620 | zero |
---|
621 | parameter(zero=0d0) |
---|
622 | |
---|
623 | xmean = zero |
---|
624 | ymean = zero |
---|
625 | x2mean = zero |
---|
626 | xymean = zero |
---|
627 | |
---|
628 | do mf = 1, nturn |
---|
629 | xmean = xmean + vectorx(mf) |
---|
630 | ymean = ymean + vectory(mf) |
---|
631 | x2mean = x2mean + vectorx(mf)**2 |
---|
632 | xymean = xymean + vectorx(mf) * vectory(mf) |
---|
633 | enddo |
---|
634 | |
---|
635 | xmean = xmean / nturn |
---|
636 | ymean = ymean / nturn |
---|
637 | x2mean = x2mean / nturn |
---|
638 | xymean = xymean / nturn |
---|
639 | |
---|
640 | slopexy = (xymean - xmean * ymean) / (x2mean - xmean**2) |
---|
641 | |
---|
642 | end function slopexy |
---|
643 | !----------------- end of slopexy function -------------------------- |
---|