[430] | 1 | subroutine soddin(ierr) |
---|
| 2 | !--------------------------------------------------------------------- |
---|
| 3 | !--------------------------------------------------------------------- |
---|
| 4 | ! The name SODD stands for "Second Order Detuning and Distortion" |
---|
| 5 | ! ---- - - - - |
---|
| 6 | !--------------------------------------------------------------------- |
---|
| 7 | !--------------------------------------------------------------------- |
---|
| 8 | ! |
---|
| 9 | ! Programm SODD consists of 3 parts: |
---|
| 10 | ! |
---|
| 11 | ! 1.) Subroutine Detune calculates the detuning function terms in |
---|
| 12 | ! first and second order in the strength of the multipoles. |
---|
| 13 | ! 2.) Subroutine Distort1 calculates the distortion function and |
---|
| 14 | ! Hamiltonian terms in first order in the strength of the |
---|
| 15 | ! multipoles. |
---|
| 16 | ! 3.) Subroutine Distort2 calculates the distortion function and |
---|
| 17 | ! Hamiltonian terms in second order in the strength of the |
---|
| 18 | ! multipoles. |
---|
| 19 | ! |
---|
| 20 | ! |
---|
| 21 | ! Theory: Analytical Calculation of Smear and Tune Shift SSC-232 |
---|
| 22 | ! Authors: J. Bengtsson and J. Irwin |
---|
| 23 | ! Date: February 1990 |
---|
| 24 | ! |
---|
| 25 | ! Program Author: Frank Schmidt, CERN SL-Group |
---|
| 26 | ! Date: November 1998 - January 1999 |
---|
| 27 | ! Maximum Order: parameter mmul=11 (22-Pole) |
---|
| 28 | ! Comments to: Frank.Schmidt@cern.ch |
---|
| 29 | ! |
---|
| 30 | !--------------------------------------------------------------------- |
---|
| 31 | ! |
---|
| 32 | ! Requirements: |
---|
| 33 | ! ------------- |
---|
| 34 | ! |
---|
| 35 | ! The program expects a file "fc.34" which contains 8 columns: |
---|
| 36 | ! - Position of Multipole [m] |
---|
| 37 | ! - Name of Multipole |
---|
| 38 | ! - Multipole Type (>0 erect Multipole; <0 skew Multipole); |
---|
| 39 | ! eg "3" is an erect sextupole, "-5" is a skew decapole |
---|
| 40 | ! - Single particle strength `a la SIXTRACK |
---|
| 41 | ! - Horizontal Betafunction |
---|
| 42 | ! - Vertical Betafunction |
---|
| 43 | ! - Horizontal Phase-advance |
---|
| 44 | ! - Vertical Phase-advance |
---|
| 45 | ! |
---|
| 46 | ! The last line in this file is special: All entries are given at |
---|
| 47 | ! the end of the accelerator with the name "END" and the artificial |
---|
| 48 | ! multipole type "100". |
---|
| 49 | ! |
---|
| 50 | ! This file can be automatically produced with SIXTRACK when the |
---|
| 51 | ! linear optics parameter are calculated with the "LINE" block |
---|
| 52 | ! (see web locaction: http://wwwslap.cern.ch/~frs/). It |
---|
| 53 | ! is also produced directly from MAD with the DOOM program by |
---|
| 54 | ! invoking the MAD-SIXTRACK convertor using the `special' flag (see |
---|
| 55 | ! web locaction: http://wwwslap.cern.ch/~hansg/doom/doom_six.html |
---|
| 56 | ! for details). |
---|
| 57 | ! |
---|
| 58 | ! The program needs roughly 50Mbytes in particular for the second |
---|
| 59 | ! order distortion function due to pairs of multipoles up to 11th |
---|
| 60 | ! order. |
---|
| 61 | ! |
---|
| 62 | !--------------------------------------------------------------------- |
---|
| 63 | ! |
---|
| 64 | ! Program Manual: |
---|
| 65 | ! --------------- |
---|
| 66 | ! |
---|
| 67 | ! a.) Program selection: |
---|
| 68 | ! Specify the program: |
---|
| 69 | ! - iprog=1,3,5,7 => Run Program Detune |
---|
| 70 | ! - iprog=2,3,6,7 => Run Program Distort1 |
---|
| 71 | ! - iprog=4,5,6,7 => Run Program Distort2 |
---|
| 72 | ! b.) Order selection: |
---|
| 73 | ! Specify low (n1) and high (n2) limit of orders to be studied |
---|
| 74 | ! erect and skew elements are denoted with positive and negative |
---|
| 75 | ! values respectively; eg.: |
---|
| 76 | ! n1=-4, n2=10 => the orders (-4,-3,-2 and 3-10) will be treated |
---|
| 77 | ! c.) Position Range: |
---|
| 78 | ! Specify a position range (etl1,etl2) in [m], the multipoles |
---|
| 79 | ! located between this range will be used for the calculation |
---|
| 80 | ! d.) Print_out Switch: |
---|
| 81 | ! Output files meant for spreadsheets, i.e. no comments. Human |
---|
| 82 | ! readable output is found in fort.6 |
---|
| 83 | ! iu_on = 0 => No printout |
---|
| 84 | ! iu_on = 1 => Output at the end of the position range |
---|
| 85 | ! iu_on = 2 => At each Multipole in the position range |
---|
| 86 | ! In Detail: |
---|
| 87 | ! |
---|
| 88 | ! Program Detune |
---|
| 89 | ! |
---|
| 90 | ! Multipole |
---|
| 91 | ! Strength |
---|
| 92 | ! Unit iu_on Order Contents in each column |
---|
| 93 | ! ---------------------------------------------------------------- |
---|
| 94 | ! 70 1 1 multipole order |
---|
| 95 | ! (hor.,ver.) plane => (1,2) |
---|
| 96 | ! hor. or ver. detuning |
---|
| 97 | ! order of horizontal Emittance |
---|
| 98 | ! order of vertical Emittance |
---|
| 99 | ! ---------------------------------------------------------------- |
---|
| 100 | ! 71 2 1 multipole order |
---|
| 101 | ! (hor.,ver.) plane => (1,2) |
---|
| 102 | ! hor. or ver. detuning |
---|
| 103 | ! order of horizontal Emittance |
---|
| 104 | ! order of vertical Emittance |
---|
| 105 | ! ---------------------------------------------------------------- |
---|
| 106 | ! 72 1 2 first multipole order |
---|
| 107 | ! second multipole order |
---|
| 108 | ! horizontal detuning |
---|
| 109 | ! order of horizontal Emittance |
---|
| 110 | ! order of vertical Emittance |
---|
| 111 | ! ---------------------------------------------------------------- |
---|
| 112 | ! 73 1 2 first multipole order |
---|
| 113 | ! second multipole order |
---|
| 114 | ! vertical detuning |
---|
| 115 | ! order of horizontal Emittance |
---|
| 116 | ! order of vertical Emittance |
---|
| 117 | ! --------------------------------------------------------------- |
---|
| 118 | ! |
---|
| 119 | ! Program Distort1 |
---|
| 120 | ! |
---|
| 121 | ! Multipole |
---|
| 122 | ! Strength |
---|
| 123 | ! Unit iu_on Order Contents in each column |
---|
| 124 | ! ---------------------------------------------------------------- |
---|
| 125 | ! 74 1 1 multipole order |
---|
| 126 | ! cosine part of distortion |
---|
| 127 | ! sine part of distortion |
---|
| 128 | ! amplitude of distortion |
---|
| 129 | ! j |
---|
| 130 | ! k |
---|
| 131 | ! l |
---|
| 132 | ! m |
---|
| 133 | ! ---------------------------------------------------------------- |
---|
| 134 | ! 75 1 1 multipole order |
---|
| 135 | ! cosine part of Hamiltonian |
---|
| 136 | ! sine part of Hamiltonian |
---|
| 137 | ! amplitude of Hamiltonian |
---|
| 138 | ! j |
---|
| 139 | ! k |
---|
| 140 | ! l |
---|
| 141 | ! m |
---|
| 142 | ! ---------------------------------------------------------------- |
---|
| 143 | ! 76 2 1 multipole order |
---|
| 144 | ! appearance number in position range |
---|
| 145 | ! number of resonance |
---|
| 146 | ! position |
---|
| 147 | ! cosine part of Hamiltonian |
---|
| 148 | ! sine part of Hamiltonian |
---|
| 149 | ! amplitude of Hamiltonian |
---|
| 150 | ! j |
---|
| 151 | ! k |
---|
| 152 | ! l |
---|
| 153 | ! m |
---|
| 154 | ! ---------------------------------------------------------------- |
---|
| 155 | ! 77 2 1 multipole order |
---|
| 156 | ! appearance number in position range |
---|
| 157 | ! number of resonance |
---|
| 158 | ! position |
---|
| 159 | ! cosine part of distortion |
---|
| 160 | ! sine part of distortion |
---|
| 161 | ! amplitude of distortion |
---|
| 162 | ! j |
---|
| 163 | ! k |
---|
| 164 | ! l |
---|
| 165 | ! m |
---|
| 166 | ! ---------------------------------------------------------------- |
---|
| 167 | ! |
---|
| 168 | ! Program Distort2 |
---|
| 169 | ! |
---|
| 170 | ! Multipole |
---|
| 171 | ! Strength |
---|
| 172 | ! Unit iu_on Order Contents in each column |
---|
| 173 | ! ---------------------------------------------------------------- |
---|
| 174 | ! 78 1 2 first multipole order |
---|
| 175 | ! second multipole order |
---|
| 176 | ! cosine part of distortion |
---|
| 177 | ! sine part of distortion |
---|
| 178 | ! amplitude of distortion |
---|
| 179 | ! j |
---|
| 180 | ! k |
---|
| 181 | ! l |
---|
| 182 | ! m |
---|
| 183 | ! ---------------------------------------------------------------- |
---|
| 184 | ! 79 1 2 first multipole order |
---|
| 185 | ! second multipole order |
---|
| 186 | ! cosine part of Hamiltonian |
---|
| 187 | ! sine part of Hamiltonian |
---|
| 188 | ! amplitude of Hamiltonian |
---|
| 189 | ! j |
---|
| 190 | ! k |
---|
| 191 | ! l |
---|
| 192 | ! m |
---|
| 193 | ! ---------------------------------------------------------------- |
---|
| 194 | ! |
---|
| 195 | !--------------------------------------------------------------------- |
---|
| 196 | |
---|
| 197 | implicit none |
---|
| 198 | |
---|
| 199 | integer ierr |
---|
| 200 | |
---|
| 201 | integer ier,iprog,mh,mmul,mmul2,mmult,mmultf,mmultw,mmultx, & |
---|
| 202 | &n3,nblz |
---|
| 203 | double precision etl20,pieni |
---|
| 204 | integer icc,icd,ifacd2,ifact4,ifacta,ihv,iplane,isig,itij,ityc, & |
---|
| 205 | &iu_on,n1,n2 |
---|
| 206 | double precision beta,betexp,bstr,cosav,det,det1,etl,etl1,etl2, & |
---|
| 207 | &fac2,fac4,facd2,fact,fact0,factb,four,ham,one, & |
---|
| 208 | &phi,pi,pi2,pihi,qx,qy,sbeta,sign,two,zero |
---|
| 209 | real tlim,time0,time1,time2,time3,time4 |
---|
| 210 | real tim1,tim2,tim3,tim4,tim5,tim6,tim7,tim8 |
---|
| 211 | integer j70,j71,j72,j73,j74,j75,j76,j77,j78,j79 |
---|
| 212 | integer table_size_70,table_size_71,table_size_72,table_size_73, & |
---|
| 213 | &table_size_74,table_size_75,table_size_76,table_size_77, & |
---|
| 214 | &table_size_78,table_size_79 |
---|
| 215 | character*16 strn |
---|
| 216 | character*18 comment |
---|
| 217 | parameter(pieni=1d-34,nblz=4000,mmul=11,mh=100) |
---|
| 218 | parameter(mmul2=mmul+1,mmult=8*mmul**2) |
---|
| 219 | parameter(mmultw=2*mmul,mmultx=mmult/7,mmultf=2*mmultx*mmult) |
---|
| 220 | |
---|
| 221 | !--- szcompar is the size of the arrays |
---|
| 222 | !--- returned by the routine comm_para |
---|
| 223 | integer szcompar |
---|
| 224 | parameter (szcompar = 100) |
---|
| 225 | |
---|
| 226 | !--- szchara is the size of the character strings char_a |
---|
| 227 | !--- returned by the routine comm_para |
---|
| 228 | integer szchara |
---|
| 229 | parameter (szchara = 400) |
---|
| 230 | |
---|
| 231 | common/c0/ pi,pi2,pihi,zero,one,two,four,comment(-mmul:mmul) |
---|
| 232 | common/c1/ tlim,time0,time1,time2,time3,time4 |
---|
| 233 | common/c2/ tim1,tim2,tim3,tim4,tim5,tim6,tim7,tim8 |
---|
| 234 | common/c3/ bstr(-mmul:mmul,nblz),strn(-mmul:mmul,nblz) |
---|
| 235 | common/c4/ beta(2,-mmul:mmul,nblz),sbeta(2,-mmul:mmul,nblz) |
---|
| 236 | common/c4q/ phi(2,-mmul:mmul,nblz) |
---|
| 237 | common/c5/ etl(-mmul:mmul,nblz),ityc(-mmul:mmul) |
---|
| 238 | common/c6/ etl1,etl2,qx,qy,iu_on,n1,n2 |
---|
| 239 | common/c7/ icc,icd(mmultx),itij(mmultx,mmult,6),isig(mmultx,mmult) |
---|
| 240 | common/c8/ betexp(mmultx,4),fac4(mmultx,mmult),fac2(mmultx) |
---|
| 241 | common/c9/ fact(0:mmul,0:mmul2),sign(4,4),cosav(0:mmul) |
---|
| 242 | common/c10/ ifacd2(4,mmultf),facd2(2,mmultf),ham(2,mmultf) |
---|
| 243 | common/c11/ det1(2,0:mmul),det(2,-mmul:mmul,0:mmul,0:mmul) |
---|
| 244 | common/c12/ ihv(mmultx),iplane(mmultx,2) |
---|
| 245 | common/c13/ fact0(-mmul:mmul,mh),factb(2,-mmul:mmul,mh) |
---|
| 246 | common/c14/ ifacta(2,-mmul:mmul,mh),ifact4(4,-mmul:mmul,mh) |
---|
| 247 | common/indeces/j70,j71,j72,j73,j74,j75,j76,j77,j78,j79 |
---|
| 248 | common/table_sizes/table_size_70,table_size_71,table_size_72, & |
---|
| 249 | &table_size_73,table_size_74,table_size_75, & |
---|
| 250 | &table_size_76,table_size_77,table_size_78, & |
---|
| 251 | &table_size_79 |
---|
| 252 | |
---|
| 253 | integer detun,disto1,disto2,prend,prblup |
---|
| 254 | integer nint, ndble, k, int_arr(szcompar), char_l(szcompar) |
---|
| 255 | double precision d_arr(szcompar) |
---|
| 256 | character * (szchara) char_a |
---|
| 257 | |
---|
| 258 | !--------------------------------------------------------------------- |
---|
| 259 | |
---|
| 260 | ierr = 0 |
---|
| 261 | j70 = 0 |
---|
| 262 | j71 = 0 |
---|
| 263 | j72 = 0 |
---|
| 264 | j73 = 0 |
---|
| 265 | j74 = 0 |
---|
| 266 | j75 = 0 |
---|
| 267 | j76 = 0 |
---|
| 268 | j77 = 0 |
---|
| 269 | j78 = 0 |
---|
| 270 | j79 = 0 |
---|
| 271 | table_size_70 = 2 |
---|
| 272 | table_size_71 = 2 |
---|
| 273 | table_size_72 = 2 |
---|
| 274 | table_size_73 = 2 |
---|
| 275 | table_size_74 = 2 |
---|
| 276 | table_size_75 = 2 |
---|
| 277 | table_size_76 = 2 |
---|
| 278 | table_size_77 = 2 |
---|
| 279 | table_size_78 = 2 |
---|
| 280 | table_size_79 = 2 |
---|
| 281 | |
---|
| 282 | !--- get detun |
---|
| 283 | |
---|
| 284 | detun = 0 |
---|
| 285 | call comm_para('detune ', nint, ndble, k, int_arr, d_arr, & |
---|
| 286 | &char_a, char_l) |
---|
| 287 | if (nint .gt. 0) detun = int_arr(1) |
---|
| 288 | |
---|
| 289 | !--- get disto1 |
---|
| 290 | |
---|
| 291 | disto1 = 0 |
---|
| 292 | call comm_para('distort1 ', nint, ndble, k, int_arr, d_arr, & |
---|
| 293 | &char_a, char_l) |
---|
| 294 | if (nint .gt. 0) disto1 = int_arr(1) |
---|
| 295 | |
---|
| 296 | !--- get disto2 |
---|
| 297 | |
---|
| 298 | disto2 = 0 |
---|
| 299 | call comm_para('distort2 ', nint, ndble, k, int_arr, d_arr, & |
---|
| 300 | &char_a, char_l) |
---|
| 301 | if (nint .gt. 0) disto2 = int_arr(1) |
---|
| 302 | |
---|
| 303 | iprog = detun + disto1*2 + disto2*4 |
---|
| 304 | if(iprog.le.0.or.iprog.gt.7) call prror(0,14,iprog) |
---|
| 305 | if(iprog.eq.1.or.iprog.eq.3.or.iprog.eq.5.or.iprog.eq.7) & |
---|
| 306 | &write(6,*) ' Program <DETUNE> will be executed ' |
---|
| 307 | if(iprog.eq.2.or.iprog.eq.3.or.iprog.eq.6.or.iprog.eq.7) & |
---|
| 308 | &write(6,*) ' Program <DISTORT1> will be executed ' |
---|
| 309 | if(iprog.eq.4.or.iprog.eq.5.or.iprog.eq.6.or.iprog.eq.7) & |
---|
| 310 | &write(6,*) ' Program <DISTORT2> will be executed ' |
---|
| 311 | |
---|
| 312 | |
---|
| 313 | char_a = ' ' |
---|
| 314 | call comm_para('multipole_order_range ', nint, ndble, k, int_arr, & |
---|
| 315 | &d_arr,char_a, char_l) |
---|
| 316 | n1 = int_arr(1) |
---|
| 317 | n2 = int_arr(2) |
---|
| 318 | if(n1.gt.n2) then |
---|
| 319 | n3=n1 |
---|
| 320 | n1=n2 |
---|
| 321 | n2=n3 |
---|
| 322 | endif |
---|
| 323 | if(abs(n1).gt.mmul.or.abs(n2).gt.mmul) call prror(0,15,mmul) |
---|
| 324 | |
---|
| 325 | !--- get etl1 & etl2 (start position & stop position) |
---|
| 326 | |
---|
| 327 | etl1 = 0. |
---|
| 328 | etl2 = 0. |
---|
| 329 | call comm_para('start_stop ', nint, ndble, k, int_arr, d_arr, & |
---|
| 330 | &char_a, char_l) |
---|
| 331 | if (ndble .gt. 0) then |
---|
| 332 | etl1 = d_arr(1) |
---|
| 333 | etl2 = d_arr(2) |
---|
| 334 | endif |
---|
| 335 | |
---|
| 336 | print *,"Start and End Position in [m]: ", etl1,etl2 |
---|
| 337 | |
---|
| 338 | if(etl1.lt.0) etl1=-etl1 |
---|
| 339 | if(etl2.lt.0) etl2=-etl2 |
---|
| 340 | if(etl2.lt.etl1) then |
---|
| 341 | etl20=etl2 |
---|
| 342 | etl1=etl2 |
---|
| 343 | etl2=etl20 |
---|
| 344 | endif |
---|
| 345 | |
---|
| 346 | !--- get no print |
---|
| 347 | |
---|
| 348 | iu_on = 0 |
---|
| 349 | call comm_para('noprint ', nint, ndble, k, int_arr, d_arr, & |
---|
| 350 | &char_a, char_l) |
---|
| 351 | |
---|
| 352 | if (int_arr(1) .eq. 0) then |
---|
| 353 | |
---|
| 354 | !--- get prend |
---|
| 355 | |
---|
| 356 | prend = 0 |
---|
| 357 | call comm_para('print_at_end ', nint, ndble, k, int_arr, d_arr, & |
---|
| 358 | &char_a, char_l) |
---|
| 359 | if (nint .gt. 0) prend = int_arr(1) |
---|
| 360 | |
---|
| 361 | !--- get prblup |
---|
| 362 | |
---|
| 363 | prblup = 0 |
---|
| 364 | call comm_para('print_all ', nint, ndble, k, int_arr, d_arr, & |
---|
| 365 | &char_a, char_l) |
---|
| 366 | if (nint .gt. 0) prblup = int_arr(1) |
---|
| 367 | |
---|
| 368 | iu_on = prend + prblup*2 |
---|
| 369 | endif |
---|
| 370 | |
---|
| 371 | if(iu_on.lt.0.or.iu_on.gt.3) iu_on=0 |
---|
| 372 | if(iu_on.eq.0) write(6,*) ' No print_out requested ' |
---|
| 373 | if(iu_on.eq.1.or.iu_on.eq.3) then |
---|
| 374 | write(6,*) ' Print_out at end of structure ' |
---|
| 375 | if(detun.eq.1) then |
---|
| 376 | open(70,file='detune_1_end',form='formatted',status='unknown',& |
---|
| 377 | &iostat=ier) |
---|
| 378 | if(ier.ne.0) call prror(0,16,70) |
---|
| 379 | write(70,10070) |
---|
| 380 | open(72,file='detune_2_hor',form='formatted',status='unknown',& |
---|
| 381 | &iostat=ier) |
---|
| 382 | if(ier.ne.0) call prror(0,16,72) |
---|
| 383 | write(72,10072) |
---|
| 384 | open(73,file='detune_2_ver',form='formatted', & |
---|
| 385 | &status='unknown',iostat=ier) |
---|
| 386 | if(ier.ne.0) call prror(0,16,73) |
---|
| 387 | write(73,10073) |
---|
| 388 | endif |
---|
| 389 | if(disto1.eq.1) then |
---|
| 390 | open(74,file='distort_1_f_end',form='formatted', & |
---|
| 391 | &status='unknown',iostat=ier) |
---|
| 392 | if(ier.ne.0) call prror(0,16,74) |
---|
| 393 | write(74,10074) |
---|
| 394 | open(75,file='distort_1_h_end',form='formatted', & |
---|
| 395 | &status='unknown',iostat=ier) |
---|
| 396 | if(ier.ne.0) call prror(0,16,75) |
---|
| 397 | write(75,10074) |
---|
| 398 | endif |
---|
| 399 | if(disto2.eq.1) then |
---|
| 400 | open(78,file='distort_2_f_end',form='formatted', & |
---|
| 401 | &status='unknown',iostat=ier) |
---|
| 402 | if(ier.ne.0) call prror(0,16,78) |
---|
| 403 | write(78,10078) |
---|
| 404 | open(79,file='distort_2_h_end',form='formatted', & |
---|
| 405 | &status='unknown',iostat=ier) |
---|
| 406 | if(ier.ne.0) call prror(0,16,79) |
---|
| 407 | write(79,10078) |
---|
| 408 | endif |
---|
| 409 | endif |
---|
| 410 | if(iu_on.eq.2.or.iu_on.eq.3) then |
---|
| 411 | write(6,*) ' Print_out at each multipole ' |
---|
| 412 | if(detun.eq.1) then |
---|
| 413 | open(71,file='detune_1_all',form='formatted', & |
---|
| 414 | &status='unknown',iostat=ier) |
---|
| 415 | if(ier.ne.0) call prror(0,16,71) |
---|
| 416 | write(71,10070) |
---|
| 417 | endif |
---|
| 418 | if(disto1.eq.1) then |
---|
| 419 | open(76,file='distort_1_f_all',form='formatted', & |
---|
| 420 | &status='unknown',iostat=ier) |
---|
| 421 | if(ier.ne.0) call prror(0,16,76) |
---|
| 422 | write(76,10076) |
---|
| 423 | open(77,file='distort_1_h_all',form='formatted', & |
---|
| 424 | &status='unknown',iostat=ier) |
---|
| 425 | if(ier.ne.0) call prror(0,16,77) |
---|
| 426 | write(77,10076) |
---|
| 427 | endif |
---|
| 428 | endif |
---|
| 429 | open(34,file='fc.34',form='formatted',status='unknown', & |
---|
| 430 | &recl=8192,iostat=ier) |
---|
| 431 | if(ier.ne.0) call prror(0,16,34) |
---|
| 432 | !----------------------------------------------------------------------- |
---|
| 433 | ! Initialisation |
---|
| 434 | !----------------------------------------------------------------------- |
---|
| 435 | call init |
---|
| 436 | call timest(tlim) |
---|
| 437 | call timex(tim1) |
---|
| 438 | !----------------------------------------------------------------------- |
---|
| 439 | ! Read Data |
---|
| 440 | !----------------------------------------------------------------------- |
---|
| 441 | call readdat |
---|
| 442 | !----------------------------------------------------------------------- |
---|
| 443 | |
---|
| 444 | call timex(tim2) |
---|
| 445 | write(6,10000) tim2-tim1 |
---|
| 446 | if(iprog.eq.1.or.iprog.eq.3.or.iprog.eq.5.or.iprog.eq.7) & |
---|
| 447 | &call detune |
---|
| 448 | call timex(tim2) |
---|
| 449 | if(iprog.eq.2.or.iprog.eq.3.or.iprog.eq.6.or.iprog.eq.7) & |
---|
| 450 | &call distort1 |
---|
| 451 | call timex(tim2) |
---|
| 452 | if(iprog.eq.4.or.iprog.eq.5.or.iprog.eq.6.or.iprog.eq.7) & |
---|
| 453 | &call distort2 |
---|
| 454 | call timex(tim6) |
---|
| 455 | write(6,10010) tim6-tim1 |
---|
| 456 | |
---|
| 457 | continue |
---|
| 458 | if(detun.eq.1) then |
---|
| 459 | close(70,status='keep') |
---|
| 460 | close(71,status='keep') |
---|
| 461 | close(72,status='keep') |
---|
| 462 | close(73,status='keep') |
---|
| 463 | endif |
---|
| 464 | if(disto1.eq.1) then |
---|
| 465 | close(74,status='keep') |
---|
| 466 | close(75,status='keep') |
---|
| 467 | close(76,status='keep') |
---|
| 468 | close(77,status='keep') |
---|
| 469 | endif |
---|
| 470 | if(disto2.eq.1) then |
---|
| 471 | close(78,status='keep') |
---|
| 472 | close(79,status='keep') |
---|
| 473 | endif |
---|
| 474 | |
---|
| 475 | 10000 format(/80('-')//'Reading Data took ',f10.3,' second(s)', & |
---|
| 476 | &' of Execution Time'//80('-')//) |
---|
| 477 | 10010 format(/80('-')//'Program <SODD> used a total of ',f10.3, & |
---|
| 478 | &' second(s) of Execution Time'//80('-')//) |
---|
| 479 | 10070 format('MulOrd ',' Plane',4x,'Hor/Vert detuning',3x,'Hinv Vinv') |
---|
| 480 | 10072 format('MulOrd1 ','MulOrd2 ',3x,'Horizontal detuning', & |
---|
| 481 | &2x,'Hinv Vinv') |
---|
| 482 | 10073 format('MulOrd1 ','MulOrd2 ',3x,'Vertical detuning ', & |
---|
| 483 | &2x,'Hinv Vinv') |
---|
| 484 | 10074 format('MulOrd',11x,'Cosine',17x,'Sine',19x,'Amplitude', & |
---|
| 485 | &8x,'j',3x,'k',3x,'l',3x,'m',1x) |
---|
| 486 | 10076 format('MulOrd',3x,'Loc',2x,'Res',6x,'Position[m]',14x,'Cosine', & |
---|
| 487 | &17x,'Sine',19x,'Amplitude', & |
---|
| 488 | &8x,'j',3x,'k',3x,'l',3x,'m',1x) |
---|
| 489 | 10078 format('MulOrd1',' MulOrd2',10x,'Cosine',17x,'Sine',19x, & |
---|
| 490 | &'Amplitude',8x,'j',3x,'k',3x,'l',3x,'m',1x) |
---|
| 491 | return |
---|
| 492 | end subroutine soddin |
---|
| 493 | |
---|
| 494 | subroutine detune |
---|
| 495 | !----------------------------------------------------------------------- |
---|
| 496 | !--Program Detune |
---|
| 497 | !----------------------------------------------------------------------- |
---|
| 498 | implicit none |
---|
| 499 | integer i,ia,iaa,iah,iamin,ib,ic,id,ih1,ii,imu,imuty,j,jj,jstep,k,& |
---|
| 500 | &k1,k2,k3,l,l0,l1,l2,l3,m,mh,mmul,mmul2,mmult,mmultf,mmultw, & |
---|
| 501 | &mmultx,nblz |
---|
| 502 | double precision betdx,betdxi,fac,facd0,facd1,phix,phiy,pieni, & |
---|
| 503 | &sigde,sinar,tij |
---|
| 504 | integer icc,icd,ifacd2,ifact4,ifacta,ihv,iplane,isig,itij,ityc, & |
---|
| 505 | &iu_on,n1,n2 |
---|
| 506 | double precision beta,betexp,bstr,cosav,det,det1,etl,etl1,etl2, & |
---|
| 507 | &fac2,fac4,facd2,fact,fact0,factb,four,ham,one,phi,pi,pi2,pihi,qx, & |
---|
| 508 | &qy,sbeta,sign,two,zero |
---|
| 509 | real tlim,time0,time1,time2,time3,time4 |
---|
| 510 | real tim1,tim2,tim3,tim4,tim5,tim6,tim7,tim8 |
---|
| 511 | character*16 strn |
---|
| 512 | character*18 comment |
---|
| 513 | parameter(pieni=1d-34,nblz=4000,mmul=11,mh=100) |
---|
| 514 | parameter(mmul2=mmul+1,mmult=8*mmul**2) |
---|
| 515 | parameter(mmultw=2*mmul,mmultx=mmult/7,mmultf=2*mmultx*mmult) |
---|
| 516 | common/c0/ pi,pi2,pihi,zero,one,two,four,comment(-mmul:mmul) |
---|
| 517 | common/c1/ tlim,time0,time1,time2,time3,time4 |
---|
| 518 | common/c2/ tim1,tim2,tim3,tim4,tim5,tim6,tim7,tim8 |
---|
| 519 | common/c3/ bstr(-mmul:mmul,nblz),strn(-mmul:mmul,nblz) |
---|
| 520 | common/c4/ beta(2,-mmul:mmul,nblz),sbeta(2,-mmul:mmul,nblz) |
---|
| 521 | common/c4q/ phi(2,-mmul:mmul,nblz) |
---|
| 522 | common/c5/ etl(-mmul:mmul,nblz),ityc(-mmul:mmul) |
---|
| 523 | common/c6/ etl1,etl2,qx,qy,iu_on,n1,n2 |
---|
| 524 | common/c7/ icc,icd(mmultx),itij(mmultx,mmult,6),isig(mmultx,mmult) |
---|
| 525 | common/c8/ betexp(mmultx,4),fac4(mmultx,mmult),fac2(mmultx) |
---|
| 526 | common/c9/ fact(0:mmul,0:mmul2),sign(4,4),cosav(0:mmul) |
---|
| 527 | common/c10/ ifacd2(4,mmultf),facd2(2,mmultf),ham(2,mmultf) |
---|
| 528 | common/c11/ det1(2,0:mmul),det(2,-mmul:mmul,0:mmul,0:mmul) |
---|
| 529 | common/c12/ ihv(mmultx),iplane(mmultx,2) |
---|
| 530 | common/c13/ fact0(-mmul:mmul,mh),factb(2,-mmul:mmul,mh) |
---|
| 531 | common/c14/ ifacta(2,-mmul:mmul,mh),ifact4(4,-mmul:mmul,mh) |
---|
| 532 | dimension tij(0:mmul2,-mmul2:mmul2) |
---|
| 533 | dimension betdxi(mmultx),betdx(mmultx) |
---|
| 534 | dimension imu(2) |
---|
| 535 | !--------------------------------------------------------------------- |
---|
| 536 | iamin=2 |
---|
| 537 | jstep=2 |
---|
| 538 | !----------------------------------------------------------------------- |
---|
| 539 | ! Read Data |
---|
| 540 | !----------------------------------------------------------------------- |
---|
| 541 | do i=0,mmul2 |
---|
| 542 | do j=-mmul2,mmul2 |
---|
| 543 | tij(i,j)=zero |
---|
| 544 | enddo |
---|
| 545 | enddo |
---|
| 546 | !--------------------------------------------------------------------- |
---|
| 547 | !-- |
---|
| 548 | !--Detuning First Order in Multipole Strength |
---|
| 549 | !-- |
---|
| 550 | !--------------------------------------------------------------------- |
---|
| 551 | do ia=2,mmul,jstep |
---|
| 552 | imu(1)=ia |
---|
| 553 | imu(2)=0 |
---|
| 554 | iah=ia/2 |
---|
| 555 | if(ityc(ia).gt.0) then |
---|
| 556 | write(6,10000) 'Structure has: ',ityc(ia),' single ', & |
---|
| 557 | &comment(ia) |
---|
| 558 | ih1=ia/2-1 |
---|
| 559 | do i=1,ityc(ia) |
---|
| 560 | fac=-two/dble(ia)*bstr(ia,i) |
---|
| 561 | sigde=-one |
---|
| 562 | do j=1,iah |
---|
| 563 | jj=j-1 |
---|
| 564 | sigde=-sigde |
---|
| 565 | det1(1,j)=sigde*fac*cosav(jj)*cosav(iah-jj)*dble(iah-jj)* & |
---|
| 566 | &fact(ia,1+jj*2)*beta(1,ia,i)**(iah-jj)*beta(2,ia,i)**jj |
---|
| 567 | det(1,ia,iah-j,jj)=det(1,ia,iah-j,jj)+det1(1,j) |
---|
| 568 | if(j.eq.iah) then |
---|
| 569 | sigde=-sigde |
---|
| 570 | det1(2,iah)=sigde*fac*dble(iah)*cosav(iah)* & |
---|
| 571 | &beta(2,ia,i)**iah |
---|
| 572 | det(2,ia,0,iah-1)=det(2,ia,0,iah-1)+det1(2,iah) |
---|
| 573 | endif |
---|
| 574 | enddo |
---|
| 575 | if(iu_on.gt.0) call detwri(0,imu,ih1) |
---|
| 576 | enddo |
---|
| 577 | !--------------------------------------------------------------------- |
---|
| 578 | !-- |
---|
| 579 | !--Printing First Order Terms |
---|
| 580 | !-- |
---|
| 581 | !--------------------------------------------------------------------- |
---|
| 582 | call detwri(1,imu,ih1) |
---|
| 583 | endif |
---|
| 584 | enddo |
---|
| 585 | call timex(tim3) |
---|
| 586 | write(6,10030) tim3-tim2 |
---|
| 587 | !--------------------------------------------------------------------- |
---|
| 588 | !-- |
---|
| 589 | !--Calculation of Detuning Poisson Brackets |
---|
| 590 | !-- |
---|
| 591 | !-----Needed for Second Order Terms |
---|
| 592 | !-- |
---|
| 593 | !--------------------------------------------------------------------- |
---|
| 594 | do imuty=1,2 |
---|
| 595 | call timex(tim3) |
---|
| 596 | !--------------------------------------------------------------------- |
---|
| 597 | !-- |
---|
| 598 | !--Detuning Second Order in Multipole Strength |
---|
| 599 | !-- |
---|
| 600 | !--------------------------------------------------------------------- |
---|
| 601 | do ia=iamin,mmul |
---|
| 602 | do ib=ia,mmul,jstep |
---|
| 603 | if(imuty.eq.1) ic=ia |
---|
| 604 | if(imuty.eq.1) id=ib |
---|
| 605 | if(imuty.eq.2) ic=-ia |
---|
| 606 | if(imuty.eq.2) id=-ib |
---|
| 607 | if(ityc(ic).gt.0.and.ityc(id).gt.0) then |
---|
| 608 | imu(1)=ic |
---|
| 609 | imu(2)=id |
---|
| 610 | do k=1,mmultx |
---|
| 611 | ihv(k)=0 |
---|
| 612 | iplane(k,1)=0 |
---|
| 613 | iplane(k,2)=0 |
---|
| 614 | betexp(k,1)=zero |
---|
| 615 | betexp(k,2)=zero |
---|
| 616 | betexp(k,3)=zero |
---|
| 617 | betexp(k,4)=zero |
---|
| 618 | fac2(k)=zero |
---|
| 619 | do m=1,mmult |
---|
| 620 | itij(k,m,1)=0 |
---|
| 621 | itij(k,m,2)=0 |
---|
| 622 | itij(k,m,3)=0 |
---|
| 623 | fac4(k,m)=zero |
---|
| 624 | enddo |
---|
| 625 | enddo |
---|
| 626 | call caldet2(ia,ib,imuty) |
---|
| 627 | do ii=1,ityc(ic) |
---|
| 628 | do k=1,icc |
---|
| 629 | betdxi(k)=(beta(1,ic,ii)**betexp(k,1))* & |
---|
| 630 | &(beta(2,ic,ii)**betexp(k,2)) |
---|
| 631 | enddo |
---|
| 632 | do jj=1,ityc(id) |
---|
| 633 | if(ii.eq.1.and.jj.eq.1) then |
---|
| 634 | if(ia.ne.ib) then |
---|
| 635 | write(6,10050) 'Pairs of: ',ityc(ic),comment(ic), & |
---|
| 636 | &' and: ',ityc(id),comment(id) |
---|
| 637 | else |
---|
| 638 | write(6,10060) 'Quadratic Effect of: ', & |
---|
| 639 | &ityc(ic),comment(ic) |
---|
| 640 | endif |
---|
| 641 | endif |
---|
| 642 | !--Phase Differences |
---|
| 643 | phix=abs(phi(1,id,jj)-phi(1,ic,ii))-qx |
---|
| 644 | phiy=abs(phi(2,id,jj)-phi(2,ic,ii))-qy |
---|
| 645 | !--Phase Factor |
---|
| 646 | if(imuty.eq.1) iaa=0 |
---|
| 647 | if(imuty.eq.2) iaa=1 |
---|
| 648 | do i=ia+iaa,0,-2 |
---|
| 649 | do j=i-ia,ia-i,2 |
---|
| 650 | if(i.ne.0.or.j.gt.0) then |
---|
| 651 | sinar=sin(dble(i)*qx+dble(j)*qy) |
---|
| 652 | if(abs(sinar).gt.pieni) then |
---|
| 653 | tij(i,j)=cos(dble(i)*phix+dble(j)*phiy)/sinar |
---|
| 654 | else |
---|
| 655 | write(6,*) ' Warning: Results for detuning ',& |
---|
| 656 | &'probably wrong; sinar= ',sinar |
---|
| 657 | tij(i,j)=zero |
---|
| 658 | endif |
---|
| 659 | endif |
---|
| 660 | enddo |
---|
| 661 | enddo |
---|
| 662 | !--Detuning Calculation |
---|
| 663 | facd0=bstr(ic,ii)*bstr(id,jj) |
---|
| 664 | do k=1,icc |
---|
| 665 | k1=iplane(k,1) |
---|
| 666 | k2=iplane(k,2) |
---|
| 667 | k3=ihv(k) |
---|
| 668 | if(k3.eq.1.or.(k3.eq.2.and.k1.eq.0)) then |
---|
| 669 | betdx(k)=betdxi(k)* & |
---|
| 670 | &(beta(1,id,jj)**betexp(k,3))* & |
---|
| 671 | &(beta(2,id,jj)**betexp(k,4)) |
---|
| 672 | facd1=zero |
---|
| 673 | l0=icd(k) |
---|
| 674 | do l=1,l0 |
---|
| 675 | l1=itij(k,l,1) |
---|
| 676 | l2=itij(k,l,2) |
---|
| 677 | l3=itij(k,l,3) |
---|
| 678 | facd1=facd1+fac4(k,l)*(tij(l2,l3)+ & |
---|
| 679 | &dble(l1)*tij(l2,-l3)) |
---|
| 680 | enddo |
---|
| 681 | det(k3,ic,k1,k2)=det(k3,ic,k1,k2)+facd0* & |
---|
| 682 | &fac2(k)*facd1*betdx(k) |
---|
| 683 | endif |
---|
| 684 | enddo |
---|
| 685 | do i=ia+iaa,0,-2 |
---|
| 686 | do j=i-ia,ia-i,2 |
---|
| 687 | if(i.ne.0.or.j.gt.0) then |
---|
| 688 | tij(i,j)=zero |
---|
| 689 | endif |
---|
| 690 | enddo |
---|
| 691 | enddo |
---|
| 692 | enddo |
---|
| 693 | enddo |
---|
| 694 | !--------------------------------------------------------------------- |
---|
| 695 | !-- |
---|
| 696 | !--Printing Second Order Terms |
---|
| 697 | !-- |
---|
| 698 | !--------------------------------------------------------------------- |
---|
| 699 | ih1=(ia+ib)/2-2 |
---|
| 700 | call detwri(2,imu,ih1) |
---|
| 701 | endif |
---|
| 702 | enddo |
---|
| 703 | enddo |
---|
| 704 | call timex(tim4) |
---|
| 705 | write(6,10090) tim4-tim3 |
---|
| 706 | enddo |
---|
| 707 | call timex(tim5) |
---|
| 708 | write(6,10100) tim5-tim2 |
---|
| 709 | !--------------------------------------------------------------------- |
---|
| 710 | 10000 format(//80('-')//5x,a,i6,2a) |
---|
| 711 | 10030 format(/'First Order Detuning Calculation took ', & |
---|
| 712 | &f10.3,' second(s)',' of Execution Time'//80('-')//) |
---|
| 713 | 10050 format(//80('-')//10x,a,i6,2a,i6,a) |
---|
| 714 | 10060 format(//80('-')//10x,a,i6,a) |
---|
| 715 | 10090 format(/80('-')/'Second Order Detuning Calculation took ', & |
---|
| 716 | &f10.3,' second(s)',' of Execution Time'//80('-')//) |
---|
| 717 | 10100 format(/80('-')//'Program <DETUNE> used ',f10.3, & |
---|
| 718 | &' second(s) of Execution Time'//80('-')//) |
---|
| 719 | return |
---|
| 720 | end subroutine detune |
---|
| 721 | subroutine distort1 |
---|
| 722 | !----------------------------------------------------------------------- |
---|
| 723 | !--Program Distort1 (first order) |
---|
| 724 | !----------------------------------------------------------------------- |
---|
| 725 | implicit none |
---|
| 726 | integer i,ia,iadd,ib,ic,id,id1,id11,id110,id12,id120,idiff1, & |
---|
| 727 | &idiff2,ie,if,ii,iia,iia2,ivar,j,mh,mmul,mmul2,mmult,mmultf,mmultw,& |
---|
| 728 | &mmultx,nblz |
---|
| 729 | double precision distc,dists,dstc,dsts,facc,facs,fact1,factb1, & |
---|
| 730 | &factb2,fpre,hamc,hams,hmc,hms,pieni,signii,sinar,tc,ts |
---|
| 731 | integer icc,icd,ifacd2,ifact4,ifacta,ihv,iplane,isig,itij,ityc, & |
---|
| 732 | &iu_on,n1,n2 |
---|
| 733 | integer j70,j71,j72,j73,j74,j75,j76,j77,j78,j79 |
---|
| 734 | integer table_size_70,table_size_71,table_size_72,table_size_73, & |
---|
| 735 | &table_size_74,table_size_75,table_size_76,table_size_77, & |
---|
| 736 | &table_size_78,table_size_79 |
---|
| 737 | integer data_size |
---|
| 738 | double precision beta,betexp,bstr,cosav,det,det1,etl,etl1,etl2, & |
---|
| 739 | &fac2,fac4,facd2,fact,fact0,factb,four,ham,one,phi,pi,pi2,pihi,qx, & |
---|
| 740 | &qy,sbeta,sign,two,zero,amp |
---|
| 741 | real tlim,time0,time1,time2,time3,time4 |
---|
| 742 | real tim1,tim2,tim3,tim4,tim5,tim6,tim7,tim8 |
---|
| 743 | character*16 strn,table_name |
---|
| 744 | character*18 comment |
---|
| 745 | parameter(pieni=1d-34,nblz=4000,mmul=11,mh=100) |
---|
| 746 | parameter(mmul2=mmul+1,mmult=8*mmul**2) |
---|
| 747 | parameter(mmultw=2*mmul,mmultx=mmult/7,mmultf=2*mmultx*mmult) |
---|
| 748 | common/c0/ pi,pi2,pihi,zero,one,two,four,comment(-mmul:mmul) |
---|
| 749 | common/c1/ tlim,time0,time1,time2,time3,time4 |
---|
| 750 | common/c2/ tim1,tim2,tim3,tim4,tim5,tim6,tim7,tim8 |
---|
| 751 | common/c3/ bstr(-mmul:mmul,nblz),strn(-mmul:mmul,nblz) |
---|
| 752 | common/c4/ beta(2,-mmul:mmul,nblz),sbeta(2,-mmul:mmul,nblz) |
---|
| 753 | common/c4q/ phi(2,-mmul:mmul,nblz) |
---|
| 754 | common/c5/ etl(-mmul:mmul,nblz),ityc(-mmul:mmul) |
---|
| 755 | common/c6/ etl1,etl2,qx,qy,iu_on,n1,n2 |
---|
| 756 | common/c7/ icc,icd(mmultx),itij(mmultx,mmult,6),isig(mmultx,mmult) |
---|
| 757 | common/c8/ betexp(mmultx,4),fac4(mmultx,mmult),fac2(mmultx) |
---|
| 758 | common/c9/ fact(0:mmul,0:mmul2),sign(4,4),cosav(0:mmul) |
---|
| 759 | common/c10/ ifacd2(4,mmultf),facd2(2,mmultf),ham(2,mmultf) |
---|
| 760 | common/c11/ det1(2,0:mmul),det(2,-mmul:mmul,0:mmul,0:mmul) |
---|
| 761 | common/c12/ ihv(mmultx),iplane(mmultx,2) |
---|
| 762 | common/c13/ fact0(-mmul:mmul,mh),factb(2,-mmul:mmul,mh) |
---|
| 763 | common/c14/ ifacta(2,-mmul:mmul,mh),ifact4(4,-mmul:mmul,mh) |
---|
| 764 | common/indeces/j70,j71,j72,j73,j74,j75,j76,j77,j78,j79 |
---|
| 765 | common/table_sizes/table_size_70,table_size_71,table_size_72, & |
---|
| 766 | &table_size_73,table_size_74,table_size_75, & |
---|
| 767 | &table_size_76,table_size_77,table_size_78, & |
---|
| 768 | &table_size_79 |
---|
| 769 | integer int_to_write(11) |
---|
| 770 | double precision double_to_write(11) |
---|
| 771 | |
---|
| 772 | dimension distc(-mmul:mmul,mh),dists(-mmul:mmul,mh) |
---|
| 773 | dimension dstc(mh),dsts(mh),ivar(-mmul:mmul),fact1(mh) |
---|
| 774 | dimension hamc(-mmul:mmul,mh),hams(-mmul:mmul,mh) |
---|
| 775 | !--------------------------------------------------------------------- |
---|
| 776 | do i=1,mh |
---|
| 777 | fact1(i)=zero |
---|
| 778 | do j=-mmul,mmul |
---|
| 779 | ivar(j)=0 |
---|
| 780 | ifacta(1,j,i)=0 |
---|
| 781 | ifacta(2,j,i)=0 |
---|
| 782 | ifact4(1,j,i)=0 |
---|
| 783 | ifact4(2,j,i)=0 |
---|
| 784 | ifact4(3,j,i)=0 |
---|
| 785 | ifact4(4,j,i)=0 |
---|
| 786 | factb(1,j,i)=zero |
---|
| 787 | factb(2,j,i)=zero |
---|
| 788 | fact0(j,i)=zero |
---|
| 789 | distc(j,i)=zero |
---|
| 790 | dists(j,i)=zero |
---|
| 791 | hamc(j,i)=zero |
---|
| 792 | hams(j,i)=zero |
---|
| 793 | enddo |
---|
| 794 | enddo |
---|
| 795 | iadd=2 |
---|
| 796 | idiff1=1 |
---|
| 797 | idiff2=2 |
---|
| 798 | ivar(1)=1 |
---|
| 799 | ivar(-1)=1 |
---|
| 800 | !--Number of Resonances |
---|
| 801 | do i=2,mmul |
---|
| 802 | ivar(i)=iadd |
---|
| 803 | ivar(-i)=iadd |
---|
| 804 | if(mod(i,2).eq.0) then |
---|
| 805 | idiff1=idiff1+idiff2 |
---|
| 806 | idiff2=idiff2+1 |
---|
| 807 | endif |
---|
| 808 | iadd=iadd+idiff1 |
---|
| 809 | enddo |
---|
| 810 | !--Calculation of the various Factors for Multipoles -mmul -> +mmul |
---|
| 811 | do ia=-mmul,mmul |
---|
| 812 | if(ityc(ia).gt.0) then |
---|
| 813 | iia=iabs(ia) |
---|
| 814 | fpre=one/dble((iia*2**(iia+1))) |
---|
| 815 | ib=0 |
---|
| 816 | signii=-one |
---|
| 817 | if(ia.lt.0) then |
---|
| 818 | iia2=iia-1 |
---|
| 819 | else |
---|
| 820 | iia2=iia |
---|
| 821 | endif |
---|
| 822 | do ic=iia2,0,-2 |
---|
| 823 | signii=signii*(-one) |
---|
| 824 | ib=ib+1 |
---|
| 825 | id=iia-ic |
---|
| 826 | id1=id+1 |
---|
| 827 | ifacta(1,ia,ib)=ic |
---|
| 828 | ifacta(2,ia,ib)=id |
---|
| 829 | ifact4(1,ia,ib)=ic |
---|
| 830 | ifact4(2,ia,ib)=0 |
---|
| 831 | ifact4(3,ia,ib)=id |
---|
| 832 | ifact4(4,ia,ib)=0 |
---|
| 833 | factb1=abs(ifacta(1,ia,ib))/two |
---|
| 834 | factb2=abs(ifacta(2,ia,ib))/two |
---|
| 835 | factb(1,ia,ib)=factb1 |
---|
| 836 | factb(2,ia,ib)=factb2 |
---|
| 837 | fact0(ia,ib)=signii*fpre*fact(iia,id1) |
---|
| 838 | if(id.gt.0.and.ic.gt.0) then |
---|
| 839 | ib=ib+1 |
---|
| 840 | ifacta(1,ia,ib)=ic |
---|
| 841 | ifacta(2,ia,ib)=-id |
---|
| 842 | ifact4(1,ia,ib)=ic |
---|
| 843 | ifact4(2,ia,ib)=0 |
---|
| 844 | ifact4(3,ia,ib)=0 |
---|
| 845 | ifact4(4,ia,ib)=id |
---|
| 846 | factb(1,ia,ib)=factb1 |
---|
| 847 | factb(2,ia,ib)=factb2 |
---|
| 848 | fact0(ia,ib)=signii*fpre*fact(iia,id1) |
---|
| 849 | endif |
---|
| 850 | id11=1 |
---|
| 851 | id110=0 |
---|
| 852 | do ie=ic,0,-2 |
---|
| 853 | if(id.ne.0.or.ie.ne.0) then |
---|
| 854 | if(ie.ne.ic) then |
---|
| 855 | id110=id110+1 |
---|
| 856 | ib=ib+1 |
---|
| 857 | id11=id11+1 |
---|
| 858 | ifacta(1,ia,ib)=ie |
---|
| 859 | ifacta(2,ia,ib)=id |
---|
| 860 | ifact4(1,ia,ib)=ie+id110 |
---|
| 861 | ifact4(2,ia,ib)=id110 |
---|
| 862 | ifact4(3,ia,ib)=id |
---|
| 863 | ifact4(4,ia,ib)=0 |
---|
| 864 | factb(1,ia,ib)=factb1 |
---|
| 865 | factb(2,ia,ib)=factb2 |
---|
| 866 | fact0(ia,ib)=signii*fpre*fact(iia,id1)* & |
---|
| 867 | &fact(ic,id11) |
---|
| 868 | if(id.gt.0.and.ie.gt.0) then |
---|
| 869 | ib=ib+1 |
---|
| 870 | ifacta(1,ia,ib)=ie |
---|
| 871 | ifacta(2,ia,ib)=-id |
---|
| 872 | ifact4(1,ia,ib)=ie+id110 |
---|
| 873 | ifact4(2,ia,ib)=id110 |
---|
| 874 | ifact4(3,ia,ib)=0 |
---|
| 875 | ifact4(4,ia,ib)=id |
---|
| 876 | factb(1,ia,ib)=factb1 |
---|
| 877 | factb(2,ia,ib)=factb2 |
---|
| 878 | fact0(ia,ib)=signii*fpre*fact(iia,id1)* & |
---|
| 879 | &fact(ic,id11) |
---|
| 880 | endif |
---|
| 881 | endif |
---|
| 882 | id12=1 |
---|
| 883 | id120=0 |
---|
| 884 | do if=id-2,0,-2 |
---|
| 885 | if(if.ne.0.or.ie.ne.0) then |
---|
| 886 | id120=id120+1 |
---|
| 887 | ib=ib+1 |
---|
| 888 | id12=id12+1 |
---|
| 889 | ifacta(1,ia,ib)=ie |
---|
| 890 | ifacta(2,ia,ib)=if |
---|
| 891 | ifact4(1,ia,ib)=ie+id110 |
---|
| 892 | ifact4(2,ia,ib)=id110 |
---|
| 893 | ifact4(3,ia,ib)=if+id120 |
---|
| 894 | ifact4(4,ia,ib)=id120 |
---|
| 895 | factb(1,ia,ib)=factb1 |
---|
| 896 | factb(2,ia,ib)=factb2 |
---|
| 897 | fact0(ia,ib)=signii*fpre*fact(iia,id1)* & |
---|
| 898 | &fact(ic,id11)*fact(id,id12) |
---|
| 899 | if(if.gt.0.and.ie.gt.0) then |
---|
| 900 | ib=ib+1 |
---|
| 901 | ifacta(1,ia,ib)=ie |
---|
| 902 | ifacta(2,ia,ib)=-if |
---|
| 903 | ifact4(1,ia,ib)=ie+id110 |
---|
| 904 | ifact4(2,ia,ib)=id110 |
---|
| 905 | ifact4(3,ia,ib)=id120 |
---|
| 906 | ifact4(4,ia,ib)=if+id120 |
---|
| 907 | factb(1,ia,ib)=factb1 |
---|
| 908 | factb(2,ia,ib)=factb2 |
---|
| 909 | fact0(ia,ib)=signii*fpre*fact(iia,id1)* & |
---|
| 910 | &fact(ic,id11)*fact(id,id12) |
---|
| 911 | endif |
---|
| 912 | endif |
---|
| 913 | enddo |
---|
| 914 | id120=0 |
---|
| 915 | endif |
---|
| 916 | enddo |
---|
| 917 | id110=0 |
---|
| 918 | enddo |
---|
| 919 | call sortres(2,ia,ib) |
---|
| 920 | endif |
---|
| 921 | enddo |
---|
| 922 | !--Calculation of I Order Distortion Function |
---|
| 923 | do iia=2,mmul |
---|
| 924 | do ii=1,2 |
---|
| 925 | if(ii.eq.1) ia=iia |
---|
| 926 | if(ii.eq.2) ia=-iia |
---|
| 927 | if(ityc(ia).gt.0) then |
---|
| 928 | write(6,10100) 'First Order Contribution of: ', & |
---|
| 929 | &ityc(ia),comment(ia) |
---|
| 930 | do i=1,ityc(ia) |
---|
| 931 | do j=1,ivar(ia) |
---|
| 932 | !--Factor of Machine Parameters |
---|
| 933 | fact1(j)=bstr(ia,i)*fact0(ia,j)* & |
---|
| 934 | &beta(1,ia,i)**factb(1,ia,j)*beta(2,ia,i)**factb(2,ia,j) |
---|
| 935 | !--Phase Factor |
---|
| 936 | tc=sin(ifacta(1,ia,j)*(phi(1,ia,i)-qx)+ & |
---|
| 937 | &ifacta(2,ia,j)*(phi(2,ia,i)-qy)) |
---|
| 938 | ts=cos(ifacta(1,ia,j)*(phi(1,ia,i)-qx)+ & |
---|
| 939 | &ifacta(2,ia,j)*(phi(2,ia,i)-qy)) |
---|
| 940 | !--Generating Function Calculation |
---|
| 941 | sinar=sin(ifacta(1,ia,j)*qx+ifacta(2,ia,j)*qy) |
---|
| 942 | if(abs(sinar).gt.pieni) then |
---|
| 943 | dstc(j)=fact1(j)*tc/sinar |
---|
| 944 | dsts(j)=fact1(j)*ts/sinar |
---|
| 945 | else |
---|
| 946 | write(6,*) ' Warning: Results for distortion ', & |
---|
| 947 | &'function probably wrong; sinar= ',sinar |
---|
| 948 | dstc(j)=zero |
---|
| 949 | dsts(j)=zero |
---|
| 950 | endif |
---|
| 951 | distc(ia,j)=distc(ia,j)+dstc(j) |
---|
| 952 | dists(ia,j)=dists(ia,j)+dsts(j) |
---|
| 953 | enddo |
---|
| 954 | !--Write Generating Function Calculation of each Element |
---|
| 955 | if(iu_on.eq.2.or.iu_on.eq.3) then |
---|
| 956 | table_name ='distort_1_f_all ' |
---|
| 957 | data_size = ivar(ia) |
---|
| 958 | call fit_table(table_name,data_size,table_size_76) |
---|
| 959 | do j=1,data_size |
---|
| 960 | amp = dsqrt(dstc(j)*dstc(j)+dsts(j)*dsts(j)) |
---|
| 961 | write(76,10000) ia,i,j,etl(ia,i),dstc(j),dsts(j),amp, & |
---|
| 962 | &ifact4(1,ia,j),ifact4(2,ia,j), & |
---|
| 963 | &ifact4(3,ia,j),ifact4(4,ia,j) |
---|
| 964 | data_size = ivar(ia) |
---|
| 965 | int_to_write(1) = ia |
---|
| 966 | int_to_write(2) = i |
---|
| 967 | int_to_write(3) = j |
---|
| 968 | double_to_write(4) = etl(ia,i) |
---|
| 969 | double_to_write(5) = dstc(j) |
---|
| 970 | double_to_write(6) = dsts(j) |
---|
| 971 | double_to_write(7) = amp |
---|
| 972 | int_to_write(8) = ifact4(1,ia,j) |
---|
| 973 | int_to_write(9) = ifact4(2,ia,j) |
---|
| 974 | int_to_write(10) = ifact4(3,ia,j) |
---|
| 975 | int_to_write(11) = ifact4(4,ia,j) |
---|
| 976 | call write_table(table_name,11,int_to_write, & |
---|
| 977 | &double_to_write) |
---|
| 978 | j76 = j76 + 1 |
---|
| 979 | call augment_count(table_name) |
---|
| 980 | enddo |
---|
| 981 | endif |
---|
| 982 | !--Calculate and write Hamiltonian of each Element |
---|
| 983 | if(iu_on.eq.2.or.iu_on.eq.3) then |
---|
| 984 | table_name ='distort_1_h_all ' |
---|
| 985 | data_size = ivar(ia) |
---|
| 986 | call fit_table(table_name,data_size,table_size_77) |
---|
| 987 | do j=1,data_size |
---|
| 988 | facc=one-cos(two*(ifacta(1,ia,j)*qx+ & |
---|
| 989 | &ifacta(2,ia,j)*qy)) |
---|
| 990 | facs=sin(two*(ifacta(1,ia,j)*qx+ifacta(2,ia,j)*qy)) |
---|
| 991 | hmc=facc*dstc(j)+facs*dsts(j) |
---|
| 992 | hms=facc*dsts(j)-facs*dstc(j) |
---|
| 993 | double_to_write(5)=hmc |
---|
| 994 | double_to_write(6)=hms |
---|
| 995 | amp = dsqrt(hmc*hmc+hms*hms) |
---|
| 996 | write(77,10000) ia,i,j,etl(ia,i),double_to_write(5), & |
---|
| 997 | &double_to_write(6),amp,ifact4(1,ia,j),ifact4(2,ia,j), & |
---|
| 998 | &ifact4(3,ia,j),ifact4(4,ia,j) |
---|
| 999 | int_to_write(1) = ia |
---|
| 1000 | int_to_write(2) = i |
---|
| 1001 | int_to_write(3) = j |
---|
| 1002 | double_to_write(4) = etl(ia,i) |
---|
| 1003 | double_to_write(7) = amp |
---|
| 1004 | int_to_write(8) = ifact4(1,ia,j) |
---|
| 1005 | int_to_write(9) = ifact4(2,ia,j) |
---|
| 1006 | int_to_write(10) = ifact4(3,ia,j) |
---|
| 1007 | int_to_write(11) = ifact4(4,ia,j) |
---|
| 1008 | call write_table(table_name,11,int_to_write, & |
---|
| 1009 | &double_to_write) |
---|
| 1010 | j77 = j77 + 1 |
---|
| 1011 | call augment_count(table_name) |
---|
| 1012 | enddo |
---|
| 1013 | endif |
---|
| 1014 | enddo |
---|
| 1015 | !--Write total Generating Function |
---|
| 1016 | if(iu_on.eq.1.or.iu_on.eq.3) then |
---|
| 1017 | write(6,10020) |
---|
| 1018 | table_name ='distort_1_f_end ' |
---|
| 1019 | data_size = ivar(ia) |
---|
| 1020 | call fit_table(table_name,data_size,table_size_74) |
---|
| 1021 | do j=1,data_size |
---|
| 1022 | amp = sqrt(distc(ia,j)*distc(ia,j)+dists(ia,j)* & |
---|
| 1023 | &dists(ia,j)) |
---|
| 1024 | write(6,10040) distc(ia,j),dists(ia,j),amp, & |
---|
| 1025 | &ifact4(1,ia,j),ifact4(2,ia,j),ifact4(3,ia,j), & |
---|
| 1026 | &ifact4(4,ia,j) |
---|
| 1027 | write(74,10010) ia,distc(ia,j),dists(ia,j),amp, & |
---|
| 1028 | &ifact4(1,ia,j),ifact4(2,ia,j),ifact4(3,ia,j), & |
---|
| 1029 | &ifact4(4,ia,j) |
---|
| 1030 | int_to_write(1) = ia |
---|
| 1031 | double_to_write(2) = distc(ia,j) |
---|
| 1032 | double_to_write(3) = dists(ia,j) |
---|
| 1033 | double_to_write(4) = amp |
---|
| 1034 | int_to_write(5) = ifact4(1,ia,j) |
---|
| 1035 | int_to_write(6) = ifact4(2,ia,j) |
---|
| 1036 | int_to_write(7) = ifact4(3,ia,j) |
---|
| 1037 | int_to_write(8) = ifact4(4,ia,j) |
---|
| 1038 | call write_table(table_name,8,int_to_write, & |
---|
| 1039 | &double_to_write) |
---|
| 1040 | j74 = j74 + 1 |
---|
| 1041 | call augment_count(table_name) |
---|
| 1042 | enddo |
---|
| 1043 | !--Calculate and write total Hamiltonian |
---|
| 1044 | write(6,10030) |
---|
| 1045 | table_name ='distort_1_h_end ' |
---|
| 1046 | data_size = ivar(ia) |
---|
| 1047 | call fit_table(table_name,data_size,table_size_75) |
---|
| 1048 | do j=1,data_size |
---|
| 1049 | facc=one-cos(two*(ifacta(1,ia,j)*qx+ifacta(2,ia,j)*qy)) |
---|
| 1050 | facs=sin(two*(ifacta(1,ia,j)*qx+ifacta(2,ia,j)*qy)) |
---|
| 1051 | hamc(ia,j)=facc*distc(ia,j)+facs*dists(ia,j) |
---|
| 1052 | hams(ia,j)=facc*dists(ia,j)-facs*distc(ia,j) |
---|
| 1053 | double_to_write(2)=hamc(ia,j) |
---|
| 1054 | double_to_write(3)=hams(ia,j) |
---|
| 1055 | amp = sqrt(hamc(ia,j)*hamc(ia,j)+hams(ia,j)*hams(ia,j)) |
---|
| 1056 | write(6,10040) double_to_write(2),double_to_write(3), & |
---|
| 1057 | &,ifact4(1,ia,j),ifact4(2,ia,j), & |
---|
| 1058 | &ifact4(3,ia,j),ifact4(4,ia,j) |
---|
| 1059 | write(75,10010) ia,double_to_write(2), & |
---|
| 1060 | &double_to_write(3),amp,ifact4(1,ia,j),ifact4(2,ia,j), & |
---|
| 1061 | &ifact4(3,ia,j),ifact4(4,ia,j) |
---|
| 1062 | int_to_write(1) = ia |
---|
| 1063 | double_to_write(4) = amp |
---|
| 1064 | int_to_write(5) = ifact4(1,ia,j) |
---|
| 1065 | int_to_write(6) = ifact4(2,ia,j) |
---|
| 1066 | int_to_write(7) = ifact4(3,ia,j) |
---|
| 1067 | int_to_write(8) = ifact4(4,ia,j) |
---|
| 1068 | call write_table(table_name,8,int_to_write, & |
---|
| 1069 | &double_to_write) |
---|
| 1070 | j75 = j75 + 1 |
---|
| 1071 | call augment_count(table_name) |
---|
| 1072 | enddo |
---|
| 1073 | endif |
---|
| 1074 | endif |
---|
| 1075 | enddo |
---|
| 1076 | enddo |
---|
| 1077 | call timex(tim5) |
---|
| 1078 | write(6,10110) tim5-tim2 |
---|
| 1079 | !--------------------------------------------------------------------- |
---|
| 1080 | 10000 format(i4,i7,i5,4(1pe23.15),4i4) |
---|
| 1081 | 10010 format(i4,4x,3(1pe23.15),6i4) |
---|
| 1082 | 10020 format(/80('-')/10x,'Distortionfunction Terms (first order)'/ & |
---|
| 1083 | &80('-')/4x,'Cosine-Term',8x,'Sine-Term',8x,'Amplitude',7x,'J', & |
---|
| 1084 | &5x,'K',5x,'L',5x,'M'/80('-')) |
---|
| 1085 | 10030 format(/80('-')/10x,'Hamiltonian Terms (first order)'/80('-')/ & |
---|
| 1086 | &4x,'Cosine-Term',8x,'Sine-Term',8x,'Amplitude',7x,'J', & |
---|
| 1087 | &5x,'K',5x,'L',5x,'M'/80('-')) |
---|
| 1088 | 10040 format(3(1pe17.9),4i6) |
---|
| 1089 | 10100 format(//80('-')//10x,a,i8,a) |
---|
| 1090 | 10110 format(/80('-')//'Program <DISTORT1> used ',f10.3, & |
---|
| 1091 | &' second(s) of Execution Time'//80('-')//) |
---|
| 1092 | return |
---|
| 1093 | end subroutine distort1 |
---|
| 1094 | subroutine distort2 |
---|
| 1095 | !----------------------------------------------------------------------- |
---|
| 1096 | !--Program Distort2 (second order) |
---|
| 1097 | !----------------------------------------------------------------------- |
---|
| 1098 | implicit none |
---|
| 1099 | integer i,ia,ia0,ib,ib0,ic,ic0,ic00,ica,ica2,id,id0,id00,ide, & |
---|
| 1100 | &ii,imax,imuty,itij0,jj,k,k1,k2,l,l0,l3,l4,ll3,ll4,m,mde,mh,mmul, & |
---|
| 1101 | &mmul2,mmult,mmultf,mmultw,mmultx,nblz,num |
---|
| 1102 | double precision arg,argij,betdx,betdxi,dl12q,dl34q,efact,ephix1, & |
---|
| 1103 | &ephiy1,esin34,etrgij,facc,facd0,facd1,facd10,facp,facs,ffff,phix1,& |
---|
| 1104 | &phix2,phix20,phiy1,phiy2,phiy20,pieni,sin34,spij,trgij |
---|
| 1105 | integer icc,icd,ifacd2,ifact4,ifacta,ihv,iplane,isig,itij,ityc, & |
---|
| 1106 | &iu_on,n1,n2 |
---|
| 1107 | integer j70,j71,j72,j73,j74,j75,j76,j77,j78,j79 |
---|
| 1108 | integer table_size_70,table_size_71,table_size_72,table_size_73, & |
---|
| 1109 | &table_size_74,table_size_75,table_size_76,table_size_77, & |
---|
| 1110 | &table_size_78,table_size_79 |
---|
| 1111 | integer data_size |
---|
| 1112 | double precision beta,betexp,bstr,cosav,det,det1,etl,etl1,etl2, & |
---|
| 1113 | &fac2,fac4,facd2,fact,fact0,factb,four,ham,one,phi,pi,pi2,pihi,qx, & |
---|
| 1114 | &qy,sbeta,sign,two,zero |
---|
| 1115 | real tlim,time0,time1,time2,time3,time4 |
---|
| 1116 | real tim1,tim2,tim3,tim4,tim5,tim6,tim7,tim8 |
---|
| 1117 | character*16 strn,table_name |
---|
| 1118 | character*18 comment |
---|
| 1119 | parameter(pieni=1d-34,nblz=4000,mmul=11,mh=100) |
---|
| 1120 | parameter(mmul2=mmul+1,mmult=8*mmul**2) |
---|
| 1121 | parameter(mmultw=2*mmul,mmultx=mmult/7,mmultf=2*mmultx*mmult) |
---|
| 1122 | common/c0/ pi,pi2,pihi,zero,one,two,four,comment(-mmul:mmul) |
---|
| 1123 | common/c1/ tlim,time0,time1,time2,time3,time4 |
---|
| 1124 | common/c2/ tim1,tim2,tim3,tim4,tim5,tim6,tim7,tim8 |
---|
| 1125 | common/c3/ bstr(-mmul:mmul,nblz),strn(-mmul:mmul,nblz) |
---|
| 1126 | common/c4/ beta(2,-mmul:mmul,nblz),sbeta(2,-mmul:mmul,nblz) |
---|
| 1127 | common/c4q/ phi(2,-mmul:mmul,nblz) |
---|
| 1128 | common/c5/ etl(-mmul:mmul,nblz),ityc(-mmul:mmul) |
---|
| 1129 | common/c6/ etl1,etl2,qx,qy,iu_on,n1,n2 |
---|
| 1130 | common/c7/ icc,icd(mmultx),itij(mmultx,mmult,6),isig(mmultx,mmult) |
---|
| 1131 | common/c8/ betexp(mmultx,4),fac4(mmultx,mmult),fac2(mmultx) |
---|
| 1132 | common/c9/ fact(0:mmul,0:mmul2),sign(4,4),cosav(0:mmul) |
---|
| 1133 | common/c10/ ifacd2(4,mmultf),facd2(2,mmultf),ham(2,mmultf) |
---|
| 1134 | common/c11/ det1(2,0:mmul),det(2,-mmul:mmul,0:mmul,0:mmul) |
---|
| 1135 | common/c12/ ihv(mmultx),iplane(mmultx,2) |
---|
| 1136 | common/c13/ fact0(-mmul:mmul,mh),factb(2,-mmul:mmul,mh) |
---|
| 1137 | common/c14/ ifacta(2,-mmul:mmul,mh),ifact4(4,-mmul:mmul,mh) |
---|
| 1138 | common/indeces/j70,j71,j72,j73,j74,j75,j76,j77,j78,j79 |
---|
| 1139 | common/table_sizes/table_size_70,table_size_71,table_size_72, & |
---|
| 1140 | &table_size_73,table_size_74,table_size_75, & |
---|
| 1141 | &table_size_76,table_size_77,table_size_78, & |
---|
| 1142 | &table_size_79 |
---|
| 1143 | dimension efact(mmultx,mmult),ll3(mmultx,mmult),ll4(mmultx,mmult) |
---|
| 1144 | dimension betdxi(mmultx),betdx(mmultx) |
---|
| 1145 | dimension facd1(mmultx,mmult,2) |
---|
| 1146 | dimension itij0(mmultx,mmult,4),facd10(mmultx,mmult,2) |
---|
| 1147 | dimension mde(mmul) |
---|
| 1148 | integer int_to_write(11) |
---|
| 1149 | double precision double_to_write(11) |
---|
| 1150 | !--------------------------------------------------------------------- |
---|
| 1151 | !--Set the rest of the Variables to zero |
---|
| 1152 | do i=1,mmult |
---|
| 1153 | do ii=1,mmultx |
---|
| 1154 | facd1(ii,i,1)=zero |
---|
| 1155 | facd1(ii,i,2)=zero |
---|
| 1156 | enddo |
---|
| 1157 | enddo |
---|
| 1158 | do i=1,mmul |
---|
| 1159 | mde(i)=0 |
---|
| 1160 | enddo |
---|
| 1161 | ii=0 |
---|
| 1162 | do i=4,mmul,2 |
---|
| 1163 | ii=ii+1 |
---|
| 1164 | mde(ii)=i |
---|
| 1165 | enddo |
---|
| 1166 | !--------------------------------------------------------------------- |
---|
| 1167 | !-- |
---|
| 1168 | !--Calculation of Driving-term Poisson Brackets |
---|
| 1169 | !-- |
---|
| 1170 | !-----Needed for Second Order Terms |
---|
| 1171 | !-- |
---|
| 1172 | !--------------------------------------------------------------------- |
---|
| 1173 | !--------------------------------------------------------------------- |
---|
| 1174 | !-- |
---|
| 1175 | !--Driving & Hamiltonian Terms Second Order in Multipole Strength |
---|
| 1176 | !-- |
---|
| 1177 | !--------------------------------------------------------------------- |
---|
| 1178 | do ic00=n1,n2 |
---|
| 1179 | ic=ic00 |
---|
| 1180 | do id00=-iabs(ic),iabs(ic) |
---|
| 1181 | id=id00 |
---|
| 1182 | ia=iabs(ic) |
---|
| 1183 | ib=iabs(id) |
---|
| 1184 | !--Difference between erect and skew |
---|
| 1185 | !--imuty=1 ==> "erect"; imuty=2 ==> "skew"; |
---|
| 1186 | !--imuty=3 ==> "erect-skew"; imuty=4 ==> "skew-erect" |
---|
| 1187 | if(ic.gt.0.and.id.gt.0) imuty=1 |
---|
| 1188 | if(ic.lt.0.and.id.lt.0) imuty=2 |
---|
| 1189 | if(ic.gt.0.and.id.lt.0) imuty=3 |
---|
| 1190 | if(ic.lt.0.and.id.gt.0) imuty=4 |
---|
| 1191 | call timex(tim3) |
---|
| 1192 | if(ityc(ic).gt.0.and.ityc(id).gt.0) then |
---|
| 1193 | icc=0 |
---|
| 1194 | imax=ia+ib-2 |
---|
| 1195 | do k=1,mmultx |
---|
| 1196 | betexp(k,1)=zero |
---|
| 1197 | betexp(k,2)=zero |
---|
| 1198 | betexp(k,3)=zero |
---|
| 1199 | betexp(k,4)=zero |
---|
| 1200 | icd(k)=0 |
---|
| 1201 | fac2(k)=zero |
---|
| 1202 | do m=1,mmult |
---|
| 1203 | isig(k,m)=0 |
---|
| 1204 | itij(k,m,1)=0 |
---|
| 1205 | itij(k,m,2)=0 |
---|
| 1206 | itij(k,m,3)=0 |
---|
| 1207 | itij(k,m,4)=0 |
---|
| 1208 | itij(k,m,5)=0 |
---|
| 1209 | itij(k,m,6)=0 |
---|
| 1210 | itij0(k,m,1)=0 |
---|
| 1211 | itij0(k,m,2)=0 |
---|
| 1212 | itij0(k,m,3)=0 |
---|
| 1213 | itij0(k,m,4)=0 |
---|
| 1214 | fac4(k,m)=zero |
---|
| 1215 | enddo |
---|
| 1216 | enddo |
---|
| 1217 | !-- |
---|
| 1218 | !--The ide(n0=0,yes=1) switch checks if the average(detuning)terms |
---|
| 1219 | !--create contributions to the distortion function in higher order |
---|
| 1220 | !--A full loop is needed |
---|
| 1221 | !-- |
---|
| 1222 | ide=0 |
---|
| 1223 | if(imuty.eq.1.and.ia.ne.ib) then |
---|
| 1224 | do k=1,mmul |
---|
| 1225 | if(ib.eq.mde(k)) ide=1 |
---|
| 1226 | enddo |
---|
| 1227 | endif |
---|
| 1228 | if(ide.eq.1) then |
---|
| 1229 | !-- |
---|
| 1230 | !--Extra (ide=1) Loop starts here |
---|
| 1231 | !-- |
---|
| 1232 | ia0=ia |
---|
| 1233 | ib0=ib |
---|
| 1234 | ic0=ic |
---|
| 1235 | id0=id |
---|
| 1236 | ia=ib0 |
---|
| 1237 | ib=ia0 |
---|
| 1238 | ic=id0 |
---|
| 1239 | id=ic0 |
---|
| 1240 | call caldt2(ia,ib,imuty) |
---|
| 1241 | do k=1,mmultx |
---|
| 1242 | do m=1,mmult |
---|
| 1243 | itij0(k,m,1)=itij(k,m,3) |
---|
| 1244 | itij0(k,m,2)=itij(k,m,4) |
---|
| 1245 | itij0(k,m,3)=itij(k,m,5) |
---|
| 1246 | itij0(k,m,4)=itij(k,m,6) |
---|
| 1247 | enddo |
---|
| 1248 | enddo |
---|
| 1249 | do ii=1,ityc(ic) |
---|
| 1250 | !--Beta term of first set of multipoles |
---|
| 1251 | do k=1,icc |
---|
| 1252 | betdxi(k)=(beta(1,ic,ii)**betexp(k,1))* & |
---|
| 1253 | &(beta(2,ic,ii)**betexp(k,2)) |
---|
| 1254 | enddo |
---|
| 1255 | do jj=1,ityc(id) |
---|
| 1256 | facd0=two*bstr(ic,ii)*bstr(id,jj) |
---|
| 1257 | !--Phase Differences (phix2,phiy2) of secondary resonances change |
---|
| 1258 | !--sign for reversed order of the position of multipoles |
---|
| 1259 | spij=one |
---|
| 1260 | if(phi(1,id,jj).lt.phi(1,ic,ii)) spij=-one |
---|
| 1261 | !--Drivingterm Calculation |
---|
| 1262 | do k=1,icc |
---|
| 1263 | !--Beta term of second set of multipoles multiplied to first |
---|
| 1264 | betdx(k)=betdxi(k)* & |
---|
| 1265 | &(beta(1,id,jj)**betexp(k,3))* & |
---|
| 1266 | &(beta(2,id,jj)**betexp(k,4)) |
---|
| 1267 | l0=icd(k) |
---|
| 1268 | !--Cosine and sine part of driving term |
---|
| 1269 | do l=1,l0 |
---|
| 1270 | if(isig(k,l).eq.1) then |
---|
| 1271 | dl12q=itij(k,l,1)*qx+itij(k,l,2)*qy |
---|
| 1272 | if(abs(dl12q).lt.pieni) dl12q=pieni |
---|
| 1273 | l3=itij(k,l,3)-itij(k,l,4) |
---|
| 1274 | l4=itij(k,l,5)-itij(k,l,6) |
---|
| 1275 | dl34q=l3*qx+l4*qy |
---|
| 1276 | if(abs(dl34q).lt.pieni) dl34q=pieni |
---|
| 1277 | sin34=sin(dl34q) |
---|
| 1278 | phix1=-qx |
---|
| 1279 | phiy1=-qy |
---|
| 1280 | trgij=sin(dl12q) |
---|
| 1281 | if(phi(1,id,jj).eq.phi(1,ic,ii)) then |
---|
| 1282 | phix2=phi(1,id,jj)+qx |
---|
| 1283 | phiy2=phi(2,id,jj)+qy |
---|
| 1284 | else |
---|
| 1285 | phix2=phi(1,id,jj)-spij*qx |
---|
| 1286 | phiy2=phi(2,id,jj)-spij*qy |
---|
| 1287 | endif |
---|
| 1288 | arg=l3*phix1+l4*phiy1+itij(k,l,1)* & |
---|
| 1289 | &phix2+itij(k,l,2)*phiy2 |
---|
| 1290 | facp=fac4(k,l)*facd0*fac2(k)* & |
---|
| 1291 | &betdx(k)/trgij/sin34 |
---|
| 1292 | facd10(k,l,1)=facd10(k,l,1)+facp*sin(arg) |
---|
| 1293 | facd10(k,l,2)=facd10(k,l,2)+facp*cos(arg) |
---|
| 1294 | endif |
---|
| 1295 | enddo |
---|
| 1296 | enddo |
---|
| 1297 | enddo |
---|
| 1298 | enddo |
---|
| 1299 | ia=ia0 |
---|
| 1300 | ib=ib0 |
---|
| 1301 | ic=ic0 |
---|
| 1302 | id=id0 |
---|
| 1303 | endif |
---|
| 1304 | !--------------------------------------------------------------------- |
---|
| 1305 | !-- |
---|
| 1306 | !--Main loop |
---|
| 1307 | !-- |
---|
| 1308 | !--------------------------------------------------------------------- |
---|
| 1309 | icc=0 |
---|
| 1310 | do k=1,mmultx |
---|
| 1311 | betexp(k,1)=zero |
---|
| 1312 | betexp(k,2)=zero |
---|
| 1313 | betexp(k,3)=zero |
---|
| 1314 | betexp(k,4)=zero |
---|
| 1315 | icd(k)=0 |
---|
| 1316 | fac2(k)=zero |
---|
| 1317 | do m=1,mmult |
---|
| 1318 | isig(k,m)=0 |
---|
| 1319 | itij(k,m,1)=0 |
---|
| 1320 | itij(k,m,2)=0 |
---|
| 1321 | itij(k,m,3)=0 |
---|
| 1322 | itij(k,m,4)=0 |
---|
| 1323 | itij(k,m,5)=0 |
---|
| 1324 | itij(k,m,6)=0 |
---|
| 1325 | fac4(k,m)=zero |
---|
| 1326 | enddo |
---|
| 1327 | enddo |
---|
| 1328 | call caldt2(ia,ib,imuty) |
---|
| 1329 | do k=1,icc |
---|
| 1330 | l0=icd(k) |
---|
| 1331 | do l=1,l0 |
---|
| 1332 | dl12q=itij(k,l,1)*qx+itij(k,l,2)*qy |
---|
| 1333 | if(abs(dl12q).lt.pieni) dl12q=pieni |
---|
| 1334 | etrgij=sin(dl12q) |
---|
| 1335 | ll3(k,l)=itij(k,l,3)-itij(k,l,4) |
---|
| 1336 | ll4(k,l)=itij(k,l,5)-itij(k,l,6) |
---|
| 1337 | dl34q=ll3(k,l)*qx+ll4(k,l)*qy |
---|
| 1338 | if(abs(dl34q).lt.pieni) dl34q=pieni |
---|
| 1339 | esin34=sin(dl34q) |
---|
| 1340 | efact(k,l)=fac4(k,l)/etrgij/esin34 |
---|
| 1341 | enddo |
---|
| 1342 | enddo |
---|
| 1343 | do ii=1,ityc(ic) |
---|
| 1344 | !--Beta term of first set of multipoles |
---|
| 1345 | do k=1,icc |
---|
| 1346 | betdxi(k)=(beta(1,ic,ii)**betexp(k,1))* & |
---|
| 1347 | &(beta(2,ic,ii)**betexp(k,2)) |
---|
| 1348 | enddo |
---|
| 1349 | do jj=1,ityc(id) |
---|
| 1350 | facd0=bstr(ic,ii)*bstr(id,jj) |
---|
| 1351 | if(ic.ne.id.and.ic.ne.-id) facd0=facd0*two |
---|
| 1352 | if(ii.eq.1.and.jj.eq.1) then |
---|
| 1353 | if(ic.ne.id) then |
---|
| 1354 | write(6,10090) 'Pairs of: ',ityc(ic),comment(ic), & |
---|
| 1355 | &' and: ',ityc(id),comment(id) |
---|
| 1356 | else |
---|
| 1357 | write(6,10100) 'Quadratic Effect of: ', & |
---|
| 1358 | &ityc(ic),comment(ic) |
---|
| 1359 | endif |
---|
| 1360 | endif |
---|
| 1361 | !--Phase Differences (phix2,phiy2) of secondary resonances change |
---|
| 1362 | !--sign for reversed order of the position of multipoles |
---|
| 1363 | spij=one |
---|
| 1364 | if(phi(1,id,jj).lt.phi(1,ic,ii)) spij=-one |
---|
| 1365 | phix20=phi(1,id,jj)-phi(1,ic,ii)-spij*qx |
---|
| 1366 | phiy20=phi(2,id,jj)-phi(2,ic,ii)-spij*qy |
---|
| 1367 | ephix1=phi(1,ic,ii)-qx |
---|
| 1368 | ephiy1=phi(2,ic,ii)-qy |
---|
| 1369 | !--Drivingterm Calculation |
---|
| 1370 | do k=1,icc |
---|
| 1371 | !--Beta term of second set of multipoles multiplied to first |
---|
| 1372 | betdx(k)=betdxi(k)* & |
---|
| 1373 | &(beta(1,id,jj)**betexp(k,3))* & |
---|
| 1374 | &(beta(2,id,jj)**betexp(k,4)) |
---|
| 1375 | ffff=facd0*fac2(k)*betdx(k) |
---|
| 1376 | l0=icd(k) |
---|
| 1377 | num=0 |
---|
| 1378 | do l=1,l0 |
---|
| 1379 | num=num+isig(k,l) |
---|
| 1380 | enddo |
---|
| 1381 | if (num.eq.0) then |
---|
| 1382 | !--Cosine and sine part of driving term |
---|
| 1383 | do l=1,l0 |
---|
| 1384 | arg=ll3(k,l)*ephix1+ll4(k,l)*ephiy1+ & |
---|
| 1385 | &itij(k,l,1)*phix20+itij(k,l,2)*phiy20 |
---|
| 1386 | facp=ffff*efact(k,l) |
---|
| 1387 | facd1(k,l,1)=facd1(k,l,1)+facp*sin(arg) |
---|
| 1388 | facd1(k,l,2)=facd1(k,l,2)+facp*cos(arg) |
---|
| 1389 | enddo |
---|
| 1390 | else |
---|
| 1391 | do l=1,l0 |
---|
| 1392 | if(isig(k,l).eq.0) then |
---|
| 1393 | arg=ll3(k,l)*ephix1+ll4(k,l)*ephiy1+ & |
---|
| 1394 | &itij(k,l,1)*phix20+itij(k,l,2)*phiy20 |
---|
| 1395 | facp=ffff*efact(k,l) |
---|
| 1396 | facd1(k,l,1)=facd1(k,l,1)+facp*sin(arg) |
---|
| 1397 | facd1(k,l,2)=facd1(k,l,2)+facp*cos(arg) |
---|
| 1398 | else |
---|
| 1399 | dl12q=itij(k,l,1)*qx+itij(k,l,2)*qy |
---|
| 1400 | if(abs(dl12q).lt.pieni) dl12q=pieni |
---|
| 1401 | l3=itij(k,l,3)-itij(k,l,4) |
---|
| 1402 | l4=itij(k,l,5)-itij(k,l,6) |
---|
| 1403 | dl34q=l3*qx+l4*qy |
---|
| 1404 | if(abs(dl34q).lt.pieni) dl34q=pieni |
---|
| 1405 | sin34=sin(dl34q) |
---|
| 1406 | phix1=-qx |
---|
| 1407 | phiy1=-qy |
---|
| 1408 | trgij=sin(dl12q) |
---|
| 1409 | phix2=phi(1,id,jj) |
---|
| 1410 | phiy2=phi(2,id,jj) |
---|
| 1411 | if(ic.eq.id.and.phi(1,id,jj).eq.phi(1,ic,ii)) & |
---|
| 1412 | &then |
---|
| 1413 | trgij=trgij/cos(dl12q) |
---|
| 1414 | else |
---|
| 1415 | phix2=phix2-spij*qx |
---|
| 1416 | phiy2=phiy2-spij*qy |
---|
| 1417 | endif |
---|
| 1418 | arg=l3*phix1+l4*phiy1+itij(k,l,1)* & |
---|
| 1419 | &phix2+itij(k,l,2)*phiy2 |
---|
| 1420 | facp=fac4(k,l)*facd0*fac2(k)* & |
---|
| 1421 | &betdx(k)/trgij/sin34 |
---|
| 1422 | !--Factor 2 comes from trigonometric functions |
---|
| 1423 | if(ic.eq.id) facp=two*facp |
---|
| 1424 | facd1(k,l,1)=facd1(k,l,1)+facp*sin(arg) |
---|
| 1425 | facd1(k,l,2)=facd1(k,l,2)+facp*cos(arg) |
---|
| 1426 | endif |
---|
| 1427 | enddo |
---|
| 1428 | endif |
---|
| 1429 | enddo |
---|
| 1430 | enddo |
---|
| 1431 | enddo |
---|
| 1432 | !--Collect a total of ica non-zero terms |
---|
| 1433 | ica=0 |
---|
| 1434 | do k1=1,mmultx |
---|
| 1435 | do k2=1,mmult |
---|
| 1436 | if(abs(facd1(k1,k2,1)).gt.pieni.or. & |
---|
| 1437 | &abs(facd1(k1,k2,2)).gt.pieni) then |
---|
| 1438 | ica=ica+1 |
---|
| 1439 | if(ica.gt.mmultf) call prror(3,10,mmultf) |
---|
| 1440 | facd2(1,ica)=facd1(k1,k2,1) |
---|
| 1441 | facd2(2,ica)=facd1(k1,k2,2) |
---|
| 1442 | facd1(k1,k2,1)=zero |
---|
| 1443 | facd1(k1,k2,2)=zero |
---|
| 1444 | do ii=1,4 |
---|
| 1445 | ifacd2(ii,ica)=itij(k1,k2,ii+2) |
---|
| 1446 | enddo |
---|
| 1447 | endif |
---|
| 1448 | enddo |
---|
| 1449 | enddo |
---|
| 1450 | !--Collect also the additional detuning induced terms |
---|
| 1451 | if(ide.eq.1) then |
---|
| 1452 | do k1=1,mmultx |
---|
| 1453 | do k2=1,mmult |
---|
| 1454 | if(abs(facd10(k1,k2,1)).gt.pieni.or. & |
---|
| 1455 | &abs(facd10(k1,k2,2)).gt.pieni) then |
---|
| 1456 | ica=ica+1 |
---|
| 1457 | if(ica.gt.mmultf) call prror(3,10,mmultf) |
---|
| 1458 | facd2(1,ica)=facd10(k1,k2,1) |
---|
| 1459 | facd2(2,ica)=facd10(k1,k2,2) |
---|
| 1460 | facd10(k1,k2,1)=zero |
---|
| 1461 | facd10(k1,k2,2)=zero |
---|
| 1462 | do ii=1,4 |
---|
| 1463 | ifacd2(ii,ica)=itij0(k1,k2,ii) |
---|
| 1464 | enddo |
---|
| 1465 | endif |
---|
| 1466 | enddo |
---|
| 1467 | enddo |
---|
| 1468 | endif |
---|
| 1469 | !--Add up all contributions related to a specific resonance |
---|
| 1470 | do k1=1,ica |
---|
| 1471 | do k2=k1+1,ica |
---|
| 1472 | if(ifacd2(1,k2).eq.ifacd2(1,k1).and. & |
---|
| 1473 | &ifacd2(2,k2).eq.ifacd2(2,k1).and. & |
---|
| 1474 | &ifacd2(3,k2).eq.ifacd2(3,k1).and. & |
---|
| 1475 | &ifacd2(4,k2).eq.ifacd2(4,k1)) then |
---|
| 1476 | facd2(1,k1)=facd2(1,k1)+facd2(1,k2) |
---|
| 1477 | facd2(2,k1)=facd2(2,k1)+facd2(2,k2) |
---|
| 1478 | facd2(1,k2)=zero |
---|
| 1479 | facd2(2,k2)=zero |
---|
| 1480 | do ii=1,4 |
---|
| 1481 | ifacd2(ii,k2)=0 |
---|
| 1482 | enddo |
---|
| 1483 | endif |
---|
| 1484 | enddo |
---|
| 1485 | enddo |
---|
| 1486 | ica2=0 |
---|
| 1487 | !--Collect the final total of ica2 non-zero resonance terms |
---|
| 1488 | do k1=1,ica |
---|
| 1489 | if(abs(facd2(1,k1)).gt.pieni.or. & |
---|
| 1490 | &abs(facd2(2,k1)).gt.pieni) then |
---|
| 1491 | ica2=ica2+1 |
---|
| 1492 | if(ica2.ne.k1) then |
---|
| 1493 | facd2(1,ica2)=facd2(1,k1) |
---|
| 1494 | facd2(2,ica2)=facd2(2,k1) |
---|
| 1495 | facd2(1,k1)=zero |
---|
| 1496 | facd2(2,k1)=zero |
---|
| 1497 | do ii=1,4 |
---|
| 1498 | ifacd2(ii,ica2)=ifacd2(ii,k1) |
---|
| 1499 | ifacd2(ii,k1)=0 |
---|
| 1500 | enddo |
---|
| 1501 | endif |
---|
| 1502 | endif |
---|
| 1503 | enddo |
---|
| 1504 | !--Order the resonance terms |
---|
| 1505 | call sortres(3,imax,ica2) |
---|
| 1506 | write(6,10060) |
---|
| 1507 | table_name ='distort_2_f_end ' |
---|
| 1508 | data_size = ica2 |
---|
| 1509 | call fit_table(table_name,data_size,table_size_78) |
---|
| 1510 | do k1=1,data_size |
---|
| 1511 | argij=two*((ifacd2(1,k1)-ifacd2(2,k1))*qx+ & |
---|
| 1512 | &(ifacd2(3,k1)-ifacd2(4,k1))*qy) |
---|
| 1513 | facc=one-cos(argij) |
---|
| 1514 | facs=sin(argij) |
---|
| 1515 | ham(1,k1)=facc*facd2(1,k1)+facs*facd2(2,k1) |
---|
| 1516 | ham(2,k1)=facc*facd2(2,k1)-facs*facd2(1,k1) |
---|
| 1517 | double_to_write(5)=dsqrt(facd2(1,k1)*facd2(1,k1)+ & |
---|
| 1518 | &facd2(2,k1)*facd2(2,k1)) |
---|
| 1519 | write(6,10080) facd2(1,k1),facd2(2,k1),double_to_write(5),& |
---|
| 1520 | &ifacd2(1,k1),ifacd2(2,k1),ifacd2(3,k1),ifacd2(4,k1) |
---|
| 1521 | if(iu_on.eq.1.or.iu_on.eq.3) & |
---|
| 1522 | &write(78,10110) ic,id,facd2(1,k1),facd2(2,k1), & |
---|
| 1523 | &double_to_write(5),ifacd2(1,k1),ifacd2(2,k1), & |
---|
| 1524 | &ifacd2(3,k1),ifacd2(4,k1) |
---|
| 1525 | int_to_write(1) = ic |
---|
| 1526 | int_to_write(2) = id |
---|
| 1527 | double_to_write(3) = facd2(1,k1) |
---|
| 1528 | double_to_write(4) = facd2(2,k1) |
---|
| 1529 | int_to_write(6) = ifacd2(1,k1) |
---|
| 1530 | int_to_write(7) = ifacd2(2,k1) |
---|
| 1531 | int_to_write(8) = ifacd2(3,k1) |
---|
| 1532 | int_to_write(9) = ifacd2(4,k1) |
---|
| 1533 | call write_table(table_name,9,int_to_write, & |
---|
| 1534 | &double_to_write) |
---|
| 1535 | j78 = j78 + 1 |
---|
| 1536 | call augment_count(table_name) |
---|
| 1537 | facd2(1,k1)=zero |
---|
| 1538 | facd2(2,k1)=zero |
---|
| 1539 | enddo |
---|
| 1540 | write(6,10070) |
---|
| 1541 | table_name ='distort_2_h_end ' |
---|
| 1542 | data_size = ica2 |
---|
| 1543 | call fit_table(table_name,data_size,table_size_79) |
---|
| 1544 | do k1=1,data_size |
---|
| 1545 | double_to_write(5)=dsqrt(ham(1,k1)*ham(1,k1)+ & |
---|
| 1546 | &ham(2,k1)*ham(2,k1)) |
---|
| 1547 | write(6,10080) ham(1,k1),ham(2,k1),double_to_write(5), & |
---|
| 1548 | &ifacd2(1,k1),ifacd2(2,k1),ifacd2(3,k1),ifacd2(4,k1) |
---|
| 1549 | if(iu_on.eq.1.or.iu_on.eq.3) & |
---|
| 1550 | &write(79,10110) ic,id,ham(1,k1),ham(2,k1), & |
---|
| 1551 | &double_to_write(5),ifacd2(1,k1),ifacd2(2,k1), & |
---|
| 1552 | &ifacd2(3,k1),ifacd2(4,k1) |
---|
| 1553 | int_to_write(1) = ic |
---|
| 1554 | int_to_write(2) = id |
---|
| 1555 | double_to_write(3) = ham(1,k1) |
---|
| 1556 | double_to_write(4) = ham(2,k1) |
---|
| 1557 | int_to_write(6) = ifacd2(1,k1) |
---|
| 1558 | int_to_write(7) = ifacd2(2,k1) |
---|
| 1559 | int_to_write(8) = ifacd2(3,k1) |
---|
| 1560 | int_to_write(9) = ifacd2(4,k1) |
---|
| 1561 | call write_table(table_name,9,int_to_write, & |
---|
| 1562 | &double_to_write) |
---|
| 1563 | j79 = j79 + 1 |
---|
| 1564 | call augment_count(table_name) |
---|
| 1565 | do ii=1,4 |
---|
| 1566 | ifacd2(ii,k1)=0 |
---|
| 1567 | ham(1,k1)=zero |
---|
| 1568 | ham(2,k1)=zero |
---|
| 1569 | enddo |
---|
| 1570 | enddo |
---|
| 1571 | call timex(tim4) |
---|
| 1572 | write(6,10030) tim4-tim3 |
---|
| 1573 | write(6,10020) tim4-tim2 |
---|
| 1574 | endif |
---|
| 1575 | enddo |
---|
| 1576 | enddo |
---|
| 1577 | !--------------------------------------------------------------------- |
---|
| 1578 | call timex(tim5) |
---|
| 1579 | write(6,10040) tim5-tim2 |
---|
| 1580 | return |
---|
| 1581 | !--------------------------------------------------------------------- |
---|
| 1582 | 10020 format(/80('-')//'Intermediate Time in <DISTORT2>: ', & |
---|
| 1583 | &f10.3,' second(s)',' of Execution Time' & |
---|
| 1584 | &//80('-')//) |
---|
| 1585 | 10030 format(/80('-')//'Second Order Drivingterm Calculation took ', & |
---|
| 1586 | &f10.3,' second(s)',' of Execution Time'//80('-')//) |
---|
| 1587 | 10040 format(/80('-')//'Program <DISTORT2> used ',f10.3, & |
---|
| 1588 | &' second(s) of Execution Time'//80('-')//) |
---|
| 1589 | 10060 format(/80('-')/10x,'Distortionfunction Terms (second order)'/ & |
---|
| 1590 | &80('-')/4x,'Cosine-Term',8x,'Sine-Term',8x,'Amplitude',7x, & |
---|
| 1591 | &'J',5x,'K',5x,'L',5x,'M'/80('-')) |
---|
| 1592 | 10070 format(/80('-')/10x,'Hamiltonian Terms (second order)'/80('-')/ & |
---|
| 1593 | &4x,'Cosine-Term',8x,'Sine-Term',8x,'Amplitude',7x, & |
---|
| 1594 | &'J',5x,'K',5x,'L',5x,'M'/80('-')) |
---|
| 1595 | 10080 format(3(1pe17.9),4i6) |
---|
| 1596 | 10090 format(//80('-')//10x,a,i8,2a,i8,a) |
---|
| 1597 | 10100 format(//80('-')//10x,a,i8,a) |
---|
| 1598 | 10110 format(i4,4x,i4,4x,3(1pe23.15),4i4) |
---|
| 1599 | end subroutine distort2 |
---|
| 1600 | subroutine caldet2(ia,ib,imuty) |
---|
| 1601 | !--------------------------------------------------------------------- |
---|
| 1602 | !-- |
---|
| 1603 | !--Calculation of the poisson bracket to derive the detuning function |
---|
| 1604 | !--Second order in the strength of the multipoles |
---|
| 1605 | !-- |
---|
| 1606 | !--------------------------------------------------------------------- |
---|
| 1607 | implicit none |
---|
| 1608 | integer ia,ib,ic,id,ido,if2,ihoff,ihw,ii,ii1,ii11,iiend,iii,iii1, & |
---|
| 1609 | &ijnu,ijnu1,imuty,ipl,iti,ivoff,jj,jj1,jj11,jjend,jjj,jjj1,k,k1,kk,& |
---|
| 1610 | &l,l1,m,mend,mexp,mh,mij,mm,mmul,mmul2,mmult,mmultf,mmultw,mmultx, & |
---|
| 1611 | &n,nblz,nend,nexp,nij,nn |
---|
| 1612 | double precision betl,betxp,fac,fac1,fac3,fcc2,fcc4,pieni,signii, & |
---|
| 1613 | &signjj |
---|
| 1614 | integer icc,icd,ifacd2,ifact4,ifacta,ihv,iplane,isig,itij,ityc, & |
---|
| 1615 | &iu_on,n1,n2 |
---|
| 1616 | double precision beta,betexp,bstr,cosav,det,det1,etl,etl1,etl2, & |
---|
| 1617 | &fac2,fac4,facd2,fact,fact0,factb,four,ham,one,phi,pi,pi2,pihi,qx, & |
---|
| 1618 | &qy,sbeta,sign,two,zero |
---|
| 1619 | real tlim,time0,time1,time2,time3,time4 |
---|
| 1620 | real tim1,tim2,tim3,tim4,tim5,tim6,tim7,tim8 |
---|
| 1621 | character*16 strn |
---|
| 1622 | character*18 comment |
---|
| 1623 | parameter(pieni=1d-34,nblz=4000,mmul=11,mh=100) |
---|
| 1624 | parameter(mmul2=mmul+1,mmult=8*mmul**2) |
---|
| 1625 | parameter(mmultw=2*mmul,mmultx=mmult/7,mmultf=2*mmultx*mmult) |
---|
| 1626 | common/c0/ pi,pi2,pihi,zero,one,two,four,comment(-mmul:mmul) |
---|
| 1627 | common/c1/ tlim,time0,time1,time2,time3,time4 |
---|
| 1628 | common/c2/ tim1,tim2,tim3,tim4,tim5,tim6,tim7,tim8 |
---|
| 1629 | common/c3/ bstr(-mmul:mmul,nblz),strn(-mmul:mmul,nblz) |
---|
| 1630 | common/c4/ beta(2,-mmul:mmul,nblz),sbeta(2,-mmul:mmul,nblz) |
---|
| 1631 | common/c4q/ phi(2,-mmul:mmul,nblz) |
---|
| 1632 | common/c5/ etl(-mmul:mmul,nblz),ityc(-mmul:mmul) |
---|
| 1633 | common/c6/ etl1,etl2,qx,qy,iu_on,n1,n2 |
---|
| 1634 | common/c7/ icc,icd(mmultx),itij(mmultx,mmult,6),isig(mmultx,mmult) |
---|
| 1635 | common/c8/ betexp(mmultx,4),fac4(mmultx,mmult),fac2(mmultx) |
---|
| 1636 | common/c9/ fact(0:mmul,0:mmul2),sign(4,4),cosav(0:mmul) |
---|
| 1637 | common/c10/ ifacd2(4,mmultf),facd2(2,mmultf),ham(2,mmultf) |
---|
| 1638 | common/c11/ det1(2,0:mmul),det(2,-mmul:mmul,0:mmul,0:mmul) |
---|
| 1639 | common/c12/ ihv(mmultx),iplane(mmultx,2) |
---|
| 1640 | common/c13/ fact0(-mmul:mmul,mh),factb(2,-mmul:mmul,mh) |
---|
| 1641 | common/c14/ ifacta(2,-mmul:mmul,mh),ifact4(4,-mmul:mmul,mh) |
---|
| 1642 | !--------------------------------------------------------------------- |
---|
| 1643 | dimension betl(mmultx,4),betxp(2,mmultx,4) |
---|
| 1644 | dimension ihw(2,mmult),ipl(2,mmult,2),iti(2,mmultx,mmult,6) |
---|
| 1645 | dimension fcc2(2,mmultx),fcc4(2,mmultx,mmult) |
---|
| 1646 | do k=1,mmultx |
---|
| 1647 | fcc2(1,k)=zero |
---|
| 1648 | fcc2(2,k)=zero |
---|
| 1649 | do m=1,mmult |
---|
| 1650 | do nn=1,2 |
---|
| 1651 | ihw(nn,k)=0 |
---|
| 1652 | ipl(nn,k,1)=0 |
---|
| 1653 | ipl(nn,k,2)=0 |
---|
| 1654 | iti(nn,k,m,1)=0 |
---|
| 1655 | iti(nn,k,m,2)=0 |
---|
| 1656 | iti(nn,k,m,3)=0 |
---|
| 1657 | betl(k,1)=zero |
---|
| 1658 | betl(k,2)=zero |
---|
| 1659 | betl(k,3)=zero |
---|
| 1660 | betl(k,4)=zero |
---|
| 1661 | betxp(nn,k,1)=zero |
---|
| 1662 | betxp(nn,k,2)=zero |
---|
| 1663 | betxp(nn,k,3)=zero |
---|
| 1664 | betxp(nn,k,4)=zero |
---|
| 1665 | fcc4(nn,k,m)=zero |
---|
| 1666 | enddo |
---|
| 1667 | enddo |
---|
| 1668 | enddo |
---|
| 1669 | fac=-two/dble(4*2**((ia+ib)/2)*ia*ib) |
---|
| 1670 | if(ia.ne.ib) fac=fac*two |
---|
| 1671 | iiend=(ia+2)/2 |
---|
| 1672 | jjend=(ib+2)/2 |
---|
| 1673 | signii=-one |
---|
| 1674 | if(imuty.eq.2.and.2*iiend.eq.ia+2) iiend=iiend-1 |
---|
| 1675 | do ii=1,iiend |
---|
| 1676 | ii1=(ii-1)*2 |
---|
| 1677 | if(imuty.eq.2) ii1=ii1+1 |
---|
| 1678 | signii=signii*(-one) |
---|
| 1679 | signjj=-one |
---|
| 1680 | if(imuty.eq.2.and.2*jjend.eq.ib+2) jjend=jjend-1 |
---|
| 1681 | do jj=1,jjend |
---|
| 1682 | if2=0 |
---|
| 1683 | jj1=(jj-1)*2 |
---|
| 1684 | if(imuty.eq.2) jj1=jj1+1 |
---|
| 1685 | signjj=signjj*(-one) |
---|
| 1686 | fac1=fac*signii*signjj*fact(ia,1+ii1)*fact(ib,1+jj1) |
---|
| 1687 | iii=ia-ii1 |
---|
| 1688 | jjj=ib-jj1 |
---|
| 1689 | ijnu=(ii-1)*jjend+jj |
---|
| 1690 | betl(ijnu,1)=dble(iii)/two |
---|
| 1691 | betl(ijnu,2)=dble(ii1)/two |
---|
| 1692 | betl(ijnu,3)=dble(jjj)/two |
---|
| 1693 | betl(ijnu,4)=dble(jj1)/two |
---|
| 1694 | ihoff=0 |
---|
| 1695 | ivoff=0 |
---|
| 1696 | if(iii.eq.0.and.jjj.ne.0) then |
---|
| 1697 | fac1=fac1*fact(jjj,jjj/2+1) |
---|
| 1698 | ihoff=1 |
---|
| 1699 | endif |
---|
| 1700 | if(iii.ne.0.and.jjj.eq.0) then |
---|
| 1701 | fac1=fac1*fact(iii,iii/2+1) |
---|
| 1702 | ihoff=1 |
---|
| 1703 | endif |
---|
| 1704 | if(ii1.eq.0.and.jj1.ne.0) then |
---|
| 1705 | fac1=fac1*fact(jj1,jj1/2+1) |
---|
| 1706 | ivoff=1 |
---|
| 1707 | endif |
---|
| 1708 | if(ii1.ne.0.and.jj1.eq.0) then |
---|
| 1709 | fac1=fac1*fact(ii1,ii1/2+1) |
---|
| 1710 | ivoff=1 |
---|
| 1711 | endif |
---|
| 1712 | k=(iii+jjj)/2 |
---|
| 1713 | k1=k-1 |
---|
| 1714 | l=(ii1+jj1)/2 |
---|
| 1715 | l1=l-1 |
---|
| 1716 | mij=min(iii,jjj) |
---|
| 1717 | nij=min(ii1,jj1) |
---|
| 1718 | mend=(mij+2)/2 |
---|
| 1719 | nend=(nij+2)/2 |
---|
| 1720 | iii1=0 |
---|
| 1721 | jjj1=0 |
---|
| 1722 | ii11=0 |
---|
| 1723 | jj11=0 |
---|
| 1724 | if(iii.gt.mij) iii1=(iii+2)/2-mend |
---|
| 1725 | if(jjj.gt.mij) jjj1=(jjj+2)/2-mend |
---|
| 1726 | if(ii1.gt.nij) ii11=(ii1+2)/2-nend |
---|
| 1727 | if(jj1.gt.nij) jj11=(jj1+2)/2-nend |
---|
| 1728 | !--------------------------------------------------------------------- |
---|
| 1729 | !-- |
---|
| 1730 | !--Deriving the Derivatives |
---|
| 1731 | !-- |
---|
| 1732 | !--------------------------------------------------------------------- |
---|
| 1733 | if(k.gt.0) then |
---|
| 1734 | if(l.gt.0.and.ivoff.eq.0.and.ihoff.eq.0) then |
---|
| 1735 | !-- |
---|
| 1736 | !--Qx Term first Derivative by Iy for the Detuning by Ix |
---|
| 1737 | !-- |
---|
| 1738 | ido=1 |
---|
| 1739 | !--Finding Double Appearance |
---|
| 1740 | do ijnu1=1,ijnu-1 |
---|
| 1741 | if( & |
---|
| 1742 | &betl(ijnu1,1).eq.betl(ijnu,3).and. & |
---|
| 1743 | &betl(ijnu1,2).eq.betl(ijnu,4).and. & |
---|
| 1744 | &betl(ijnu1,3).eq.betl(ijnu,1).and. & |
---|
| 1745 | &betl(ijnu1,4).eq.betl(ijnu,2)) then |
---|
| 1746 | ido=0 |
---|
| 1747 | if2=if2+1 |
---|
| 1748 | if(if2.gt.2) call prror(1,13,if2) |
---|
| 1749 | fcc2(if2,ijnu1)=two*fcc2(if2,ijnu1) |
---|
| 1750 | endif |
---|
| 1751 | enddo |
---|
| 1752 | if(ido.eq.1) then |
---|
| 1753 | if2=if2+1 |
---|
| 1754 | if(if2.gt.2) call prror(1,13,if2) |
---|
| 1755 | ihw(if2,ijnu)=1 |
---|
| 1756 | ipl(if2,ijnu,1)=k1 |
---|
| 1757 | ipl(if2,ijnu,2)=l1 |
---|
| 1758 | betxp(if2,ijnu,1)=betl(ijnu,1) |
---|
| 1759 | betxp(if2,ijnu,2)=betl(ijnu,2) |
---|
| 1760 | betxp(if2,ijnu,3)=betl(ijnu,3) |
---|
| 1761 | betxp(if2,ijnu,4)=betl(ijnu,4) |
---|
| 1762 | fcc2(if2,ijnu)=fac1*dble(k*ii1*jj1)/dble(2**(k1+l1)) |
---|
| 1763 | do m=1,mend |
---|
| 1764 | mexp=mij-(m-1)*2 |
---|
| 1765 | fac3=fact(iii,iii1+m)*fact(jjj,jjj1+m) |
---|
| 1766 | do n=1,nend |
---|
| 1767 | nexp=nij-(n-1)*2 |
---|
| 1768 | nn=(m-1)*nend+n |
---|
| 1769 | fcc4(if2,ijnu,nn)= & |
---|
| 1770 | &fac3*fact(ii1,ii11+n)*fact(jj1,jj11+n)* & |
---|
| 1771 | &dble(nexp)*(one/dble(ii1)+one/dble(jj1))/two |
---|
| 1772 | if(mexp.ne.0) then |
---|
| 1773 | iti(if2,ijnu,nn,1)=-1 |
---|
| 1774 | iti(if2,ijnu,nn,2)=mexp |
---|
| 1775 | iti(if2,ijnu,nn,3)=nexp |
---|
| 1776 | else |
---|
| 1777 | iti(if2,ijnu,nn,1)=0 |
---|
| 1778 | iti(if2,ijnu,nn,2)=0 |
---|
| 1779 | iti(if2,ijnu,nn,3)=nexp |
---|
| 1780 | endif |
---|
| 1781 | enddo |
---|
| 1782 | enddo |
---|
| 1783 | endif |
---|
| 1784 | endif |
---|
| 1785 | if(k1.gt.0.and.ihoff.eq.0) then |
---|
| 1786 | !-- |
---|
| 1787 | !--Qx Term first Derivative by Ix for the Detuning by Ix |
---|
| 1788 | !-- |
---|
| 1789 | ido=1 |
---|
| 1790 | !--Finding Double Appearance |
---|
| 1791 | do ijnu1=1,ijnu-1 |
---|
| 1792 | if( & |
---|
| 1793 | &betl(ijnu1,1).eq.betl(ijnu,3).and. & |
---|
| 1794 | &betl(ijnu1,2).eq.betl(ijnu,4).and. & |
---|
| 1795 | &betl(ijnu1,3).eq.betl(ijnu,1).and. & |
---|
| 1796 | &betl(ijnu1,4).eq.betl(ijnu,2)) then |
---|
| 1797 | ido=0 |
---|
| 1798 | if2=if2+1 |
---|
| 1799 | if(if2.gt.2) call prror(1,13,if2) |
---|
| 1800 | fcc2(if2,ijnu1)=two*fcc2(if2,ijnu1) |
---|
| 1801 | endif |
---|
| 1802 | enddo |
---|
| 1803 | if(ido.eq.1) then |
---|
| 1804 | if2=if2+1 |
---|
| 1805 | if(if2.gt.2) call prror(1,13,if2) |
---|
| 1806 | ihw(if2,ijnu)=1 |
---|
| 1807 | ipl(if2,ijnu,1)=k1-1 |
---|
| 1808 | ipl(if2,ijnu,2)=l |
---|
| 1809 | betxp(if2,ijnu,1)=betl(ijnu,1) |
---|
| 1810 | betxp(if2,ijnu,2)=betl(ijnu,2) |
---|
| 1811 | betxp(if2,ijnu,3)=betl(ijnu,3) |
---|
| 1812 | betxp(if2,ijnu,4)=betl(ijnu,4) |
---|
| 1813 | fcc2(if2,ijnu)=fac1*dble(k1*iii*jjj)/ & |
---|
| 1814 | &dble(2**(k1+l1)) |
---|
| 1815 | !--Terms Iy**n where n>0 |
---|
| 1816 | if(ii1.ne.0.and.jj1.ne.0) then |
---|
| 1817 | do n=1,nend |
---|
| 1818 | nexp=nij-(n-1)*2 |
---|
| 1819 | fac3=fact(ii1,ii11+n)*fact(jj1,jj11+n) |
---|
| 1820 | do m=1,mend |
---|
| 1821 | mexp=mij-(m-1)*2 |
---|
| 1822 | mm=(n-1)*mend+m |
---|
| 1823 | fcc4(if2,ijnu,mm)= & |
---|
| 1824 | &fac3*fact(iii,iii1+m)*fact(jjj,jjj1+m)* & |
---|
| 1825 | &dble(mexp)*(one/dble(iii)+one/dble(jjj))/two |
---|
| 1826 | if(abs(fcc4(if2,ijnu,mm)).gt.pieni) then |
---|
| 1827 | if(nexp.ne.0) then |
---|
| 1828 | iti(if2,ijnu,mm,1)=1 |
---|
| 1829 | iti(if2,ijnu,mm,2)=mexp |
---|
| 1830 | iti(if2,ijnu,mm,3)=nexp |
---|
| 1831 | else if(mexp.ne.0) then |
---|
| 1832 | iti(if2,ijnu,mm,1)=0 |
---|
| 1833 | iti(if2,ijnu,mm,2)=mexp |
---|
| 1834 | iti(if2,ijnu,mm,3)=0 |
---|
| 1835 | endif |
---|
| 1836 | endif |
---|
| 1837 | enddo |
---|
| 1838 | enddo |
---|
| 1839 | !--Terms Iy**0=1 |
---|
| 1840 | else |
---|
| 1841 | do m=1,mend |
---|
| 1842 | mexp=mij-(m-1)*2 |
---|
| 1843 | fcc4(if2,ijnu,m)= & |
---|
| 1844 | &fact(iii,iii1+m)*fact(jjj,jjj1+m)* & |
---|
| 1845 | &dble(mexp)*(one/dble(iii)+one/dble(jjj))/two |
---|
| 1846 | if(abs(fcc4(if2,ijnu,m)).gt.pieni) then |
---|
| 1847 | iti(if2,ijnu,m,1)=0 |
---|
| 1848 | iti(if2,ijnu,m,2)=mexp |
---|
| 1849 | iti(if2,ijnu,m,3)=0 |
---|
| 1850 | endif |
---|
| 1851 | enddo |
---|
| 1852 | endif |
---|
| 1853 | endif |
---|
| 1854 | else if(ihoff.eq.1.and.ivoff.eq.0) then |
---|
| 1855 | !-- |
---|
| 1856 | !--Horizontal Term only one x, first Derivative by Iy and then Ix |
---|
| 1857 | !-- |
---|
| 1858 | ido=1 |
---|
| 1859 | !--Finding Double Appearance |
---|
| 1860 | do ijnu1=1,ijnu-1 |
---|
| 1861 | if( & |
---|
| 1862 | &betl(ijnu1,1).eq.betl(ijnu,3).and. & |
---|
| 1863 | &betl(ijnu1,2).eq.betl(ijnu,4).and. & |
---|
| 1864 | &betl(ijnu1,3).eq.betl(ijnu,1).and. & |
---|
| 1865 | &betl(ijnu1,4).eq.betl(ijnu,2)) then |
---|
| 1866 | ido=0 |
---|
| 1867 | if2=if2+1 |
---|
| 1868 | if(if2.gt.2) call prror(1,13,if2) |
---|
| 1869 | fcc2(if2,ijnu1)=two*fcc2(if2,ijnu1) |
---|
| 1870 | endif |
---|
| 1871 | enddo |
---|
| 1872 | if(ido.eq.1) then |
---|
| 1873 | if2=if2+1 |
---|
| 1874 | if(if2.gt.2) call prror(1,13,if2) |
---|
| 1875 | ihw(if2,ijnu)=1 |
---|
| 1876 | ipl(if2,ijnu,1)=k1 |
---|
| 1877 | ipl(if2,ijnu,2)=l1 |
---|
| 1878 | betxp(if2,ijnu,1)=betl(ijnu,1) |
---|
| 1879 | betxp(if2,ijnu,2)=betl(ijnu,2) |
---|
| 1880 | betxp(if2,ijnu,3)=betl(ijnu,3) |
---|
| 1881 | betxp(if2,ijnu,4)=betl(ijnu,4) |
---|
| 1882 | fcc2(if2,ijnu)=fac1*dble(k*ii1*jj1)/dble(2**(k1+l1)) |
---|
| 1883 | do n=1,nend |
---|
| 1884 | nexp=nij-(n-1)*2 |
---|
| 1885 | fcc4(if2,ijnu,n)= & |
---|
| 1886 | &fact(ii1,ii11+n)*fact(jj1,jj11+n)* & |
---|
| 1887 | &dble(nexp)*(one/dble(ii1)+one/dble(jj1))/two |
---|
| 1888 | iti(if2,ijnu,n,1)=0 |
---|
| 1889 | iti(if2,ijnu,n,2)=0 |
---|
| 1890 | iti(if2,ijnu,n,3)=nexp |
---|
| 1891 | enddo |
---|
| 1892 | endif |
---|
| 1893 | else if(iii.eq.1.and.jjj.eq.1) then |
---|
| 1894 | !-- |
---|
| 1895 | !--Vertical Term for uneven Multipoles |
---|
| 1896 | !-- |
---|
| 1897 | if2=if2+1 |
---|
| 1898 | if(if2.gt.2) call prror(1,13,if2) |
---|
| 1899 | ihw(if2,ijnu)=2 |
---|
| 1900 | ipl(if2,ijnu,1)=0 |
---|
| 1901 | ipl(if2,ijnu,2)=l-1 |
---|
| 1902 | betxp(if2,ijnu,1)=betl(ijnu,1) |
---|
| 1903 | betxp(if2,ijnu,2)=betl(ijnu,2) |
---|
| 1904 | betxp(if2,ijnu,3)=betl(ijnu,3) |
---|
| 1905 | betxp(if2,ijnu,4)=betl(ijnu,4) |
---|
| 1906 | fcc2(if2,ijnu)=fac1*dble(l*iii*jjj)/dble(2**(l-1)) |
---|
| 1907 | do n=1,nend |
---|
| 1908 | fac3=fact(ii1,ii11+n)*fact(jj1,jj11+n) |
---|
| 1909 | nexp=nij-(n-1)*2 |
---|
| 1910 | fcc4(if2,ijnu,n)=fac3* & |
---|
| 1911 | &(one/dble(iii)+one/dble(jjj))/two |
---|
| 1912 | if(nexp.ne.0) then |
---|
| 1913 | iti(if2,ijnu,n,1)=1 |
---|
| 1914 | iti(if2,ijnu,n,2)=1 |
---|
| 1915 | iti(if2,ijnu,n,3)=nexp |
---|
| 1916 | else |
---|
| 1917 | iti(if2,ijnu,n,1)=0 |
---|
| 1918 | iti(if2,ijnu,n,2)=1 |
---|
| 1919 | iti(if2,ijnu,n,3)=nexp |
---|
| 1920 | endif |
---|
| 1921 | enddo |
---|
| 1922 | endif |
---|
| 1923 | !-- |
---|
| 1924 | !--Vertical Term for even Multipoles |
---|
| 1925 | !-- |
---|
| 1926 | else if(ivoff.eq.0) then |
---|
| 1927 | if2=if2+1 |
---|
| 1928 | if(if2.gt.2) call prror(1,13,if2) |
---|
| 1929 | ihw(if2,ijnu)=2 |
---|
| 1930 | ipl(if2,ijnu,1)=0 |
---|
| 1931 | ipl(if2,ijnu,2)=l1-1 |
---|
| 1932 | betxp(if2,ijnu,1)=betl(ijnu,1) |
---|
| 1933 | betxp(if2,ijnu,2)=betl(ijnu,2) |
---|
| 1934 | betxp(if2,ijnu,3)=betl(ijnu,3) |
---|
| 1935 | betxp(if2,ijnu,4)=betl(ijnu,4) |
---|
| 1936 | fcc2(if2,ijnu)=fac1*dble(l1*ii1*jj1)/dble(2**(l-2)) |
---|
| 1937 | do n=1,nend |
---|
| 1938 | nexp=nij-(n-1)*2 |
---|
| 1939 | fcc4(if2,ijnu,n)= & |
---|
| 1940 | &fact(ii1,ii11+n)*fact(jj1,jj11+n)* & |
---|
| 1941 | &dble(nexp)*(one/dble(ii1)+one/dble(jj1))/two |
---|
| 1942 | iti(if2,ijnu,n,1)=0 |
---|
| 1943 | iti(if2,ijnu,n,2)=0 |
---|
| 1944 | iti(if2,ijnu,n,3)=nexp |
---|
| 1945 | enddo |
---|
| 1946 | endif |
---|
| 1947 | enddo |
---|
| 1948 | enddo |
---|
| 1949 | !--Write out Parameters properly |
---|
| 1950 | ic=0 |
---|
| 1951 | do nn=1,2 |
---|
| 1952 | do kk=1,mmultx |
---|
| 1953 | if(abs(fcc2(nn,kk)).gt.pieni) then |
---|
| 1954 | id=0 |
---|
| 1955 | ic=ic+1 |
---|
| 1956 | if(ic.gt.mmultx) call prror(1,6,mmultx) |
---|
| 1957 | ihv(ic)=ihw(nn,kk) |
---|
| 1958 | iplane(ic,1)=ipl(nn,kk,1) |
---|
| 1959 | iplane(ic,2)=ipl(nn,kk,2) |
---|
| 1960 | betexp(ic,1)=betxp(nn,kk,1) |
---|
| 1961 | betexp(ic,2)=betxp(nn,kk,2) |
---|
| 1962 | betexp(ic,3)=betxp(nn,kk,3) |
---|
| 1963 | betexp(ic,4)=betxp(nn,kk,4) |
---|
| 1964 | fac2(ic)=fcc2(nn,kk) |
---|
| 1965 | do mm=1,mmult |
---|
| 1966 | if(abs(fcc4(nn,kk,mm)).gt.pieni) then |
---|
| 1967 | id=id+1 |
---|
| 1968 | if(id.gt.mmult) call prror(1,7,mmult) |
---|
| 1969 | itij(ic,id,1)=iti(nn,kk,mm,1) |
---|
| 1970 | itij(ic,id,2)=iti(nn,kk,mm,2) |
---|
| 1971 | itij(ic,id,3)=iti(nn,kk,mm,3) |
---|
| 1972 | fac4(ic,id)=fcc4(nn,kk,mm) |
---|
| 1973 | endif |
---|
| 1974 | enddo |
---|
| 1975 | icd(ic)=id |
---|
| 1976 | endif |
---|
| 1977 | enddo |
---|
| 1978 | icc=ic |
---|
| 1979 | enddo |
---|
| 1980 | return |
---|
| 1981 | end subroutine caldet2 |
---|
| 1982 | subroutine caldt2(ia,ib,imuty) |
---|
| 1983 | !--------------------------------------------------------------------- |
---|
| 1984 | !-- |
---|
| 1985 | !--Calculation of the poisson bracket to derive the distortion function |
---|
| 1986 | !--Second order in the strength of the multipoles |
---|
| 1987 | !-- |
---|
| 1988 | !--------------------------------------------------------------------- |
---|
| 1989 | implicit none |
---|
| 1990 | integer ia,ib,ic,ica,icol1,icol2,icol3,icol4,icol5,icol6,id,ido1, & |
---|
| 1991 | &ido2,ido3,ido4,ihoff,ihoff1,ihoff2,ihvall,ii,ii1,iiend,iii,ijnu, & |
---|
| 1992 | &imuty,isiga,iti,ivoff,ivoff1,ivoff2,jj,jj1,jjend,jjj,k,k1,k2,k3, & |
---|
| 1993 | &k4,k41,k412s,k412s1,k412s2,k41s,k42,k42s,k43,k434s,k434s1,k434s2, & |
---|
| 1994 | &k44,k5,kk,l1,m,mh,mm,mmul,mmul2,mmult,mmultf,mmultw,mmultx,nblz, & |
---|
| 1995 | &nn,nnanf,nnend,nsig |
---|
| 1996 | double precision betl,diab,fac11,fac12,fac21,fac22,faca,facc, & |
---|
| 1997 | &facc1,facc2,fcc2,fcc4,pieni,signii,signjj |
---|
| 1998 | integer icc,icd,ifacd2,ifact4,ifacta,ihv,iplane,isig,itij,ityc, & |
---|
| 1999 | &iu_on,n1,n2 |
---|
| 2000 | double precision beta,betexp,bstr,cosav,det,det1,etl,etl1,etl2, & |
---|
| 2001 | &fac2,fac4,facd2,fact,fact0,factb,four,ham,one,phi,pi,pi2,pihi,qx, & |
---|
| 2002 | &qy,sbeta,sign,two,zero |
---|
| 2003 | real tlim,time0,time1,time2,time3,time4 |
---|
| 2004 | real tim1,tim2,tim3,tim4,tim5,tim6,tim7,tim8 |
---|
| 2005 | character*16 strn |
---|
| 2006 | character*18 comment |
---|
| 2007 | parameter(pieni=1d-34,nblz=4000,mmul=11,mh=100) |
---|
| 2008 | parameter(mmul2=mmul+1,mmult=8*mmul**2) |
---|
| 2009 | parameter(mmultw=2*mmul,mmultx=mmult/7,mmultf=2*mmultx*mmult) |
---|
| 2010 | common/c0/ pi,pi2,pihi,zero,one,two,four,comment(-mmul:mmul) |
---|
| 2011 | common/c1/ tlim,time0,time1,time2,time3,time4 |
---|
| 2012 | common/c2/ tim1,tim2,tim3,tim4,tim5,tim6,tim7,tim8 |
---|
| 2013 | common/c3/ bstr(-mmul:mmul,nblz),strn(-mmul:mmul,nblz) |
---|
| 2014 | common/c4/ beta(2,-mmul:mmul,nblz),sbeta(2,-mmul:mmul,nblz) |
---|
| 2015 | common/c4q/ phi(2,-mmul:mmul,nblz) |
---|
| 2016 | common/c5/ etl(-mmul:mmul,nblz),ityc(-mmul:mmul) |
---|
| 2017 | common/c6/ etl1,etl2,qx,qy,iu_on,n1,n2 |
---|
| 2018 | common/c7/ icc,icd(mmultx),itij(mmultx,mmult,6),isig(mmultx,mmult) |
---|
| 2019 | common/c8/ betexp(mmultx,4),fac4(mmultx,mmult),fac2(mmultx) |
---|
| 2020 | common/c9/ fact(0:mmul,0:mmul2),sign(4,4),cosav(0:mmul) |
---|
| 2021 | common/c10/ ifacd2(4,mmultf),facd2(2,mmultf),ham(2,mmultf) |
---|
| 2022 | common/c11/ det1(2,0:mmul),det(2,-mmul:mmul,0:mmul,0:mmul) |
---|
| 2023 | common/c12/ ihv(mmultx),iplane(mmultx,2) |
---|
| 2024 | common/c13/ fact0(-mmul:mmul,mh),factb(2,-mmul:mmul,mh) |
---|
| 2025 | common/c14/ ifacta(2,-mmul:mmul,mh),ifact4(4,-mmul:mmul,mh) |
---|
| 2026 | dimension iti(2,mmultx,mmult,6),isiga(2,mmultx,mmult) |
---|
| 2027 | dimension fcc2(2,mmultx),fcc4(2,mmultx,mmult) |
---|
| 2028 | dimension betl(mmultx,4) |
---|
| 2029 | !--------------------------------------------------------------------- |
---|
| 2030 | do k=1,mmultx |
---|
| 2031 | fcc2(1,k)=zero |
---|
| 2032 | fcc2(2,k)=zero |
---|
| 2033 | do m=1,mmult |
---|
| 2034 | do nn=1,2 |
---|
| 2035 | iti(nn,k,m,1)=0 |
---|
| 2036 | iti(nn,k,m,2)=0 |
---|
| 2037 | iti(nn,k,m,3)=0 |
---|
| 2038 | iti(nn,k,m,4)=0 |
---|
| 2039 | iti(nn,k,m,5)=0 |
---|
| 2040 | iti(nn,k,m,6)=0 |
---|
| 2041 | isiga(nn,k,m)=0 |
---|
| 2042 | betl(k,1)=zero |
---|
| 2043 | betl(k,2)=zero |
---|
| 2044 | betl(k,3)=zero |
---|
| 2045 | betl(k,4)=zero |
---|
| 2046 | fcc4(nn,k,m)=zero |
---|
| 2047 | enddo |
---|
| 2048 | enddo |
---|
| 2049 | enddo |
---|
| 2050 | diab=dble(ia+ib)/two |
---|
| 2051 | facc=one/four/(two**diab)/dble(ia*ib) |
---|
| 2052 | !--Prefactor I (General factors) |
---|
| 2053 | iiend=(ia+2)/2 |
---|
| 2054 | jjend=(ib+2)/2 |
---|
| 2055 | !--(iiend,jjend)=>Total number of terms in (x+iy)**n |
---|
| 2056 | !--(i,j) is on left or right side of the Poisson bracket respectively |
---|
| 2057 | signii=-one |
---|
| 2058 | if((imuty.eq.4.or.imuty.eq.2).and.2*iiend.eq.ia+2) iiend=iiend-1 |
---|
| 2059 | do ii=1,iiend |
---|
| 2060 | ii1=(ii-1)*2 |
---|
| 2061 | if(imuty.eq.4.or.imuty.eq.2) ii1=ii1+1 |
---|
| 2062 | signii=signii*(-one) |
---|
| 2063 | signjj=-one |
---|
| 2064 | if((imuty.eq.3.or.imuty.eq.2).and.2*jjend.eq.ib+2) jjend=jjend-1 |
---|
| 2065 | do jj=1,jjend |
---|
| 2066 | jj1=(jj-1)*2 |
---|
| 2067 | !--(ii1,jj1)=>Increment of powers (0->{(iiend,jjend)-1}*2) |
---|
| 2068 | !-- goes in steps of 2 |
---|
| 2069 | !--These are the exponents of the second variable y |
---|
| 2070 | if(imuty.eq.3.or.imuty.eq.2) jj1=jj1+1 |
---|
| 2071 | signjj=signjj*(-one) |
---|
| 2072 | facc1=facc*signii*signjj*fact(ia,1+ii1)*fact(ib,1+jj1) |
---|
| 2073 | !--Prefactor II (signs, binominal coefficients) |
---|
| 2074 | iii=ia-ii1 |
---|
| 2075 | jjj=ib-jj1 |
---|
| 2076 | !--(iii,jjj)=>Exponent of the first variable x |
---|
| 2077 | ijnu=(ii-1)*jjend+jj |
---|
| 2078 | if(ijnu.gt.mmultx) call prror(3,6,mmultx) |
---|
| 2079 | !--inju=># of current case treated of a total of (iiend*jjend) |
---|
| 2080 | betl(ijnu,1)=dble(iii)/two |
---|
| 2081 | betl(ijnu,2)=dble(ii1)/two |
---|
| 2082 | betl(ijnu,3)=dble(jjj)/two |
---|
| 2083 | betl(ijnu,4)=dble(jj1)/two |
---|
| 2084 | ihoff=0 |
---|
| 2085 | ihoff1=0 |
---|
| 2086 | ihoff2=0 |
---|
| 2087 | ivoff=0 |
---|
| 2088 | ivoff1=0 |
---|
| 2089 | ivoff2=0 |
---|
| 2090 | if(iii.eq.0.and.jjj.eq.0) ihoff=1 |
---|
| 2091 | if(iii.eq.0.and.jjj.ne.0.and.ihoff.eq.0) ihoff1=1 |
---|
| 2092 | if(iii.ne.0.and.jjj.eq.0.and.ihoff.eq.0) ihoff2=1 |
---|
| 2093 | if(ii1.eq.0.and.jj1.eq.0) ivoff=1 |
---|
| 2094 | if(ii1.eq.0.and.jj1.ne.0.and.ivoff.eq.0) ivoff1=1 |
---|
| 2095 | if(ii1.ne.0.and.jj1.eq.0.and.ivoff.eq.0) ivoff2=1 |
---|
| 2096 | ihvall=ihoff+ihoff1+ihoff2+ivoff+ivoff1+ivoff2 |
---|
| 2097 | !--(ihoff,ivoff)=0/1 (x,y) on (at least on one side/on no sides) of |
---|
| 2098 | !-- Poisson bracket |
---|
| 2099 | !--(ihoff1,ivoff1)=0/1 (x,y) on (both sides/on left side) of |
---|
| 2100 | !-- Poisson bracket |
---|
| 2101 | !--(ihoff2,ivoff2)=0/1 (x,y) on (both sides/on right side) of |
---|
| 2102 | !-- Poisson bracket |
---|
| 2103 | facc2=-facc1/(two**(diab-one))/four |
---|
| 2104 | !--Prefactor III (Invariant->emittance) and -1/4d0 from evaluation of |
---|
| 2105 | !-- trigonometric functions |
---|
| 2106 | !--------------------------------------------------------------------- |
---|
| 2107 | !-- |
---|
| 2108 | !--Deriving the Derivatives |
---|
| 2109 | !-- |
---|
| 2110 | !--------------------------------------------------------------------- |
---|
| 2111 | nnanf=1 |
---|
| 2112 | nnend=2 |
---|
| 2113 | if(ihoff.eq.1.or.ihoff1.eq.1.or.ihoff2.eq.1) then |
---|
| 2114 | nnanf=2 |
---|
| 2115 | nnend=2 |
---|
| 2116 | endif |
---|
| 2117 | if(ivoff.eq.1.or.ivoff1.eq.1.or.ivoff2.eq.1) nnend=1 |
---|
| 2118 | if(ihvall.eq.0) then |
---|
| 2119 | nnanf=1 |
---|
| 2120 | nnend=2 |
---|
| 2121 | endif |
---|
| 2122 | !--Differentiation nn=1 => horizontal nn=2 => vertical |
---|
| 2123 | !-- ==> exchange planes (ido1,ido2 => ido3,ido4) |
---|
| 2124 | !-- ==> exchange terms of secondary resonances (icol1 => icol2) |
---|
| 2125 | !-- ==> exchange primary resonance terms (icol3,icol4 => icol5,icol6) |
---|
| 2126 | do nn=nnanf,nnend |
---|
| 2127 | if(nn.eq.1) then |
---|
| 2128 | ido1=iii |
---|
| 2129 | ido2=jjj |
---|
| 2130 | ido3=ii1 |
---|
| 2131 | ido4=jj1 |
---|
| 2132 | icol1=1 |
---|
| 2133 | icol2=2 |
---|
| 2134 | icol3=3 |
---|
| 2135 | icol4=4 |
---|
| 2136 | icol5=5 |
---|
| 2137 | icol6=6 |
---|
| 2138 | !--nsig see below |
---|
| 2139 | nsig=3 |
---|
| 2140 | endif |
---|
| 2141 | if(nn.eq.2) then |
---|
| 2142 | ido1=ii1 |
---|
| 2143 | ido2=jj1 |
---|
| 2144 | ido3=iii |
---|
| 2145 | ido4=jjj |
---|
| 2146 | icol1=2 |
---|
| 2147 | icol2=1 |
---|
| 2148 | icol3=5 |
---|
| 2149 | icol4=6 |
---|
| 2150 | icol5=3 |
---|
| 2151 | icol6=4 |
---|
| 2152 | !--nsig is set to 4 for the vertical differentiation as the |
---|
| 2153 | !--horizontal stays always positive |
---|
| 2154 | nsig=4 |
---|
| 2155 | endif |
---|
| 2156 | fcc2(nn,ijnu)=facc2*dble(ido1*ido2) |
---|
| 2157 | ica=0 |
---|
| 2158 | !--(k1,k2) always differentiation nn=1 horizontal nn=2 vertical |
---|
| 2159 | do k1=0,ido1,2 |
---|
| 2160 | k41=ido1-k1 |
---|
| 2161 | !--k41 see k44 |
---|
| 2162 | fac11=dble(k41)/dble(ido1)*fact(ido1,1+k1/2) |
---|
| 2163 | fac12=fact(ido1,1+k1/2) |
---|
| 2164 | do k2=0,ido2,2 |
---|
| 2165 | k42=ido2-k2 |
---|
| 2166 | !--k42 see k44 |
---|
| 2167 | fac21=fact(ido2,1+k2/2) |
---|
| 2168 | fac22=dble(k42)/dble(ido2)*fact(ido2,1+k2/2) |
---|
| 2169 | !--(l1) sums over the 3-4 possible sign combinations |
---|
| 2170 | do l1=1,nsig |
---|
| 2171 | faca=sign(l1,1)*fac11*fac21- & |
---|
| 2172 | &sign(l1,2)*fac12*fac22 |
---|
| 2173 | if(k41.eq.0) faca=-sign(l1,2)*fac12*fac22 |
---|
| 2174 | if(k42.eq.0) faca=sign(l1,1)*fac11*fac21 |
---|
| 2175 | !--Condition 1 & 5 clear |
---|
| 2176 | !--2 Condition avoids duplication of l1 by l2 |
---|
| 2177 | !--3 Condition avoids duplication of l1 by l3 |
---|
| 2178 | !--4 Condition avoids duplication of l1 with any other |
---|
| 2179 | if((abs(faca).gt.pieni).and. & |
---|
| 2180 | &(l1.ne.2.or.k41.ne.0).and. & |
---|
| 2181 | &(l1.ne.3.or.k42.ne.0).and. & |
---|
| 2182 | &(l1.ne.4.or.(k41.ne.0.and.k42.ne.0)).and. & |
---|
| 2183 | &((k41+k42).ge.1)) then |
---|
| 2184 | !--(k3,k4) plane with no differentiation nn=1 vertical nn=2 horizontal |
---|
| 2185 | !--((ido3+3),(ido4+3)) allow for constant term in the |
---|
| 2186 | !--undifferentiated part |
---|
| 2187 | do k3=0,2*ido3,2 |
---|
| 2188 | do k4=0,2*ido4,2 |
---|
| 2189 | !--k43,k44 control the secondary resonance |
---|
| 2190 | !--i.e. the same that appear in first order |
---|
| 2191 | k43=ido3-k3 |
---|
| 2192 | k44=ido4-k4 |
---|
| 2193 | !--(k41s,k42s) takes care of the sign of k41 and k42 |
---|
| 2194 | k41s=sign(l1,1)*k41 |
---|
| 2195 | k42s=sign(l1,2)*k42 |
---|
| 2196 | !--(k412s) sum of the 2 terms with differentiation |
---|
| 2197 | k412s=k41s+k42s |
---|
| 2198 | k412s1=(ido1+ido2+k412s)/2-1 |
---|
| 2199 | k412s2=k412s1-k412s |
---|
| 2200 | !--(k434s) sum of the 2 terms without differentiation |
---|
| 2201 | k434s=k43+k44 |
---|
| 2202 | k434s1=(ido3+ido4+k434s)/2 |
---|
| 2203 | k434s2=k434s1-k434s |
---|
| 2204 | !--1. (0,0) Secondary Resonance has to be excluded |
---|
| 2205 | !--2. Select resonances in the desired nomenclature |
---|
| 2206 | !--at least the sum must be larger than 0 |
---|
| 2207 | !--and the resonance should be written in descending order starting |
---|
| 2208 | !--with the horizontal one |
---|
| 2209 | if( & |
---|
| 2210 | &(k42s.ne.0.or.k44.ne.0).and. & |
---|
| 2211 | &(k412s.gt.0.or.k434s.gt.0).and. & |
---|
| 2212 | &((nn.eq.1.and.k412s1.ge.k412s2).or. & |
---|
| 2213 | &(nn.eq.2.and.k434s1.ge.k434s2))) then |
---|
| 2214 | !--Program check: finds duplicated cases |
---|
| 2215 | do k5=1,ica |
---|
| 2216 | if( & |
---|
| 2217 | &iti(nn,ijnu,k5,icol1).eq.k42s.and. & |
---|
| 2218 | &iti(nn,ijnu,k5,icol2).eq.k44.and. & |
---|
| 2219 | &iti(nn,ijnu,k5,icol3).eq.k412s1.and. & |
---|
| 2220 | &iti(nn,ijnu,k5,icol4).eq.k412s2.and. & |
---|
| 2221 | &iti(nn,ijnu,k5,icol5).eq.k434s1.and. & |
---|
| 2222 | &iti(nn,ijnu,k5,icol6).eq.k434s2 & |
---|
| 2223 | &) call prror(3,9,2) |
---|
| 2224 | enddo |
---|
| 2225 | ica=ica+1 |
---|
| 2226 | if(ica.gt.mmult) call prror(3,5,mmult) |
---|
| 2227 | fcc4(nn,ijnu,ica)=faca* & |
---|
| 2228 | &fact(ido3,1+k3/2)*fact(ido4,1+k4/2) |
---|
| 2229 | iti(nn,ijnu,ica,icol1)=k42s |
---|
| 2230 | iti(nn,ijnu,ica,icol2)=k44 |
---|
| 2231 | iti(nn,ijnu,ica,icol3)=k412s1 |
---|
| 2232 | iti(nn,ijnu,ica,icol4)=k412s2 |
---|
| 2233 | iti(nn,ijnu,ica,icol5)=k434s1 |
---|
| 2234 | iti(nn,ijnu,ica,icol6)=k434s2 |
---|
| 2235 | if(k41.eq.0.and.k43.eq.0) isiga(nn,ijnu,ica)=1 |
---|
| 2236 | endif |
---|
| 2237 | enddo |
---|
| 2238 | enddo |
---|
| 2239 | endif |
---|
| 2240 | enddo |
---|
| 2241 | enddo |
---|
| 2242 | enddo |
---|
| 2243 | enddo |
---|
| 2244 | enddo |
---|
| 2245 | enddo |
---|
| 2246 | !--Write out Parameters properly |
---|
| 2247 | ic=0 |
---|
| 2248 | do nn=1,2 |
---|
| 2249 | do kk=1,mmultx |
---|
| 2250 | if(abs(fcc2(nn,kk)).gt.pieni) then |
---|
| 2251 | id=0 |
---|
| 2252 | ic=ic+1 |
---|
| 2253 | if(ic.gt.mmultx) call prror(3,6,mmultx) |
---|
| 2254 | betexp(ic,1)=betl(kk,1) |
---|
| 2255 | betexp(ic,2)=betl(kk,2) |
---|
| 2256 | betexp(ic,3)=betl(kk,3) |
---|
| 2257 | betexp(ic,4)=betl(kk,4) |
---|
| 2258 | fac2(ic)=fcc2(nn,kk) |
---|
| 2259 | do mm=1,mmult |
---|
| 2260 | if(abs(fcc4(nn,kk,mm)).gt.pieni) then |
---|
| 2261 | id=id+1 |
---|
| 2262 | if(id.gt.mmult) call prror(3,7,mmult) |
---|
| 2263 | if(isiga(nn,kk,mm).eq.1) then |
---|
| 2264 | isig(ic,id)=1 |
---|
| 2265 | isiga(nn,kk,mm)=0 |
---|
| 2266 | endif |
---|
| 2267 | itij(ic,id,1)=iti(nn,kk,mm,1) |
---|
| 2268 | itij(ic,id,2)=iti(nn,kk,mm,2) |
---|
| 2269 | itij(ic,id,3)=iti(nn,kk,mm,3) |
---|
| 2270 | itij(ic,id,4)=iti(nn,kk,mm,4) |
---|
| 2271 | itij(ic,id,5)=iti(nn,kk,mm,5) |
---|
| 2272 | itij(ic,id,6)=iti(nn,kk,mm,6) |
---|
| 2273 | fac4(ic,id)=fcc4(nn,kk,mm) |
---|
| 2274 | endif |
---|
| 2275 | enddo |
---|
| 2276 | icd(ic)=id |
---|
| 2277 | endif |
---|
| 2278 | enddo |
---|
| 2279 | icc=ic |
---|
| 2280 | enddo |
---|
| 2281 | return |
---|
| 2282 | end subroutine caldt2 |
---|
| 2283 | subroutine init |
---|
| 2284 | !--------------------------------------------------------------------- |
---|
| 2285 | !--Initilisation |
---|
| 2286 | !--------------------------------------------------------------------- |
---|
| 2287 | implicit none |
---|
| 2288 | integer i,ii,j,k,mh,mmul,mmul2,mmult,mmultf,mmultw,mmultx,nblz |
---|
| 2289 | double precision pieni |
---|
| 2290 | integer icc,icd,ifacd2,ifact4,ifacta,ihv,iplane,isig,itij,ityc, & |
---|
| 2291 | &iu_on,n1,n2 |
---|
| 2292 | double precision beta,betexp,bstr,cosav,det,det1,etl,etl1,etl2, & |
---|
| 2293 | &fac2,fac4,facd2,fact,fact0,factb,four,ham,one,phi,pi,pi2,pihi,qx, & |
---|
| 2294 | &qy,sbeta,sign,two,zero |
---|
| 2295 | real tlim,time0,time1,time2,time3,time4 |
---|
| 2296 | real tim1,tim2,tim3,tim4,tim5,tim6,tim7,tim8 |
---|
| 2297 | character*16 strn |
---|
| 2298 | character*18 comment |
---|
| 2299 | parameter(pieni=1d-34,nblz=4000,mmul=11,mh=100) |
---|
| 2300 | parameter(mmul2=mmul+1,mmult=8*mmul**2) |
---|
| 2301 | parameter(mmultw=2*mmul,mmultx=mmult/7,mmultf=2*mmultx*mmult) |
---|
| 2302 | common/c0/ pi,pi2,pihi,zero,one,two,four,comment(-mmul:mmul) |
---|
| 2303 | common/c1/ tlim,time0,time1,time2,time3,time4 |
---|
| 2304 | common/c2/ tim1,tim2,tim3,tim4,tim5,tim6,tim7,tim8 |
---|
| 2305 | common/c3/ bstr(-mmul:mmul,nblz),strn(-mmul:mmul,nblz) |
---|
| 2306 | common/c4/ beta(2,-mmul:mmul,nblz),sbeta(2,-mmul:mmul,nblz) |
---|
| 2307 | common/c4q/ phi(2,-mmul:mmul,nblz) |
---|
| 2308 | common/c5/ etl(-mmul:mmul,nblz),ityc(-mmul:mmul) |
---|
| 2309 | common/c6/ etl1,etl2,qx,qy,iu_on,n1,n2 |
---|
| 2310 | common/c7/ icc,icd(mmultx),itij(mmultx,mmult,6),isig(mmultx,mmult) |
---|
| 2311 | common/c8/ betexp(mmultx,4),fac4(mmultx,mmult),fac2(mmultx) |
---|
| 2312 | common/c9/ fact(0:mmul,0:mmul2),sign(4,4),cosav(0:mmul) |
---|
| 2313 | common/c10/ ifacd2(4,mmultf),facd2(2,mmultf),ham(2,mmultf) |
---|
| 2314 | common/c11/ det1(2,0:mmul),det(2,-mmul:mmul,0:mmul,0:mmul) |
---|
| 2315 | common/c12/ ihv(mmultx),iplane(mmultx,2) |
---|
| 2316 | common/c13/ fact0(-mmul:mmul,mh),factb(2,-mmul:mmul,mh) |
---|
| 2317 | common/c14/ ifacta(2,-mmul:mmul,mh),ifact4(4,-mmul:mmul,mh) |
---|
| 2318 | !--------------------------------------------------------------------- |
---|
| 2319 | tlim=1e7 |
---|
| 2320 | zero=dble(0d0) |
---|
| 2321 | one=1d0 |
---|
| 2322 | two=2d0 |
---|
| 2323 | four=4d0 |
---|
| 2324 | pi=atan(one)*four |
---|
| 2325 | pi2=pi*two |
---|
| 2326 | pihi=two/pi |
---|
| 2327 | comment(2)=' Erect Quadrupoles' |
---|
| 2328 | comment(-2)=' Skew Quadrupoles ' |
---|
| 2329 | comment(3)=' Erect Sextupoles ' |
---|
| 2330 | comment(-3)=' Skew Sextupoles ' |
---|
| 2331 | comment(4)=' Erect Octupoles ' |
---|
| 2332 | comment(-4)=' Skew Octupoles ' |
---|
| 2333 | comment(5)=' Erect Decapoles ' |
---|
| 2334 | comment(-5)=' Skew Decapoles ' |
---|
| 2335 | comment(6)=' Erect Dodecapoles' |
---|
| 2336 | comment(-6)=' Skew Dodecapoles ' |
---|
| 2337 | comment(7)=' Erect 14th-Pole ' |
---|
| 2338 | comment(-7)=' Skew 14th-Pole ' |
---|
| 2339 | comment(8)=' Erect 16th-Pole ' |
---|
| 2340 | comment(-8)=' Skew 16th-Pole ' |
---|
| 2341 | comment(9)=' Erect 18th-Pole ' |
---|
| 2342 | comment(-9)=' Skew 18th-Pole ' |
---|
| 2343 | comment(10)=' Erect 20th-Pole ' |
---|
| 2344 | comment(-10)=' Skew 20th-Pole ' |
---|
| 2345 | comment(11)=' Erect 22th-Pole ' |
---|
| 2346 | comment(-11)=' Skew 22th-Pole ' |
---|
| 2347 | !-- |
---|
| 2348 | !--Binominals |
---|
| 2349 | !-- |
---|
| 2350 | cosav(0)=one |
---|
| 2351 | do i=2,mmul,2 |
---|
| 2352 | ii=i/2 |
---|
| 2353 | cosav(ii)=cosav(ii-1)*dble(i-1)/dble(i) |
---|
| 2354 | enddo |
---|
| 2355 | do j=0,mmul |
---|
| 2356 | do i=0,mmul2 |
---|
| 2357 | fact(j,i)=one |
---|
| 2358 | enddo |
---|
| 2359 | enddo |
---|
| 2360 | do i=2,mmul |
---|
| 2361 | do j=2,mmul |
---|
| 2362 | if(j.ge.i) then |
---|
| 2363 | fact(j,i)=fact(j-1,i)+fact(j-1,i-1) |
---|
| 2364 | endif |
---|
| 2365 | enddo |
---|
| 2366 | enddo |
---|
| 2367 | do i=-mmul,mmul |
---|
| 2368 | do j=0,mmul |
---|
| 2369 | do k=0,mmul |
---|
| 2370 | det(1,i,j,k)=zero |
---|
| 2371 | det(2,i,j,k)=zero |
---|
| 2372 | enddo |
---|
| 2373 | enddo |
---|
| 2374 | enddo |
---|
| 2375 | do i=0,mmul |
---|
| 2376 | det1(1,i)=zero |
---|
| 2377 | det1(2,i)=zero |
---|
| 2378 | enddo |
---|
| 2379 | !-- |
---|
| 2380 | !--sign array |
---|
| 2381 | !-- |
---|
| 2382 | sign(1,1)=one |
---|
| 2383 | sign(1,2)=one |
---|
| 2384 | sign(2,1)=-one |
---|
| 2385 | sign(2,2)=one |
---|
| 2386 | sign(3,1)=one |
---|
| 2387 | sign(3,2)=-one |
---|
| 2388 | sign(4,1)=-one |
---|
| 2389 | sign(4,2)=-one |
---|
| 2390 | return |
---|
| 2391 | end subroutine init |
---|
| 2392 | subroutine readdat |
---|
| 2393 | !--------------------------------------------------------------------- |
---|
| 2394 | !-- |
---|
| 2395 | !--Reading Lattice Information from File # 34 |
---|
| 2396 | !-- |
---|
| 2397 | !--------------------------------------------------------------------- |
---|
| 2398 | implicit none |
---|
| 2399 | integer i,ich1,ich2,ii,isize,itype,j,jj,mh,mmul,mmul2,mmult, & |
---|
| 2400 | &mmultf,mmultw,mmultx,nblz,ndum,num |
---|
| 2401 | double precision betx,bety,et,phix0,phiy0,pieni,str |
---|
| 2402 | character*16 name |
---|
| 2403 | character*200 ch |
---|
| 2404 | character*205 ch1 |
---|
| 2405 | integer icc,icd,ifacd2,ifact4,ifacta,ihv,iplane,isig,itij,ityc, & |
---|
| 2406 | &iu_on,n1,n2 |
---|
| 2407 | double precision beta,betexp,bstr,cosav,det,det1,etl,etl1,etl2, & |
---|
| 2408 | &fac2,fac4,facd2,fact,fact0,factb,four,ham,one,phi,pi,pi2,pihi,qx, & |
---|
| 2409 | &qy,sbeta,sign,two,zero, factor |
---|
| 2410 | real tlim,time0,time1,time2,time3,time4 |
---|
| 2411 | real tim1,tim2,tim3,tim4,tim5,tim6,tim7,tim8 |
---|
| 2412 | character*16 strn |
---|
| 2413 | character*18 comment |
---|
| 2414 | parameter(pieni=1d-34,nblz=4000,mmul=11,mh=100) |
---|
| 2415 | parameter(mmul2=mmul+1,mmult=8*mmul**2) |
---|
| 2416 | parameter(mmultw=2*mmul,mmultx=mmult/7,mmultf=2*mmultx*mmult) |
---|
| 2417 | common/c0/ pi,pi2,pihi,zero,one,two,four,comment(-mmul:mmul) |
---|
| 2418 | common/c1/ tlim,time0,time1,time2,time3,time4 |
---|
| 2419 | common/c2/ tim1,tim2,tim3,tim4,tim5,tim6,tim7,tim8 |
---|
| 2420 | common/c3/ bstr(-mmul:mmul,nblz),strn(-mmul:mmul,nblz) |
---|
| 2421 | common/c4/ beta(2,-mmul:mmul,nblz),sbeta(2,-mmul:mmul,nblz) |
---|
| 2422 | common/c4q/ phi(2,-mmul:mmul,nblz) |
---|
| 2423 | common/c5/ etl(-mmul:mmul,nblz),ityc(-mmul:mmul) |
---|
| 2424 | common/c6/ etl1,etl2,qx,qy,iu_on,n1,n2 |
---|
| 2425 | common/c7/ icc,icd(mmultx),itij(mmultx,mmult,6),isig(mmultx,mmult) |
---|
| 2426 | common/c8/ betexp(mmultx,4),fac4(mmultx,mmult),fac2(mmultx) |
---|
| 2427 | common/c9/ fact(0:mmul,0:mmul2),sign(4,4),cosav(0:mmul) |
---|
| 2428 | common/c10/ ifacd2(4,mmultf),facd2(2,mmultf),ham(2,mmultf) |
---|
| 2429 | common/c11/ det1(2,0:mmul),det(2,-mmul:mmul,0:mmul,0:mmul) |
---|
| 2430 | common/c12/ ihv(mmultx),iplane(mmultx,2) |
---|
| 2431 | common/c13/ fact0(-mmul:mmul,mh),factb(2,-mmul:mmul,mh) |
---|
| 2432 | common/c14/ ifacta(2,-mmul:mmul,mh),ifact4(4,-mmul:mmul,mh) |
---|
| 2433 | !--------------------------------------------------------------------- |
---|
| 2434 | |
---|
| 2435 | do i=1,nblz |
---|
| 2436 | do j=-mmul,mmul |
---|
| 2437 | ityc(j)=0 |
---|
| 2438 | bstr(j,i)=zero |
---|
| 2439 | strn(j,i)=" " |
---|
| 2440 | beta(1,j,i)=zero |
---|
| 2441 | beta(2,j,i)=zero |
---|
| 2442 | sbeta(1,j,i)=zero |
---|
| 2443 | sbeta(2,j,i)=zero |
---|
| 2444 | phi(1,j,i)=zero |
---|
| 2445 | phi(2,j,i)=zero |
---|
| 2446 | etl(j,i)=zero |
---|
| 2447 | enddo |
---|
| 2448 | enddo |
---|
| 2449 | do ii=1,10000000 |
---|
| 2450 | read(34,'(A150)',end=900) ch |
---|
| 2451 | ich1=0 |
---|
| 2452 | ich2=0 |
---|
| 2453 | do jj=1,200 |
---|
| 2454 | if(ich1.eq.0.and.ch(jj:jj).ne.' ') ich1=1 |
---|
| 2455 | if(ich1.eq.1.and.ch(jj:jj).eq.' ') ich1=2 |
---|
| 2456 | if(ich1.eq.2.and.ch(jj:jj).ne.' ') then |
---|
| 2457 | ich1=3 |
---|
| 2458 | ich2=jj |
---|
| 2459 | endif |
---|
| 2460 | if(ich1.eq.3.and.ch(jj:jj).eq.' ') then |
---|
| 2461 | ch1(1:200)=ch(1:ich2-1)//''''//ch(ich2:jj-1)//''''// & |
---|
| 2462 | &ch(jj:200)//' / ' |
---|
| 2463 | goto 210 |
---|
| 2464 | endif |
---|
| 2465 | enddo |
---|
| 2466 | 210 continue |
---|
| 2467 | read(ch1,*) et,name,itype,str,betx,bety,phix0,phiy0 |
---|
| 2468 | if(itype.eq.100) then |
---|
| 2469 | qx=phix0*pi |
---|
| 2470 | qy=phiy0*pi |
---|
| 2471 | goto 5 |
---|
| 2472 | endif |
---|
| 2473 | if(et.ge.etl1.and.et.le.etl2) then |
---|
| 2474 | if(abs(str).ge.pieni.and.itype.ge.n1.and.itype.le.n2) then |
---|
| 2475 | ityc(itype)=ityc(itype)+1 |
---|
| 2476 | isize=ityc(itype) |
---|
| 2477 | if(isize.gt.nblz) call prror(3,3,nblz) |
---|
| 2478 | etl(itype,isize)=et |
---|
| 2479 | strn(itype,isize)=name |
---|
| 2480 | factor = 10**(3*(iabs(itype)-2)) |
---|
| 2481 | bstr(itype,isize)=str*factor |
---|
| 2482 | beta(1,itype,isize)=betx |
---|
| 2483 | beta(2,itype,isize)=bety |
---|
| 2484 | sbeta(1,itype,isize)=dsqrt(betx) |
---|
| 2485 | sbeta(2,itype,isize)=dsqrt(bety) |
---|
| 2486 | phi(1,itype,isize)=phix0*two*pi |
---|
| 2487 | phi(2,itype,isize)=phiy0*two*pi |
---|
| 2488 | endif |
---|
| 2489 | endif |
---|
| 2490 | enddo |
---|
| 2491 | call prror(3,4,num) |
---|
| 2492 | 5 continue |
---|
| 2493 | write(6,*) ' Qx= ',qx/pi,' Qy= ',qy/pi |
---|
| 2494 | close(34) |
---|
| 2495 | return |
---|
| 2496 | 900 close(34) |
---|
| 2497 | call prror(3,1,ndum) |
---|
| 2498 | end subroutine readdat |
---|
| 2499 | subroutine detwri(icase,imu,ih1) |
---|
| 2500 | !--------------------------------------------------------------------- |
---|
| 2501 | !--Organize Writing for Detune1 |
---|
| 2502 | !--------------------------------------------------------------------- |
---|
| 2503 | implicit none |
---|
| 2504 | integer i,ic,icase,id,ih1,imu,j,k,mh,mmul,mmul2,mmult,mmultf, & |
---|
| 2505 | &mmultw,mmultx,nblz |
---|
| 2506 | double precision pieni |
---|
| 2507 | integer icc,icd,ifacd2,ifact4,ifacta,ihv,iplane,isig,itij,ityc, & |
---|
| 2508 | &iu_on,n1,n2 |
---|
| 2509 | double precision beta,betexp,bstr,cosav,det,det1,etl,etl1,etl2, & |
---|
| 2510 | &fac2,fac4,facd2,fact,fact0,factb,four,ham,one,phi,pi,pi2, & |
---|
| 2511 | &pihi,qx,qy,sbeta,sign,two,zero |
---|
| 2512 | real tlim,time0,time1,time2,time3,time4 |
---|
| 2513 | real tim1,tim2,tim3,tim4,tim5,tim6,tim7,tim8 |
---|
| 2514 | integer j70,j71,j72,j73,j74,j75,j76,j77,j78,j79 |
---|
| 2515 | integer table_size_70,table_size_71,table_size_72,table_size_73, & |
---|
| 2516 | &table_size_74,table_size_75,table_size_76,table_size_77, & |
---|
| 2517 | &table_size_78,table_size_79 |
---|
| 2518 | character*16 strn,table_name |
---|
| 2519 | character*18 comment |
---|
| 2520 | parameter(pieni=1d-34,nblz=4000,mmul=11,mh=100) |
---|
| 2521 | parameter(mmul2=mmul+1,mmult=8*mmul**2) |
---|
| 2522 | parameter(mmultw=2*mmul,mmultx=mmult/7,mmultf=2*mmultx*mmult) |
---|
| 2523 | common/c0/ pi,pi2,pihi,zero,one,two,four,comment(-mmul:mmul) |
---|
| 2524 | common/c1/ tlim,time0,time1,time2,time3,time4 |
---|
| 2525 | common/c2/ tim1,tim2,tim3,tim4,tim5,tim6,tim7,tim8 |
---|
| 2526 | common/c3/ bstr(-mmul:mmul,nblz),strn(-mmul:mmul,nblz) |
---|
| 2527 | common/c4/ beta(2,-mmul:mmul,nblz),sbeta(2,-mmul:mmul,nblz) |
---|
| 2528 | common/c4q/ phi(2,-mmul:mmul,nblz) |
---|
| 2529 | common/c5/ etl(-mmul:mmul,nblz),ityc(-mmul:mmul) |
---|
| 2530 | common/c6/ etl1,etl2,qx,qy,iu_on,n1,n2 |
---|
| 2531 | common/c7/ icc,icd(mmultx),itij(mmultx,mmult,6),isig(mmultx,mmult) |
---|
| 2532 | common/c8/ betexp(mmultx,4),fac4(mmultx,mmult),fac2(mmultx) |
---|
| 2533 | common/c9/ fact(0:mmul,0:mmul2),sign(4,4),cosav(0:mmul) |
---|
| 2534 | common/c10/ ifacd2(4,mmultf),facd2(2,mmultf),ham(2,mmultf) |
---|
| 2535 | common/c11/ det1(2,0:mmul),det(2,-mmul:mmul,0:mmul,0:mmul) |
---|
| 2536 | common/c12/ ihv(mmultx),iplane(mmultx,2) |
---|
| 2537 | common/c13/ fact0(-mmul:mmul,mh),factb(2,-mmul:mmul,mh) |
---|
| 2538 | common/c14/ ifacta(2,-mmul:mmul,mh),ifact4(4,-mmul:mmul,mh) |
---|
| 2539 | common/indeces/j70,j71,j72,j73,j74,j75,j76,j77,j78,j79 |
---|
| 2540 | common/table_sizes/table_size_70,table_size_71,table_size_72, & |
---|
| 2541 | &table_size_73,table_size_74,table_size_75, & |
---|
| 2542 | &table_size_76,table_size_77,table_size_78, & |
---|
| 2543 | &table_size_79 |
---|
| 2544 | integer int_to_write(11) |
---|
| 2545 | double precision double_to_write(11) |
---|
| 2546 | dimension imu(2) |
---|
| 2547 | integer data_size |
---|
| 2548 | !--------------------------------------------------------------------- |
---|
| 2549 | ic=imu(1) |
---|
| 2550 | id=imu(2) |
---|
| 2551 | if(icase.eq.0) then |
---|
| 2552 | |
---|
| 2553 | !--- Zero order multiple strength order |
---|
| 2554 | |
---|
| 2555 | !--- handle unit 71 (detune_1_all) |
---|
| 2556 | |
---|
| 2557 | if(iu_on.eq.2.or.iu_on.eq.3) then |
---|
| 2558 | table_name = 'detune_1_all ' |
---|
| 2559 | data_size = 2*ih1+2+j71 |
---|
| 2560 | call fit_table(table_name,data_size,table_size_71) |
---|
| 2561 | do i=ih1,0,-1 |
---|
| 2562 | write(71,10040) ic,1,det(1,ic,i,ih1-i)/pi2,i,ih1-i |
---|
| 2563 | int_to_write(1) = ic |
---|
| 2564 | int_to_write(2) = 1 |
---|
| 2565 | double_to_write(3) = det(1,ic,i,ih1-i)/pi2 |
---|
| 2566 | int_to_write(4) = i |
---|
| 2567 | int_to_write(5) = ih1-i |
---|
| 2568 | call write_table(table_name,5,int_to_write,double_to_write) |
---|
| 2569 | j71 = j71 + 1 |
---|
| 2570 | call augment_count(table_name) |
---|
| 2571 | enddo |
---|
| 2572 | do i=ih1-1,0,-1 |
---|
| 2573 | write(71,10040) ic,2, & |
---|
| 2574 | &det(1,ic,i,ih1-i)*dble(ih1-i)/dble(i+1)/pi2,i+1,ih1-i-1 |
---|
| 2575 | int_to_write(2) = 2 |
---|
| 2576 | double_to_write(3) = det(1,ic,i,ih1-i)*dble(ih1-i)/dble(i+1)& |
---|
| 2577 | &/pi2 |
---|
| 2578 | int_to_write(4) = i+1 |
---|
| 2579 | int_to_write(5) = ih1-i-1 |
---|
| 2580 | call write_table(table_name,5,int_to_write,double_to_write) |
---|
| 2581 | j71 = j71 + 1 |
---|
| 2582 | call augment_count(table_name) |
---|
| 2583 | enddo |
---|
| 2584 | write(71,10040) ic,2,det(2,ic,0,ih1)/pi2,0,ih1 |
---|
| 2585 | double_to_write(3) = det(2,ic,0,ih1)/pi2 |
---|
| 2586 | int_to_write(4) = 0 |
---|
| 2587 | int_to_write(5) = ih1 |
---|
| 2588 | call write_table(table_name,5,int_to_write,double_to_write) |
---|
| 2589 | j71 = j71 + 1 |
---|
| 2590 | call augment_count(table_name) |
---|
| 2591 | endif |
---|
| 2592 | |
---|
| 2593 | !--- First order multiple strength order |
---|
| 2594 | |
---|
| 2595 | !--- handle unit 70 (detune_1_end) |
---|
| 2596 | |
---|
| 2597 | else if(icase.eq.1) then |
---|
| 2598 | write(6,10000) |
---|
| 2599 | table_name = 'detune_1_end ' |
---|
| 2600 | data_size = 2*ih1+2+j70 |
---|
| 2601 | call fit_table(table_name,data_size,table_size_70) |
---|
| 2602 | do i=ih1,0,-1 |
---|
| 2603 | write(6,10010) det(1,ic,i,ih1-i)/pi2,i,ih1-i |
---|
| 2604 | if(iu_on.eq.1.or.iu_on.eq.3) then |
---|
| 2605 | write(70,10040) ic,1,det(1,ic,i,ih1-i)/pi2,i,ih1-i |
---|
| 2606 | int_to_write(1) = ic |
---|
| 2607 | int_to_write(2) = 1 |
---|
| 2608 | double_to_write(3) = det(1,ic,i,ih1-i)/pi2 |
---|
| 2609 | int_to_write(4) = i |
---|
| 2610 | int_to_write(5) = ih1-i |
---|
| 2611 | call write_table(table_name,5,int_to_write,double_to_write) |
---|
| 2612 | j70 = j70 + 1 |
---|
| 2613 | call augment_count(table_name) |
---|
| 2614 | endif |
---|
| 2615 | enddo |
---|
| 2616 | write(6,10020) |
---|
| 2617 | do i=ih1-1,0,-1 |
---|
| 2618 | write(6,10010) det(1,ic,i,ih1-i)*dble(ih1-i)/dble(i+1)/pi2, & |
---|
| 2619 | &i+1,ih1-i-1 |
---|
| 2620 | if(iu_on.eq.1.or.iu_on.eq.3) then |
---|
| 2621 | write(70,10040) ic,2,det(1,ic,i,ih1-i)*dble(ih1-i)/dble(i+1)& |
---|
| 2622 | &/pi2,i+1,ih1-i-1 |
---|
| 2623 | int_to_write(2) = 2 |
---|
| 2624 | double_to_write(3) = det(1,ic,i,ih1-i)*dble(ih1-i)/dble(i+1)& |
---|
| 2625 | &/pi2 |
---|
| 2626 | int_to_write(4) = i+1 |
---|
| 2627 | int_to_write(5) = ih1-i-1 |
---|
| 2628 | call write_table(table_name,5,int_to_write,double_to_write) |
---|
| 2629 | j70 = j70 + 1 |
---|
| 2630 | call augment_count(table_name) |
---|
| 2631 | endif |
---|
| 2632 | enddo |
---|
| 2633 | write(6,10010) det(2,ic,0,ih1)/pi2,0,ih1 |
---|
| 2634 | if(iu_on.eq.1.or.iu_on.eq.3) then |
---|
| 2635 | write(70,10040) ic,2,det(2,ic,0,ih1)/pi2,0,ih1 |
---|
| 2636 | double_to_write(3) = det(2,ic,0,ih1)/pi2 |
---|
| 2637 | int_to_write(4) = 0 |
---|
| 2638 | int_to_write(5) = ih1 |
---|
| 2639 | call write_table(table_name,5,int_to_write,double_to_write) |
---|
| 2640 | j70 = j70 + 1 |
---|
| 2641 | call augment_count(table_name) |
---|
| 2642 | endif |
---|
| 2643 | else if(icase.eq.2) then |
---|
| 2644 | |
---|
| 2645 | !--- Second order multiple strength order |
---|
| 2646 | |
---|
| 2647 | !--- handle unit 72 (detune_2_hor) |
---|
| 2648 | |
---|
| 2649 | if(iu_on.eq.1.or.iu_on.eq.3) then |
---|
| 2650 | write(6,10000) |
---|
| 2651 | table_name = 'detune_2_hor ' |
---|
| 2652 | data_size = ih1+1+j72 |
---|
| 2653 | call fit_table(table_name,data_size,table_size_72) |
---|
| 2654 | endif |
---|
| 2655 | do i=ih1,0,-1 |
---|
| 2656 | write(6,10010) det(1,ic,i,ih1-i)/pi2,i,ih1-i |
---|
| 2657 | if(iu_on.eq.1.or.iu_on.eq.3) then |
---|
| 2658 | write(72,10050) ic,id,det(1,ic,i,ih1-i)/pi2,i,ih1-i |
---|
| 2659 | int_to_write(1) = ic |
---|
| 2660 | int_to_write(2) = id |
---|
| 2661 | double_to_write(3) = det(1,ic,i,ih1-i)/pi2 |
---|
| 2662 | int_to_write(4) = i |
---|
| 2663 | int_to_write(5) = ih1-i |
---|
| 2664 | call write_table(table_name,5,int_to_write,double_to_write) |
---|
| 2665 | j72 = j72 + 1 |
---|
| 2666 | call augment_count(table_name) |
---|
| 2667 | endif |
---|
| 2668 | enddo |
---|
| 2669 | |
---|
| 2670 | !--- handle unit 73 (detune_2_ver) |
---|
| 2671 | |
---|
| 2672 | if(iu_on.eq.1.or.iu_on.eq.3) then |
---|
| 2673 | write(6,10020) |
---|
| 2674 | table_name = 'detune_2_ver ' |
---|
| 2675 | data_size = ih1+2+j73 |
---|
| 2676 | call fit_table(table_name,data_size,table_size_73) |
---|
| 2677 | endif |
---|
| 2678 | do i=ih1-1,0,-1 |
---|
| 2679 | write(6,10010) det(1,ic,i,ih1-i)*dble(ih1-i)/ & |
---|
| 2680 | &dble(i+1)/pi2,i+1,ih1-i-1 |
---|
| 2681 | if(iu_on.eq.1.or.iu_on.eq.3) then |
---|
| 2682 | write(73,10050) ic,id, & |
---|
| 2683 | &det(1,ic,i,ih1-i)*dble(ih1-i)/dble(i+1)/pi2,i+1,ih1-i-1 |
---|
| 2684 | int_to_write(1) = ic |
---|
| 2685 | int_to_write(2) = id |
---|
| 2686 | double_to_write(3) = det(1,ic,i,ih1-i)*dble(ih1-i)/dble(i+1)& |
---|
| 2687 | &/pi2 |
---|
| 2688 | int_to_write(4) = i+1 |
---|
| 2689 | int_to_write(5) = ih1-i-1 |
---|
| 2690 | call write_table(table_name,5,int_to_write,double_to_write) |
---|
| 2691 | j73 = j73 + 1 |
---|
| 2692 | call augment_count(table_name) |
---|
| 2693 | endif |
---|
| 2694 | enddo |
---|
| 2695 | write(6,10010) det(2,ic,0,ih1)/pi2,0,ih1 |
---|
| 2696 | if(iu_on.eq.1.or.iu_on.eq.3) then |
---|
| 2697 | write(73,10050) ic,id, & |
---|
| 2698 | &det(2,ic,0,ih1)/pi2,0,ih1 |
---|
| 2699 | double_to_write(3) = det(2,ic,0,ih1)/pi2 |
---|
| 2700 | int_to_write(4) = 0 |
---|
| 2701 | int_to_write(5) = ih1 |
---|
| 2702 | call write_table(table_name,5,int_to_write,double_to_write) |
---|
| 2703 | j73 = j73 + 1 |
---|
| 2704 | call augment_count(table_name) |
---|
| 2705 | endif |
---|
| 2706 | do i=-mmul,mmul |
---|
| 2707 | do j=0,mmul |
---|
| 2708 | do k=0,mmul |
---|
| 2709 | det(1,i,j,k)=zero |
---|
| 2710 | det(2,i,j,k)=zero |
---|
| 2711 | enddo |
---|
| 2712 | enddo |
---|
| 2713 | enddo |
---|
| 2714 | do i=0,mmul |
---|
| 2715 | det1(1,i)=zero |
---|
| 2716 | det1(2,i)=zero |
---|
| 2717 | enddo |
---|
| 2718 | endif |
---|
| 2719 | return |
---|
| 2720 | !--------------------------------------------------------------------- |
---|
| 2721 | 10000 format(/80('-')/8x,'Horizontal Coefficient',4x,'E-x',3x,'E-y'/ & |
---|
| 2722 | &80('-')) |
---|
| 2723 | 10010 format(10x,1pe20.12,2i6) |
---|
| 2724 | 10020 format(/80('-')/8x,' Vertical Coefficient',4x,'E-x',3x,'E-y'/ & |
---|
| 2725 | &80('-')) |
---|
| 2726 | 10040 format(i4,3x,i4,2x,1pe23.15,i4,1x,i4) |
---|
| 2727 | 10050 format(i4,4x,i4,3x,1pe23.15,1x,i4,1x,i4) |
---|
| 2728 | |
---|
| 2729 | end subroutine detwri |
---|
| 2730 | subroutine sortres(icase,im1,ic) |
---|
| 2731 | !--------------------------------------------------------------------- |
---|
| 2732 | !-- |
---|
| 2733 | !--Order the resonance terms |
---|
| 2734 | !-- |
---|
| 2735 | !--------------------------------------------------------------------- |
---|
| 2736 | implicit none |
---|
| 2737 | integer ic,icase,idum1,idum2,idum3,idum4,idum5,idum6,im,im1,is, & |
---|
| 2738 | &is0,k0,k1,k2,k3,k4,l2,l4,m,mh,mmul,mmul2,mmult,mmultf,mmultw, & |
---|
| 2739 | &mmultx,nblz |
---|
| 2740 | double precision dum1,dum2,dum3,dum4,pieni |
---|
| 2741 | integer icc,icd,ifacd2,ifact4,ifacta,ihv,iplane,isig,itij,ityc, & |
---|
| 2742 | &iu_on,n1,n2 |
---|
| 2743 | double precision beta,betexp,bstr,cosav,det,det1,etl,etl1,etl2, & |
---|
| 2744 | &fac2,fac4,facd2,fact,fact0,factb,four,ham,one,phi,pi,pi2,pihi,qx, & |
---|
| 2745 | &qy,sbeta,sign,two,zero |
---|
| 2746 | real tlim,time0,time1,time2,time3,time4 |
---|
| 2747 | real tim1,tim2,tim3,tim4,tim5,tim6,tim7,tim8 |
---|
| 2748 | character*16 strn |
---|
| 2749 | character*18 comment |
---|
| 2750 | parameter(pieni=1d-34,nblz=4000,mmul=11,mh=100) |
---|
| 2751 | parameter(mmul2=mmul+1,mmult=8*mmul**2) |
---|
| 2752 | parameter(mmultw=2*mmul,mmultx=mmult/7,mmultf=2*mmultx*mmult) |
---|
| 2753 | common/c0/ pi,pi2,pihi,zero,one,two,four,comment(-mmul:mmul) |
---|
| 2754 | common/c1/ tlim,time0,time1,time2,time3,time4 |
---|
| 2755 | common/c2/ tim1,tim2,tim3,tim4,tim5,tim6,tim7,tim8 |
---|
| 2756 | common/c3/ bstr(-mmul:mmul,nblz),strn(-mmul:mmul,nblz) |
---|
| 2757 | common/c4/ beta(2,-mmul:mmul,nblz),sbeta(2,-mmul:mmul,nblz) |
---|
| 2758 | common/c4q/ phi(2,-mmul:mmul,nblz) |
---|
| 2759 | common/c5/ etl(-mmul:mmul,nblz),ityc(-mmul:mmul) |
---|
| 2760 | common/c6/ etl1,etl2,qx,qy,iu_on,n1,n2 |
---|
| 2761 | common/c7/ icc,icd(mmultx),itij(mmultx,mmult,6),isig(mmultx,mmult) |
---|
| 2762 | common/c8/ betexp(mmultx,4),fac4(mmultx,mmult),fac2(mmultx) |
---|
| 2763 | common/c9/ fact(0:mmul,0:mmul2),sign(4,4),cosav(0:mmul) |
---|
| 2764 | common/c10/ ifacd2(4,mmultf),facd2(2,mmultf),ham(2,mmultf) |
---|
| 2765 | common/c11/ det1(2,0:mmul),det(2,-mmul:mmul,0:mmul,0:mmul) |
---|
| 2766 | common/c12/ ihv(mmultx),iplane(mmultx,2) |
---|
| 2767 | common/c13/ fact0(-mmul:mmul,mh),factb(2,-mmul:mmul,mh) |
---|
| 2768 | common/c14/ ifacta(2,-mmul:mmul,mh),ifact4(4,-mmul:mmul,mh) |
---|
| 2769 | !--------------------------------------------------------------------- |
---|
| 2770 | im=iabs(im1) |
---|
| 2771 | is0=0 |
---|
| 2772 | is=0 |
---|
| 2773 | do k0=im,0,-1 |
---|
| 2774 | do k1=k0,0,-1 |
---|
| 2775 | do k2=k1,0,-1 |
---|
| 2776 | l2=k1-k2 |
---|
| 2777 | do k3=k0-k1,k0-k1,-1 |
---|
| 2778 | do k4=k3,0,-1 |
---|
| 2779 | l4=k3-k4 |
---|
| 2780 | do m=1,ic |
---|
| 2781 | if(icase.eq.2.and. & |
---|
| 2782 | &k2.eq.ifact4(1,im1,m).and.l2.eq.ifact4(2,im1,m) & |
---|
| 2783 | &.and.k4.eq.ifact4(3,im1,m).and. & |
---|
| 2784 | &l4.eq.ifact4(4,im1,m)) then |
---|
| 2785 | is0=is0+1 |
---|
| 2786 | dum1=fact0(im1,is0) |
---|
| 2787 | dum2=factb(1,im1,is0) |
---|
| 2788 | dum3=factb(2,im1,is0) |
---|
| 2789 | idum1=ifacta(1,im1,is0) |
---|
| 2790 | idum2=ifacta(2,im1,is0) |
---|
| 2791 | idum3=ifact4(1,im1,is0) |
---|
| 2792 | idum4=ifact4(2,im1,is0) |
---|
| 2793 | idum5=ifact4(3,im1,is0) |
---|
| 2794 | idum6=ifact4(4,im1,is0) |
---|
| 2795 | fact0(im1,is0)=fact0(im1,m) |
---|
| 2796 | factb(1,im1,is0)=factb(1,im1,m) |
---|
| 2797 | factb(2,im1,is0)=factb(2,im1,m) |
---|
| 2798 | ifacta(1,im1,is0)=ifacta(1,im1,m) |
---|
| 2799 | ifacta(2,im1,is0)=ifacta(2,im1,m) |
---|
| 2800 | ifact4(1,im1,is0)=ifact4(1,im1,m) |
---|
| 2801 | ifact4(2,im1,is0)=ifact4(2,im1,m) |
---|
| 2802 | ifact4(3,im1,is0)=ifact4(3,im1,m) |
---|
| 2803 | ifact4(4,im1,is0)=ifact4(4,im1,m) |
---|
| 2804 | fact0(im1,m)=dum1 |
---|
| 2805 | factb(1,im1,m)=dum2 |
---|
| 2806 | factb(2,im1,m)=dum3 |
---|
| 2807 | ifacta(1,im1,m)=idum1 |
---|
| 2808 | ifacta(2,im1,m)=idum2 |
---|
| 2809 | ifact4(1,im1,m)=idum3 |
---|
| 2810 | ifact4(2,im1,m)=idum4 |
---|
| 2811 | ifact4(3,im1,m)=idum5 |
---|
| 2812 | ifact4(4,im1,m)=idum6 |
---|
| 2813 | endif |
---|
| 2814 | if(icase.eq.3.and. & |
---|
| 2815 | &k2.eq.ifacd2(1,m).and.l2.eq.ifacd2(2,m) & |
---|
| 2816 | &.and.k4.eq.ifacd2(3,m).and.l4.eq.ifacd2(4,m)) then |
---|
| 2817 | is=is+1 |
---|
| 2818 | dum1=facd2(1,is) |
---|
| 2819 | dum2=facd2(2,is) |
---|
| 2820 | dum3=ham(1,is) |
---|
| 2821 | dum4=ham(2,is) |
---|
| 2822 | idum1=ifacd2(1,is) |
---|
| 2823 | idum2=ifacd2(2,is) |
---|
| 2824 | idum3=ifacd2(3,is) |
---|
| 2825 | idum4=ifacd2(4,is) |
---|
| 2826 | facd2(1,is)=facd2(1,m) |
---|
| 2827 | facd2(2,is)=facd2(2,m) |
---|
| 2828 | ham(1,is)=ham(1,m) |
---|
| 2829 | ham(2,is)=ham(2,m) |
---|
| 2830 | ifacd2(1,is)=ifacd2(1,m) |
---|
| 2831 | ifacd2(2,is)=ifacd2(2,m) |
---|
| 2832 | ifacd2(3,is)=ifacd2(3,m) |
---|
| 2833 | ifacd2(4,is)=ifacd2(4,m) |
---|
| 2834 | facd2(1,m)=dum1 |
---|
| 2835 | facd2(2,m)=dum2 |
---|
| 2836 | ham(1,m)=dum3 |
---|
| 2837 | ham(2,m)=dum4 |
---|
| 2838 | ifacd2(1,m)=idum1 |
---|
| 2839 | ifacd2(2,m)=idum2 |
---|
| 2840 | ifacd2(3,m)=idum3 |
---|
| 2841 | ifacd2(4,m)=idum4 |
---|
| 2842 | endif |
---|
| 2843 | enddo |
---|
| 2844 | enddo |
---|
| 2845 | enddo |
---|
| 2846 | enddo |
---|
| 2847 | enddo |
---|
| 2848 | enddo |
---|
| 2849 | return |
---|
| 2850 | end subroutine sortres |
---|
| 2851 | subroutine prror(ipro,ier,num) |
---|
| 2852 | implicit none |
---|
| 2853 | integer ier,ipro,num |
---|
| 2854 | !----------------------------------------------------------------------- |
---|
| 2855 | !--Error Output |
---|
| 2856 | !----------------------------------------------------------------------- |
---|
| 2857 | |
---|
| 2858 | if(ipro.eq.0) write(6,*) ' Error occurred in Reading User Input ' |
---|
| 2859 | if(ipro.eq.1) write(6,*) ' Error occurred in Program Detune, ', & |
---|
| 2860 | &' which calculates the first and second Order Detuning.' |
---|
| 2861 | if(ipro.eq.2) write(6,*) ' Error occurred in Program Distort1, ',& |
---|
| 2862 | &' which calculates the Distortion Function in first Order.' |
---|
| 2863 | if(ipro.eq.3) write(6,*) ' Error occurred in Program Distort2, ',& |
---|
| 2864 | &' which calculates the Distortion Function in second Order.' |
---|
| 2865 | goto(10,20,30,40,50,60,70,80,90,100,110,120,130,140,150,160),ier |
---|
| 2866 | 10 write(6,*) ' File 34 empty or corrupted - Program stops' |
---|
| 2867 | stop 1 |
---|
| 2868 | 20 write(6,*) ' File 34 cannot be read - Program stops' |
---|
| 2869 | stop 2 |
---|
| 2870 | 30 write(6,*) ' Number of Elements nblz: ',num,' too small!' |
---|
| 2871 | stop 3 |
---|
| 2872 | 40 write(6,*) ' Input File fc.34 corrupted ' |
---|
| 2873 | stop 4 |
---|
| 2874 | 50 write(6,*) ' Too many Cases: in Derivation mmult:',num |
---|
| 2875 | stop 5 |
---|
| 2876 | 60 write(6,*) ' Too many Cases: ic terms, increase mmult: ',num |
---|
| 2877 | stop 6 |
---|
| 2878 | 70 write(6,*) ' Too many Cases: id terms, increase mmultx: ',num |
---|
| 2879 | stop 7 |
---|
| 2880 | 80 write(6,*) ' Number of Resonance Terms exceeded, increase ', & |
---|
| 2881 | &'mmultx: ',num |
---|
| 2882 | stop 8 |
---|
| 2883 | 90 write(6,*) ' Program Error unphysical Terms in Case: ',num |
---|
| 2884 | stop 9 |
---|
| 2885 | 100 write(6,*) ' Maximum Number of Resonances: ',num,' too large ', & |
---|
| 2886 | &'==> increase mmultf' |
---|
| 2887 | stop 10 |
---|
| 2888 | 110 write(6,*) ' Inconsistent Number: ',num,' of Resonances for ', & |
---|
| 2889 | &'the mixed Case' |
---|
| 2890 | stop 11 |
---|
| 2891 | 120 write(6,*) ' Final Ordering of Resonances is not possible as ', & |
---|
| 2892 | &'their Number: ',num,' is too large.' |
---|
| 2893 | stop 12 |
---|
| 2894 | 130 write(6,*) ' There are only 2 derivatives possible - program', & |
---|
| 2895 | &' error' |
---|
| 2896 | stop 13 |
---|
| 2897 | 140 write(6,*) ' Program Specifier iprog: ',num,' must lie between', & |
---|
| 2898 | &' 1 and 7 ' |
---|
| 2899 | stop 14 |
---|
| 2900 | 150 write(6,*) ' Analysis is restricted to order: ',num |
---|
| 2901 | stop 15 |
---|
| 2902 | 160 write(6,*) ' Unit: ',num,' could not be opened ' |
---|
| 2903 | stop 16 |
---|
| 2904 | end subroutine prror |
---|
| 2905 | |
---|
| 2906 | subroutine fit_table(table_name,data_size,tab_size) |
---|
| 2907 | implicit none |
---|
| 2908 | integer data_size,tab_size,ll |
---|
| 2909 | character table_name*(*) |
---|
| 2910 | character*16 tab_name |
---|
| 2911 | |
---|
| 2912 | tab_name = table_name |
---|
| 2913 | ll = len(tab_name) |
---|
| 2914 | tab_name(ll:ll) = ' ' |
---|
| 2915 | 10 if(data_size .gt. tab_size) then |
---|
| 2916 | call double_table(tab_name) |
---|
| 2917 | tab_size = tab_size*2 |
---|
| 2918 | goto 10 |
---|
| 2919 | endif |
---|
| 2920 | |
---|
| 2921 | return |
---|
| 2922 | end subroutine fit_table |
---|
| 2923 | |
---|
| 2924 | subroutine ertab(k,table_name,column,row) |
---|
| 2925 | implicit none |
---|
| 2926 | integer k,row |
---|
| 2927 | character table_name*(*),column*(*) |
---|
| 2928 | if (k .eq.0) return |
---|
| 2929 | if (k .eq. -1) print *,"table ",table_name," does not exist" |
---|
| 2930 | if (k .eq. -2) print *," in table ",table_name,"column ",column, & |
---|
| 2931 | &" does not exist" |
---|
| 2932 | if (k .eq. -3) print *,"in table ",table_name,"row",row, & |
---|
| 2933 | &" does not exist" |
---|
| 2934 | |
---|
| 2935 | return |
---|
| 2936 | end subroutine ertab |
---|
| 2937 | |
---|
| 2938 | subroutine write_table(table_name,table_type,int_to_write, & |
---|
| 2939 | &double_to_write) |
---|
| 2940 | implicit none |
---|
| 2941 | integer table_type,int_to_write(11) |
---|
| 2942 | integer k,table_type_index(11) |
---|
| 2943 | integer j70,j71,j72,j73,j74,j75,j76,j77,j78,j79 |
---|
| 2944 | integer tab_types_5(5),tab_types_8(8),tab_types_9(9), & |
---|
| 2945 | &tab_types_11(11) |
---|
| 2946 | double precision double_to_write(11) |
---|
| 2947 | character table_name*(*) |
---|
| 2948 | character*16 name_5(5),name_8(8),name_9(9),name_11(11) |
---|
| 2949 | common/indeces/j70,j71,j72,j73,j74,j75,j76,j77,j78,j79 |
---|
| 2950 | data table_type_index/0,0,0,0,1,0,0,2,3,0,4/ |
---|
| 2951 | data tab_types_5/1,1,2,1,1/ |
---|
| 2952 | data tab_types_8/1,2,2,2,1,1,1,1/ |
---|
| 2953 | data tab_types_9/1,1,2,2,2,1,1,1,1/ |
---|
| 2954 | data tab_types_11/1,1,1,2,2,2,2,1,1,1,1/ |
---|
| 2955 | data name_5/'multipoleorder ','plane ','detuning ','h_inv_order ',& |
---|
| 2956 | &'v_inv_order '/ |
---|
| 2957 | data name_8/'multipoleorder ','cosine ','sine ','amplitude ', & |
---|
| 2958 | &'j ','k ','l ','m '/ |
---|
| 2959 | data name_9/'multipoleorder1 ','multipoleorder2 ','cosine ', & |
---|
| 2960 | &'sine ','amplitude ','j ','k ','l ','m '/ |
---|
| 2961 | data name_11/'multipoleorder ','location ','resonance ', & |
---|
| 2962 | &'position[m] ','cosine ','sine ','amplitude ','j ','k ','l ','m '/ |
---|
| 2963 | |
---|
| 2964 | goto(50,80,90,110) table_type_index(table_type) |
---|
| 2965 | 50 do k = 1,5 |
---|
| 2966 | if(tab_types_5(k) .eq. 2) then |
---|
| 2967 | call double_to_table_curr(table_name,name_5(k),double_to_write(k)) |
---|
| 2968 | else |
---|
| 2969 | call double_to_table_curr(table_name,name_5(k), & |
---|
| 2970 | &dble(int_to_write(k))) |
---|
| 2971 | endif |
---|
| 2972 | enddo |
---|
| 2973 | go to 1000 |
---|
| 2974 | |
---|
| 2975 | 80 do k = 1,8 |
---|
| 2976 | if(tab_types_8(k) .eq. 2) then |
---|
| 2977 | call double_to_table_curr(table_name,name_8(k),double_to_write(k)) |
---|
| 2978 | else |
---|
| 2979 | call double_to_table_curr(table_name,name_8(k), & |
---|
| 2980 | &dble(int_to_write(k))) |
---|
| 2981 | endif |
---|
| 2982 | enddo |
---|
| 2983 | go to 1000 |
---|
| 2984 | |
---|
| 2985 | 90 do k = 1,9 |
---|
| 2986 | if(tab_types_9(k) .eq. 2) then |
---|
| 2987 | call double_to_table_curr(table_name,name_9(k),double_to_write(k)) |
---|
| 2988 | else |
---|
| 2989 | call double_to_table_curr(table_name,name_9(k), & |
---|
| 2990 | &dble(int_to_write(k))) |
---|
| 2991 | endif |
---|
| 2992 | enddo |
---|
| 2993 | goto 1000 |
---|
| 2994 | |
---|
| 2995 | 110 do k = 1,11 |
---|
| 2996 | if(tab_types_11(k) .eq. 2) then |
---|
| 2997 | call double_to_table_curr(table_name,name_11(k),double_to_write(k)) |
---|
| 2998 | else |
---|
| 2999 | call double_to_table_curr(table_name,name_11(k), & |
---|
| 3000 | &dble(int_to_write(k))) |
---|
| 3001 | endif |
---|
| 3002 | enddo |
---|
| 3003 | |
---|
| 3004 | 1000 return |
---|
| 3005 | end subroutine write_table |
---|