source: PSPA/madxPSPA/src/ibsdb.f90 @ 430

Last change on this file since 430 was 430, checked in by touze, 11 years ago

import madx-5.01.00
File size: 32.0 KB
Line 
1! **************************************************************************
2! Note by F. Antoniou and F. Zimmermann March 2012
3! Note that in this version the dispersion is corrected (multiplied by beta)
4! within the module as in the twiss table the disperion is given in the pt
5! frame (dx/dpt) while for the ibs calculations the dx/dp is needed.
6! ***************************************************************************
7
8subroutine enprem
9  use ibsdbfi
10  implicit none
11  !----------------------------------------------------------------------*
12  ! Purpose:                                                             *
13  !   Print emittances and sigmas.                                       *
14  !----------------------------------------------------------------------*
15  !---- Double precision version.
16  !----------------------------------------------------------------------*
17  double precision ten3p,ten6p
18  parameter(ten6p=1d6,ten3p=1d3)
19
20  !---- Emittances and sigmas.
21  write (*, 910) ten6p * ex, ten3p * sigx,                          &
22       ten6p * ey, ten3p * sigy,                                         &
23       ten6p * et, ten3p * sigt, ten3p * sige
24
25910 format(' '/' Emittances:'/' '/                                    &
26       t6,'Ex',t16,e16.6,' pi*mm*mrad',t48,'sigx',t58,f14.6,' mm'/       &
27       t6,'Ey',t16,e16.6,' pi*mm*mrad',t48,'sigy',t58,f14.6,' mm'/       &
28       t6,'Et',t16,e16.6,' pi*mm*mrad',t48,'sigt',t58,f14.6,' mm',       &
29       t88,'sigE',t96,f14.6,' 1/1000')
30
31end subroutine enprem
32
33! **********************************************************
34
35subroutine enprgl
36  use ibsdbfi
37  implicit none
38
39
40  !----------------------------------------------------------------------*
41  ! Purpose:                                                             *
42  !   Print global data for machine.                                     *
43  !----------------------------------------------------------------------*
44  !---- Double precision version.
45  !----------------------------------------------------------------------*
46  logical frad
47  integer n1
48  double precision eta,gamtr,t0,get_value,zero,one
49  parameter(zero=0d0,one=1d0)
50
51  n1 = get_value('probe ','radiate ')
52  frad = n1.ne.0
53
54  !---- Global parameters.
55
56  if (alfa .gt. zero) then
57     gamtr = sqrt(one / alfa)
58  else if (alfa .eq. zero) then
59     gamtr = zero
60  else
61     gamtr = - sqrt(-one / alfa)
62  endif
63  t0 = one / freq0
64  eta = alfa - one / gamma**2
65
66  write (*, 910) frad, circ, freq0, t0, alfa,                       &
67       eta, gamtr, currnt, bunch, parnum, en0, gamma, beta
68   
69910 format(' '/' Global parameters for the machine: ',//,             &
70       'radiate = ',l1,':'/' '/                                          &
71       t6,'C',t16,f14.6,' m',t46,'f0',t56,f14.6,' MHz',                  &
72       t86,'T0',t96,f14.6,' microseconds'/                               &
73       t6,'alfa',t16,e18.6,t46,'eta',t56,e18.6,                          &
74       t86,'gamma(tr)',t96,f14.6/                                        &
75       t6,'Bcurrent',t16,f14.6,' A/bunch',t46,'Kbunch',t56,I6,           &
76       t86,'Npart',t96,e18.6,' per bunch'/                               &
77       t6,'E',t16,f14.6,' GeV',t46,'gamma',t56,f14.6,                    &
78       t86,'beta',t96,f14.6)
79
80end subroutine enprgl
81! *************************************************************
82subroutine cavprt()
83
84  use name_lenfi
85  implicit none
86
87  integer i,lg,code,get_string,restart_sequ,advance_node
88  double precision get_value,node_value,el,rfv,rff,rfl,deltap
89  character(name_len) sequ_name,el_name
90
91  lg = get_string('sequence ', 'name ', sequ_name)
92  if (lg .gt. 0) print *, 'sequence name: ', sequ_name(:lg)
93  i = restart_sequ()
9410 continue
95  code = node_value('mad8_type ')
96  if(code.eq.39) code=15
97  if(code.eq.38) code=24
98  if (code .eq. 10) then      ! cavity
99     lg = get_string('element ', 'name ', el_name)
100     el = node_value('l ')
101     rfv = node_value('volt ')
102     rff = node_value('freq ')
103     rfl = node_value('lag ')
104     deltap = get_value('probe ','deltap ')
105     print '(a,5g14.6)', el_name(:lg), el, rfv, rff, rfl, deltap
106  endif
107  if (advance_node().ne.0)  goto 10
108end subroutine cavprt
109
110! *********************************************************************
111subroutine twclog(bxbar, bybar,dxbar,dybar, const)
112  use ibsdbfi
113  use physconsfi
114  implicit none
115
116  !----------------------------------------------------------------------*
117  ! Purpose:                                                             *
118  !   Calculation of Coulomb logarithm (and print)                       *
119  !   based on the formulae in AIP physics vade mecum p.264 (1981)       *
120  ! Input:                                                               *
121  !   BXBAR     (real)    Average horizontal beta.                       *
122  !   BYBAR     (real)    Average vertical beta.                         *
123  ! Output:                                                              *
124  !   CONST     (real)    Constant in eq. (IV.9.1), ZAP user's manual.   *
125  !----------------------------------------------------------------------*
126  !---- Double precision version.
127  !----------------------------------------------------------------------*
128  logical fbch
129  integer n
130  double precision get_value,bgam,bxbar,bybar,dxbar,dybar,cbunch,const,coulog,  &
131       debyel,densty,etrans,pnbtot,qion,rmax,rmin,rmincl,rminqm,sigtcm,  &
132       sigxcm,sigycm,tempev,vol,pi,get_variable,zero,two,four,eight,ot2, &
133       ft8,ot5,ttm3,fac1,fac2
134  parameter(zero=0d0,two=2d0,four=4d0,eight=8d0,ot2=1d2,ft8=5d8,    &
135       ot5=1d5,ttm3=2d-3,fac1=743.4d0,fac2=1.44d-7)
136
137  pi=get_variable('pi ')
138 
139  ! **************************** DB *********************
140  n = get_value('probe ', 'bunched ')
141  fbch = n.ne.0
142  ! *****************************************************
143
144  !---- Calculate transverse temperature as 2*P*X',
145  !     i.e., assume the transverse energy is temperature/2.
146  qion   = abs(charge)
147  etrans = ft8 * (gammas * en0 - amass) * (ex / bxbar)
148  tempev = two * etrans
149
150  !---- Calculate beam volume to get density (in cm**-3).
151
152  sigxcm = ot2 * sqrt(ex * bxbar + (dxbar * sige)**2)
153  sigycm = ot2 * sqrt(ey * bybar + (dybar * sige)**2)
154  sigtcm = ot2 * sigt
155  if (fbch) then
156     vol    = eight * sqrt(pi**3) * sigxcm * sigycm * sigtcm
157     densty = parnum / vol
158  else
159     vol    = four * pi * sigxcm * sigycm * ot2 * circ
160     pnbtot = currnt * circ / (qion * qelect * betas * clight)
161     densty = pnbtot / vol
162  endif
163
164  !---- Calculate RMAX as smaller of SIGXCM and DEBYE length.
165  debyel = fac1 * sqrt(tempev/densty) / qion
166  rmax   = min(sigxcm,debyel)
167
168  !---- Calculate RMIN as larger of classical distance of closest approach
169  !     or quantum mechanical diffraction limit from nuclear radius.
170  rmincl = fac2 * qion**2 / tempev
171  rminqm = hbar*clight*ot5 / (two*sqrt(ttm3*etrans*amass))
172  rmin   = max(rmincl,rminqm)
173  coulog = log(rmax/rmin)
174  bgam = betas * gammas
175  qion   = abs(charge)
176  if (fbch) then
177     const = parnum * coulog * arad**2 * clight / (eight * pi * betas&
178          **3 * gammas**4 * ex * ey * sige * sigt)
179     cbunch = qion * parnum * qelect * betas * clight / circ
180  else
181     const = currnt * coulog * arad**2 /                             &
182          (four * sqrt(pi) * qion * qelect * bgam**4 * ex * ey * sige)
183  endif
184  write (*, 910) const
185
186  write (*, 920) en0, betas, gammas, coulog
187
188  !---- Print warning here if Coulomb logarithm gave bad results.
189  !     Usually this error is due to a starting guess far from
190  !     the equilibrium value.
191  if (coulog .lt. zero) then
192     call aawarn('TWCLOG', 'Coulomb logarithm gives invalid'         &
193          // ' result --- check input parameters.')
194  endif
195
196  write (*, 940) ex, ey
197
198  if (fbch) then
199     write (*, 950) sige, sigt, parnum, cbunch
200  else
201     write (*, 960) sige, currnt
202  endif
203
204910 format(' '/5x,'CONST               = ',1p,e14.6)
205920 format(' '/5x,'ENERGY              = ',f14.6,' GeV'/              &
206       5x,'BETA                = ',f14.6/                                &
207       5x,'GAMMA               = ',f14.3/                                &
208       5x,'COULOMB LOG         = ',f14.3)
209940 format(' '/5x,'X-emittance         = ',1p,e14.6,' m*rad'/         &
210       5x,'Y-emittance         = ',   e14.6,' m*rad')
211950 format(' '/5x,'Momentum spread     = ',1p,e14.6/                  &
212       5x,'Bunch length        = ',0p,f14.6,' m'/' '/                    &
213       5x,'Particles per bunch = ',1p,e14.6/                             &
214       5x,'Bunch current       = ',1p,e14.6,' A')
215960 format(' '/5x,'Momentum spread     = ',1p,e14.6/' '/              &
216       5x,'Current             = ',0p,f14.6,' A'/' ')
217
218end subroutine twclog
219! *********************************************************************
220subroutine ibs
221
222  use ibsdbfi
223  use physconsfi
224  use name_lenfi
225  implicit none
226
227
228  !----------------------------------------------------------------------*
229  ! Purpose:                                                             *
230  !   INTRABEAM SCATTERING, IBS Command                                  *
231  !   These routines are a much reduced version of IBS as taken          *
232  !   from the program ZAP, written by M. Zisman.                        *
233  !   One should refer to the ZAP USERS MANUAL LBL-21270 UC-28.          *
234  ! Attribute:                                                           *
235  !   TABLE     (name)    Name of Twiss table.                           *
236  !----------------------------------------------------------------------*
237  integer step,i,j,flag,testtype,range(2),n,get_option,double_from_table_row,    &
238       restart_sequ,advance_to_pos
239  double precision tol,alx,alxbar,alxwtd,aly,alybar,ax1,ax2,ay1,ay2,&
240       betax,betay,beteff,bx1,bx2,bxbar,bxinv,by1,by2,bybar,byinv,bywtd, &
241       const,dels,dpx,dpx1,dpx2,dpxbr,dpxwtd,dx,dx1,dx2,dxbar,dxwtd,     &
242       hscrpt,hscwtd,s1,s2,ss2,l1,l2,ll2,salxb,salyb,sbxb,sbxinv,sbyb,sbyinv,sdpxb,    &
243       sdxb,taul,taux,tauy,tavl,tavlc,tavx,tavxc,tavy,tavyc,tlbar,tlidc, &
244       tlwtd,txbar,txidc,txwtd,tybar,tyidc,tywtd,wnorm,sdum,get_value,   &
245       get_variable,zero,one,two,half,dy,dy1,dy2,dybar,dywtd,hscrpty,    &
246       hscwtdy,sdpyb,sdyb,dpy,dpy1,dpy2,dpybr,dpywtd,beteffy,alywtd
247  parameter(zero=0d0,one=1d0,two=2d0,half=0.5d0)
248 
249  !---- Universal physical constants.
250
251  !     Permeability of vacuum [V*s/A*m]:
252  amu0 = get_variable('amu0 ')
253  !     Permittivity of vaccum [A*S/V*m]:
254  eps0 = get_variable('eps0 ')
255  !     Reduced Plack's constant [GeV*s]:
256  hbar = get_variable('hbar ')
257
258  !---- Electromagnetic constants.
259  !     Elementary charge [A*s]:
260  qelect = get_variable('qelect ')
261
262  !---- Electron.
263  !     Rest mass [GeV]:
264  emass = get_variable('emass ')
265  !     Classical radius [m]:
266  erad = get_variable('erad ')
267  !     Reduced Compton wavelength [m]:
268  elamda = get_variable('elamda ')
269
270  !---- Proton.
271  !     Rest mass [GeV]:
272  pmass = get_variable('pmass ')
273  !     Classical radius [m]:
274  prad = get_variable('prad ')
275  !     Reduced Compton wavelength [m]:
276  plamda = get_variable('plamda ')
277
278  !---- Muon.
279  !     Rest mass [GeV]:
280  mumass = get_variable('mumass ')
281
282  ! ************* Get the parameters for the common blocks *************
283  ! *************         /machin/ and /beamdb/            *************
284
285
286  charge   = get_value('probe ', 'charge ')
287  gammas   = get_value('probe ', 'gamma ')
288  gamma    = get_value('probe ', 'gamma ')
289  en0      = get_value('probe ', 'energy ')
290  amass    = get_value('probe ', 'mass ')
291  ex       = get_value('probe ', 'ex ')
292  ey       = get_value('probe ', 'ey ')
293  et       = get_value('probe ', 'et ')
294  sigt     = get_value('probe ', 'sigt ')
295  sige     = get_value('probe ', 'sige ')
296  parnum   = get_value('probe ', 'npart ')
297  circ     = get_value('probe ', 'circ ')
298  currnt   = get_value('probe ', 'bcurrent ')
299  betas    = get_value('probe ', 'beta ')
300  beta     = get_value('probe ', 'beta ')
301  clight   = get_variable('clight ')
302  arad     = get_value('probe ', 'arad ')
303  alfa     = get_value('probe ', 'alfa ')
304  freq0    = get_value('probe ', 'freq0 ')
305  bunch   = get_value('probe ', 'kbunch ')
306 
307  ! NOTE:
308  !****************************************************************
309  ! Sige is the dE/E. dp/p needed as input for the IBS calculations
310  ! dp/p= (dE/E)/beta**2
311  !*****************************************************************
312  sige    = sige/beta/beta
313  print *, 'sige ', sige
314 
315  !  ****************** Test print ********
316  !     print *, 'Charge ', charge
317  !     print *, 'gammas ', gammas
318  !     print *, 'gamma ', gamma
319  !     print *, 'Energy ', en0
320  !     print *, 'Mass ', amass
321  !     print *, 'Ex ', ex
322  !     print *, 'Ey ', ey
323  !     print *, 'Et ', et
324  !     print *, 'sigt ', sigt
325  !     print *, 'sige ', sige
326  !     print *, 'parnum ', parnum
327  !     print *, 'circ ', circ
328  !     print *, 'currnt ', currnt
329  !     print *, 'betas ', betas
330  !     print *, 'beat ', beta
331  !     print *, 'clight ', clight
332  !     print *, 'arad ', arad
333  !     print *, 'alfa ', alfa
334  !     print *, 'freq0 ', freq0
335  !     print *, 'kbunch ', bunch
336  ! ***************************************
337
338  !---- Initialize variables to accumulate weighted average lifetimes.
339  tavlc  = zero
340  tavxc  = zero
341  tavyc  = zero
342  dxwtd  = zero
343  dpxwtd = zero
344  dywtd  = zero
345  dpywtd = zero
346  bywtd  = zero
347  alxwtd = zero
348  alywtd = zero
349  hscwtd = zero
350  hscwtdy= zero
351  wnorm  = zero
352  sbxb   = zero
353  sbyb   = zero
354  salxb  = zero
355  salyb  = zero
356  sdxb   = zero
357  sdpxb  = zero
358  sdyb   = zero
359  sdpyb  = zero
360  sbxinv = zero
361  sbyinv = zero
362
363  ! ****** Start new Twiss Table reading *****************
364  !
365  step = get_value('ibs ', 'steps ')
366  tol = get_value('ibs ', 'tolerance ')
367  !     print *, 'steps: ', step, '  tolerance: ', tol
368  call table_range('twiss ', '#s/#e ', range)
369  !     print *, 'Range for Table ', range
370  flag = double_from_table_row('twiss ', 's ', range(1), s1)
371  if (flag .ne. 0)  goto 102
372  flag = double_from_table_row('twiss ', 'l ', range(1), l1)
373  if (flag .ne. 0)  goto 102
374  flag = double_from_table_row('twiss ', 'betx ', range(1), bx1)
375  if (flag .ne. 0)  goto 102
376  flag = double_from_table_row('twiss ', 'bety ', range(1), by1)
377  if (flag .ne. 0)  goto 102
378  flag = double_from_table_row('twiss ', 'alfx ', range(1), ax1)
379  if (flag .ne. 0)  goto 102
380  flag = double_from_table_row('twiss ', 'alfy ', range(1), ay1)
381  if (flag .ne. 0)  goto 102
382  flag = double_from_table_row('twiss ', 'dx ', range(1), dx1)
383  if (flag .ne. 0)  goto 102
384  flag = double_from_table_row('twiss ', 'dpx ', range(1), dpx1)
385  if (flag .ne. 0)  goto 102
386  flag = double_from_table_row('twiss ', 'dy ', range(1), dy1)
387  if (flag .ne. 0)  goto 102
388  flag = double_from_table_row('twiss ', 'dpy ', range(1), dpy1)
389  if (flag .ne. 0)  goto 102
390 
391  j = restart_sequ()
392  do i=range(1)+1, range(2)     
393     j = advance_to_pos('twiss ', i)
394     flag = double_from_table_row('twiss ', 's ', i, ss2)
395     if (flag .ne. 0)  goto 102
396     flag = double_from_table_row('twiss ', 'l ', i, ll2)
397     if (flag .gt. 0)  goto 102
398     if (ll2 .gt. 0.0001) goto 103
399  enddo
400  103 continue
401 
402  ! NOTE by F.A & F.Z
403  ! ************************************************************************************
404  ! Added 16.01.2012 to check if the twiss is taken at the center (testtype=2) or the
405  ! exit (testtype=1) of the elements.
406  ! I testtype=1 linear interpolation will be used to
407  ! calculate the twiss at the center of the elements.
408  !*************************************************************************************
409 
410  if ((ss2-s1) .eq. ll2) then
411        testtype = 1
412        print *, 'Twiss was calculated at the exit of the elements. Twiss functions at the center &
413                      & of the elements are calculated through linear interpolation'
414  else if ((ss2-s1) .eq. (l1+ll2)/2) then
415        testtype = 2
416        print *, 'Twiss was calculated at the center of the elements. No interpolation is used'
417  endif
418
419  ! ************** Check if "ibs_table" required  ****************
420
421  n = get_option('ibs_table ')
422
423  ! ********** Start Do loop ***************
424  !
425  j = restart_sequ()
426  do i = range(1)+1, range(2)
427     j = advance_to_pos('twiss ', i)
428     flag = double_from_table_row('twiss ', 's ', i, s2)
429     if (flag .ne. 0)  goto 102
430     flag = double_from_table_row('twiss ', 'l ', i, l2)
431     if (flag .ne. 0)  goto 102
432     flag = double_from_table_row('twiss ', 'betx ', i, bx2)
433     if (flag .ne. 0)  goto 102
434     flag = double_from_table_row('twiss ', 'bety ', i, by2)
435     if (flag .ne. 0)  goto 102
436     flag = double_from_table_row('twiss ', 'alfx ', i, ax2)
437     if (flag .ne. 0)  goto 102
438     flag = double_from_table_row('twiss ', 'alfy ', i, ay2)
439     if (flag .ne. 0)  goto 102
440     flag = double_from_table_row('twiss ', 'dx ', i, dx2)
441     if (flag .ne. 0)  goto 102
442     flag = double_from_table_row('twiss ', 'dpx ', i, dpx2)
443     if (flag .ne. 0)  goto 102
444     flag = double_from_table_row('twiss ', 'dy ', i, dy2)
445     if (flag .ne. 0)  goto 102
446     flag = double_from_table_row('twiss ', 'dpy ', i, dpy2)
447     if (flag .ne. 0)  goto 102
448
449   
450    !  NOTE by F.A & F.Z
451    ! ************************************************************************************
452    ! Dispersion and Dispersion prime is multiplied by beta, in order to be in the deltap
453    ! and not the pt frame. This correction is necessary for non-relativistic beams
454    !*************************************************************************************
455       
456     if (testtype .eq. 1) then
457        dels = s2-s1
458        sdum = half * (s2 + s1)
459        betax  = half * (bx2 + bx1)
460        betay  = half * (by2 + by1)
461        alx    = half * (ax2 + ax1)
462        aly    = half * (ay2 + ay1)
463        dx     = beta * half * (dx2 + dx1)
464        dpx    = beta * half * (dpx2 + dpx1)
465        dy     = beta * half * (dy2 + dy1)
466        dpy    = beta * half * (dpy2 + dpy1)
467
468     else if (testtype .eq. 2) then
469        dels = l2
470        sdum = s2
471        betax  = bx2
472        betay  = by2
473        alx    = ax2
474        aly    = ay2
475        dx     = beta * dx2
476        dpx    = beta * dpx2
477        dy     = beta * dy2
478        dpy    = beta * dpy2
479    endif
480
481        sbxb   = sbxb + betax * dels
482        sbxinv = sbxinv + dels / betax
483        sbyb   = sbyb + betay * dels
484        sbyinv = sbyinv + dels / betay
485        salxb  = salxb + alx * dels
486        salyb  = salyb + aly * dels
487        sdxb   = sdxb + dx * dels
488        sdpxb  = sdpxb + dpx * dels
489        sdyb   = sdyb + dy * dels
490        sdpyb  = sdpyb + dpy * dels
491
492
493 
494     !*---- Calculate weighted average in region of non-zero DX's.
495     !     These values are used to calculate "average" ring lifetimes
496     !     in TWSINT.
497     if (dx .gt. zero) then
498        wnorm  = wnorm + dels
499        dxwtd  = dxwtd + dels * dx
500        dpxwtd = dpxwtd + dels * dpx
501        dywtd  = dywtd + dels * dy
502        dpywtd = dpywtd + dels * dpy
503        bywtd  = bywtd + dels / sqrt(betay)
504        alxwtd = alxwtd + dels * alx
505        alywtd = alywtd + dels * aly
506        hscrpt = betax * dpx**2 + two * alx * dx * dpx +              &
507             (one + alx**2) * dx**2 / betax
508        hscrpty = betay * dpy**2 + two * aly * dy * dpy +             &
509             (one + aly**2) * dy**2 / betay
510        hscwtd = hscwtd + dels * sqrt(hscrpt)
511        hscwtdy = hscwtdy + dels * sqrt(hscrpty)
512     endif
513
514     !---- TWSINT calculates the Bjorken/Mtingwa integral.
515     call twsint(betax, betay, alx, aly, dx, dpx, dy, dpy,           &
516          txidc, tyidc, tlidc)
517
518     !---- Accumulate contributions.
519     tavlc = tavlc + tlidc * dels
520     tavxc = tavxc + txidc * dels
521     tavyc = tavyc + tyidc * dels
522
523     ! *************** Fill "ibs_table" if required *********************
524
525     if(n.ne.0) then
526        call string_to_table_curr('ibs ', 'name ', 'name ')
527        call double_to_table_curr('ibs ','s ', sdum)
528        call double_to_table_curr('ibs ','dels ', dels)
529        call double_to_table_curr('ibs ','tli ', tlidc)
530        call double_to_table_curr('ibs ','txi ', txidc)
531        call double_to_table_curr('ibs ','tyi ', tyidc)       
532        call double_to_table_curr('ibs ','betx ', betax)
533        call double_to_table_curr('ibs ','alfx ', alx)
534        call double_to_table_curr('ibs ','dx ', dx)
535        call double_to_table_curr('ibs ','dpx ', dpx)
536        call double_to_table_curr('ibs ','bety ', betay)
537        call double_to_table_curr('ibs ','alfy ', aly)
538        call double_to_table_curr('ibs ','dy ', dy)
539        call double_to_table_curr('ibs ','dpy ', dpy)
540        call augment_count('ibs ')
541     endif
542
543     ! *********** Make sure the following lines are not moved ******
544     ! *********** by the compiler ***************
545     s1   = s2
546     bx1  = bx2
547     by1  = by2
548     ax1  = ax2
549     ay1  = ay2
550     dx1  = dx2
551     dpx1 = dpx2
552     dy1  = dy2
553     dpy1 = dpy2
554
555  enddo
556  goto 101
557102 continue
558  call aawarn('IBS ', 'table value not found, rest skipped ')
559  stop
560101 continue
561
562
563  !---- We have finished reading the lattice from MAD
564  bxbar  = sbxb / s2
565  bybar  = sbyb / s2
566  alxbar = salxb / s2
567  alybar = salyb / s2
568  dxbar  = sdxb / s2
569  dpxbr  = sdpxb / s2
570  dybar  = sdyb / s2
571  dpybr  = sdpyb / s2
572  bxinv  = sbxinv / s2
573  byinv  = sbyinv / s2
574
575  dxwtd  = dxwtd / wnorm
576  dpxwtd = dpxwtd / wnorm
577  dywtd  = dywtd / wnorm
578  dpywtd = dpywtd / wnorm
579  bywtd  = bywtd / wnorm
580  bywtd  = one / bywtd**2
581  alxwtd = alxwtd / wnorm
582  alywtd = alywtd / wnorm
583  hscwtd = (hscwtd/wnorm)**2
584  beteff = dxwtd**2 / hscwtd
585  if (hscwtdy.ne.0.d0) then
586     beteffy = dywtd**2 / hscwtdy
587  else
588     beteffy = bywtd
589  endif
590  ! ********** Compute beam sizes with average betas ************
591
592  sigx = sqrt(ex * bxbar + (dx * sige)**2)
593  sigy = sqrt(ey * bybar + (dy * sige)**2)
594
595  call enprgl
596  call enprem
597  call cavprt()
598  ! ************************************************************
599
600  !---- Integral for averaged quantities.
601  call twsint(bxbar,bybar,alxbar,alybar, dxbar,dpxbr,               &
602       dybar,dpybr,txbar,tybar,tlbar)
603
604  !---- Integral for effective quantities.
605  call twsint(beteff,beteffy,alxwtd,alywtd,dxwtd,dpxwtd,            &
606       dywtd,dpywtd,txwtd,tywtd,tlwtd)
607
608  !---- Calculate the Coulomb logarithm.
609  call twclog(bxbar, bybar, dxbar, dybar, const)
610
611  !---- Output (weighted) average values.
612  write (*, 940) bxbar, bybar, dxbar, dybar, alxbar, alybar,        &
613       dpxbr, dpybr,                                                     &
614       bxinv, byinv
615
616  !---- Output averaged values.
617  tavl   = tavlc * const / s2
618  tavx   = tavxc * const / s2
619  tavy   = tavyc * const / s2
620
621  taul   = one / tavl
622  taux   = one / tavx
623  tauy   = one / tavy
624 
625  call set_variable('ibs.tx ',taux)
626  call set_variable('ibs.ty ',tauy)
627  call set_variable('ibs.tl ',taul)
628 
629  write (*, 950) tavl, tavx, tavy, taul, taux, tauy
630
631910 format(' '/' Particle beam: ',a,10x,a,'bunched.')
632920 format(' '/' Individual lattice point lifetimes'/' '/             &
633       26x,'TLI/const',10x,'TXI/const',10x,'TYI/const'/                  &
634       27x,'(1/sec)',12x,'(1/sec)',12x,'(1/sec)'/' ')
635930 format(1x,i8,2x,a8,3x,3(1pe15.6,3x))
636940 format(' '/' Ring average values (m)'/' '/ 5x,'betx   = ',        &
637       1pe13.5,4x, 'bety   = ',1pe13.5,4x,'Dx  = ',1pe12.5,              &
638       4x,'Dy  = ',1pe12.5/                                              &
639       5x,'alfx   = ',1pe13.5,4x,'alfy   = ',1pe13.5,4x,'Dpx = ',        &
640       1pe12.5/5x,'Dpy = ',                                              &
641       1pe12.5/5x,'1/betx = ',1pe13.5,4x,'1/bety = ',1pe13.5)
642950 format(' '/5x,'(Weighted) average rates (1/sec):'/                &
643       5x,'Longitudinal= ',1p,e15.6/                                     &
644       5x,'Horizontal  = ',   e15.6/                                     &
645       5x,'Vertical    = ',   e15.6/                                     &
646       ' '/5x,'(Weighted) average lifetimes (sec):'/                     &
647       5x,'Longitudinal= ',1p,e15.6/                                     &
648       5x,'Horizontal  = ',   e15.6/                                     &
649       5x,'Vertical    = ',   e15.6/' ')
650       
651end subroutine ibs
652! *********************************************************************
653subroutine twsint(betax, betay, alx, aly, dx, dpx, dy, dpy,       &
654     txi, tyi, tli)
655
656  use ibsdbfi
657  use physconsfi
658  implicit none
659
660  !----------------------------------------------------------------------*
661  ! Purpose:                                                             *
662  !   Subroutine uses Simpson's rule integration                         *
663  !   to calculate Bjorken/Mtingwa integrals (eqn. 3.4)                  *
664  !   Particle Accelerators 13, 115 (1983)                               *
665  !                                                                      *
666  !   The expressions found in Conte/Martini                             *
667  !   Particle Accelerators 17, 1 (1985) contain two false               *
668  !   terms in the expression for tau_x which have been corrected        *
669  !   in this version of MADX;                                           *
670  !   contributions from vertical dispersion were also added;            *
671  !   AB Note by Frank Zimmermann be published (2005)                    *
672  !                                                                      *
673  !                                                                      *
674  !   Integrals are broken into decades to optimize speed.               *
675  !                                                                      *
676  !   For the VAX, values may not exceed 10**33, therefore TSTLOG=33     *
677  !   For the IBM, values may not exceed 10**74, therefore TSTLOG=74     *
678  !   (PMG, March 1988)                                                  *
679  !                                                                      *
680  !   The integral is split into MAXDEC decades with NS steps /decade.   *
681  !   TEST is used for testing convergence of the integral               *
682  ! Input:                                                               *
683  !   BETAX     (real)    Horizontal beta.                               *
684  !   BETAY     (real)    Vertical beta.                                 *
685  !   ALX       (real)    Horizontal alpha.                              *
686  !   ALY       (real)    Vertical alpha.                                *
687  !   DX        (real)    Horizontal dispersion.                         *
688  !   DPX       (real)    Derivative of horizontal dispersion.           *
689  !   DY        (real)    Vertical dispersion.                           *
690  !   DPY       (real)    Derivative of vertical dispersion.             *
691  ! Output:                                                              *
692  !   TXI       (real)    Horizontal rate / const.                       *
693  !   TYI       (real)    Vertical rate / const.                         *
694  !   TLI       (real)    Longitudinal rate / const.                     *
695  !----------------------------------------------------------------------*
696  integer iiz,iloop,maxdec,ns
697  parameter(maxdec=30,ns=50)
698  double precision a,al(31),alam,aloop,alx,am,b,betax,betay,bl(30), &
699       c1,c2,c3,ccy,chklog,cl,coeff(2),cof,cprime,cscale,cx,cy,dpx,dx,f, &
700       func,h,phi,polyl,polyx,polyy,r1,suml,sumx,sumy,td1,td2,term,tl1,  &
701       tl2,tli,tmpl,tmpx,tmpy,tx1,tx2,txi,ty1,ty2,tyi,zintl,zintx,zinty, &
702       zero,one,two,three,six,tstlog,power,ten,test,dy,dpy,aly,phiy,c1y, &
703       c2y,chy,four,onetominus20
704  parameter(zero=0d0,one=1d0,two=2d0,three=3d0,six=6d0,tstlog=74d0, &
705       power=-two/three,ten=1d1,test=1d-7,four=4d0,onetominus20=1d-20)
706  data coeff / 2d0, 4d0 /
707
708  phi    = dpx + (alx * dx / betax)
709  phiy   = dpy + (aly * dy / betay)
710  am     = one
711  c1     = (gammas * dx)**2 / (ex * betax)
712  c1y    = (gammas * dy)**2 / (ey * betay)
713  c3     = betax / ex
714  c2     = c3 * (gammas*phi)**2
715  cx     = c1 + c2
716  cl     = am * (gammas/sige)**2
717  cy     = betay / ey
718  c2y    = cy * (gammas*phiy)**2
719  chy    = c1y + c2y
720  r1     = three / cy
721  a      = cx + cl + chy + c3 + cy
722!--- corrected b 23.02.2011
723  b      = (c3 + cy) * (c1 + cl+c1y) + cy * c2 + c3 * c2y+ c3 * cy
724
725  !---- Define CPRIME=C*CSCALE to try to keep the value.
726  !     small enough for the VAX in single precision or
727  !     IBM in double precision.
728  !     Test LOG(C) to see if it needs scaling
729  cscale = one
730  chklog = log10(c3) + log10(cy) + log10(c1 + cl)
731  if (chklog .gt. tstlog) cscale = ten**(tstlog-chklog)
732  cprime = c3 * cy * cscale * (c1 + cl + c1y)
733
734  !---- Split integral into decades, with NS steps per decade.
735  !     variables to save integral segments
736  zintl  = zero
737  zintx  = zero
738  zinty  = zero
739
740  !---- Constants for integration loop.
741  !     To keep the numbers reasonable, the numerator is
742  !     scaled by 1/CPRIME and the denominator by 1/CPRIME**2.
743  !     The extra factor of CPRIME is accounted for after integrating
744  ccy    = cprime**power
745  td1    = (a - cy) * ccy
746  td2    = one / (sqrt(ccy) * cscale * cy)
747  tl1    = (two * a - three *cy - three * c3) / cprime
748  tl2    = (b - three * c3 * cy ) / cprime
749!--- corrected ty1 23.02.2011
750  ty1    = (- a + three * cy -chy - chy/cy*(c3 -                 &
751       two*gammas**2/sige**2) + two * chy * (cx +chy)/cy + six * c2y)    &
752       / cprime
753!--- corrected ty2 23.02.2011
754  ty2    = (b - c1y * (c3+cy) +chy*(cy+chy)+chy*ey*(one/ey+&
755            betax/(betay*ex))                                            &
756       *gammas**2/sige**2-chy*betax/ex*four+(one+(betax*ey)/             &
757       (betay*ex))*                                                      &
758       cx*chy+(chy**2)*(betax*ey)/(betay*ex)-chy*ey*c2*c3/betay          &
759       -c2y*(cy+c3+chy)+three*c3*(two*c2y+c1y)) / cprime - r1 / cscale
760
761!--- corrected tx1 23.02.2011
762  tx1    = (two * (a-c3-cy) * (cx - c3) - cy * cx +                         &
763       c3 * (c1y + six * c2 + c2y + two * c3 + cl - cy)) / cprime
764!--- corrected tx2 23.02.2011
765  tx2    = (c3 + cx) * ((b - c1y * (c3 + cy)) / cprime)-        &
766       six / cscale + three * c3 * cy * (cl / cprime)                    &
767       + ( six*c3*cy*c1y                               &
768       + (betay/ey+betax/ex)*chy*cx +                                    &
769       chy*(c3**2-two*cy*c3)-c2y*cx*(cy+c3)+(two*cy*c3-c3*c3)*           &
770       c2y ) / cprime
771
772  al(1)  = zero
773
774  do iloop = 1, maxdec
775     bl(iloop) = ten**iloop
776     al(iloop+1) = bl(iloop)
777     h = (bl(iloop) - al(iloop)) / ns
778     aloop = al(iloop)
779
780     !---- Evaluate Simpson's rule summation for one interval.
781     !     The integrand is calculated in the loop itself
782     if (abs(cy+aloop).gt.onetominus20) then
783        term = sqrt((cy+aloop)*ccy)*sqrt(                             &
784             (aloop*ccy*aloop+td1*aloop+td2)+aloop*c2y*(c3-cy)*ccy/(cy+aloop))
785     else
786        term = sqrt((cy+aloop)*ccy)*sqrt(                             &
787             (aloop*ccy*aloop+td1*aloop+td2))
788     endif
789     func = sqrt(aloop) / term**3
790     polyl = tl1 * aloop + tl2
791     polyx = tx1 * aloop + tx2
792     polyy = ty1 * aloop + ty2
793     suml = func * polyl
794     sumx = func * polyx
795     sumy = func * polyy
796
797
798     do iiz = 1, ns
799        alam = aloop + iiz * h
800        cof = coeff(mod(iiz,2)+1)
801        if (abs(cy+alam).gt.onetominus20) then
802           term = sqrt((cy+alam)*ccy)*sqrt(                            &
803                (alam*ccy*alam+td1*alam+td2)+alam*c2y*(c3-cy)*ccy/(cy+alam))
804        else
805           term = sqrt((cy+alam)*ccy)*sqrt(                            &
806                (alam*ccy*alam+td1*alam+td2))
807        endif
808        f = sqrt(alam) / term**3
809        polyl = tl1 * alam + tl2
810        polyx = tx1 * alam + tx2
811        polyy = ty1 * alam + ty2
812
813        suml = suml + cof * f * polyl
814        sumx = sumx + cof * f * polyx
815        sumy = sumy + cof * f * polyy
816     enddo
817
818     suml = suml - f * polyl
819     sumx = sumx - f * polyx
820     sumy = sumy - f * polyy
821     tmpl = (suml / three) * h
822     tmpx = (sumx / three) * h
823     tmpy = (sumy / three) * h
824     zintl = zintl + tmpl
825     zintx = zintx + tmpx
826     zinty = zinty + tmpy
827
828     !---- Test to see if integral has converged.
829     if (abs(tmpl/zintl) .lt. test .and.                             &
830          abs(tmpx/zintx) .lt. test .and.                                   &
831          abs(tmpy/zinty) .lt. test) goto 100
832  enddo
833  write (*, *) tmpl,zintl, tmpx,zintx, tmpy,zinty, test
834  write (*, 910) maxdec
835910 format('Bjorken/Mtingwa integrals did not converge in ',          &
836       i3,' decades.')
837  call aawarn('TWSINT: ', 'Problem with TWSINT, program stopped ')
838  stop
839100 continue
840
841  !---- Divide answers by cprime to account for scaling.
842  txi    =      (zintx / cprime)
843  tli    = cl * (zintl / cprime)
844  tyi    = cy * (zinty / cprime) 
845 
846end subroutine twsint
847! ***************************************************************
Note: See TracBrowser for help on using the repository browser.