| 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 -------------------------- |
|---|
![(please configure the [header_logo] section in trac.ini)](/PSPA/chrome/site/pspa.gif)
