source: PSPA/madxPSPA/src/touschek.f90 @ 476

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

import madx-5.01.00
File size: 30.1 KB
Line 
1! *********************************************************************
2subroutine cavtouschek (um,uloss,iflag)
3
4  use name_lenfi
5  use touschekfi
6  implicit none
7
8  integer i,lg,code,get_string,restart_sequ,advance_node,           &
9       double_from_table_row,flag,iflag
10  double precision get_value,node_value,el,rfv,rff,rfl,             &
11       um,harmonl,get_variable,pi,                                       &
12       phirf, c0, vrf, pc, omega, orbit5,                                &
13       twopi, ten3m,                                                     &
14       ten6p, zero, one, two, half, eta,                                 &
15       qover, fq, uloss, vrfsum, harmonlm,                               &
16       umt,synch_2
17
18
19  character(name_len) sequ_name,el_name
20  parameter(zero=0d0,one=1d0,two=2d0,half=5d-1,ten6p=1d6,           &
21       ten3m=1d-3)
22
23  !---- Initialize.
24  qover = zero
25  vrfsum = zero
26  harmonlm = ten6p
27  umt = 0.d0
28  iflag = 0
29
30  flag = double_from_table_row('summ ','synch_2 ',1,synch_2)
31  if (synch_2 .eq. 0) then
32     iflag = 1
33     uloss = 0d0
34  else
35     uloss = 2d0/3d0*arad*en0**4*beta**3*synch_2*1d3/(amass)**3
36  endif
37  twopi=get_variable('twopi ')
38  pi=get_variable('pi ')
39  um=0.d0
40
41  lg = get_string('sequence ', 'name ', sequ_name)
42
43  i = restart_sequ()
4410 continue
45  code = node_value('mad8_type ')
46  if(code.eq.39) code=15
47  if(code.eq.38) code=24
48  ! cavity
49  if (code .eq. 10) then
50     lg = get_string('element ', 'name ', el_name)
51     el = node_value('l ')
52     rfv = node_value('volt ')
53     rff = node_value('freq ')
54     if (rff.eq.zero.or.rfv.eq.zero) goto 11
55     rfl = node_value('lag ')
56
57     harmonl = 1.E+06*rff*circ/clight
58
59     pc = get_value('probe ','pc ')
60     !        print *, ' pc ', pc
61     omega = rff * ten6p * twopi / clight
62     vrf   = rfv * ten3m / (pc * (one + deltap))
63
64     flag = double_from_table_row('twiss ','t ',1,orbit5)
65     phirf = rfl * twopi - omega * orbit5
66     c0 =   vrf*charge
67     if (cos(phirf).lt.0) vrf=-vrf
68     eta = alfa - one / gammas**2
69     if (uloss.ne.zero) then
70        qover = qover+rfv/uloss*charge
71        vrfsum = vrfsum+rfv/harmonl*charge
72        harmonlm=min(harmonl,harmonlm)
73     else
74        umt = umt+2d0/(harmonl*eta*pi)*c0
75     endif
76
77  endif
7811 continue
79
80  if (advance_node().ne.0)  goto 10
81
82  if (uloss.ne.zero) then
83     fq = two*(sqrt(one-one/qover**2)*vrfsum*harmonlm                &
84          -uloss*acos(one/qover))
85     um = ten3m/(harmonlm*eta*pi)*fq/(pc*(one+deltap))
86  else
87     um = umt
88  endif
89
90  um=abs(um)*beta**2
91
92  RETURN
93END subroutine cavtouschek
94
95! *********************************************************************
96subroutine touschek
97
98  use name_lenfi
99  use physconsfi
100  use touschekfi
101  implicit none
102
103  !----------------------------------------------------------------------*
104  ! Purpose:                                                             *
105  !   TOUSCHEK SCATTERING, TOUSCHEK Command                              *
106  !   These routines implement the formalism from Piwinski               *
107  !   (DESY 98-179 & A. Chao/M. Tigner, Handbook of Acc. Physics)        *
108  ! Attribute:                                                           *
109  !   TABLE     (name)    Name of Twiss table.                           *
110  !----------------------------------------------------------------------*
111  integer i,j,flag,range(2),n,get_option,double_from_table_row,lp,      &
112       restart_sequ,string_from_table_row,advance_to_pos,get_string,iflag
113  double precision get_value,get_variable,                          &
114       ccost,rr,bx,by,ax,ay,dx,dpx,dy,dpy,pi,                            &
115       sigx2,sigy2,ddx2,ddy2,disigh2,sigh2,fact,                         &
116       um,DGAUSS,piwint,litousch,tlitouschek,litouschw,                  &
117       dels,s1,s2,dx1,dx2,dpx1,dpx2,ax1,                                 &
118       ax2,bx1,bx2,beta2,gamma2,tol,fa0,fa1,fa2,                         &
119       dy1,dy2,dpy1,dpy2,ay1,ay2,by1,by2,sdum, half,                     &
120       tltouschek, tltouschekhr,pi2,km,ftousch,uloss
121
122
123  external ftousch,dgauss
124
125  parameter(half=5d-1)
126  character(name_len) name,sequ_name
127
128  pi=get_variable('pi ')
129  tlitouschek=0d0
130
131  ! ************* Get the parameters for the common blocks *************
132  ! *************         /machin/ and /beamdb/            *************
133
134  lp = get_string('beam ', 'particle ', sequ_name)
135
136  charge   = get_value('probe ', 'charge ')
137  gammas   = get_value('probe ', 'gamma ')
138  gamma    = get_value('probe ', 'gamma ')
139  en0      = get_value('probe ', 'energy ')
140  amass    = get_value('probe ', 'mass ')
141  ex       = get_value('probe ', 'ex ')
142  ey       = get_value('probe ', 'ey ')
143  et       = get_value('probe ', 'et ')
144  sigt     = get_value('probe ', 'sigt ')
145  sige     = get_value('probe ', 'sige ')
146  parnum   = get_value('probe ', 'npart ')
147  circ     = get_value('probe ', 'circ ')
148  currnt   = get_value('probe ', 'bcurrent ')
149  betas    = get_value('probe ', 'beta ')
150  beta     = get_value('probe ', 'beta ')
151  clight   = get_variable('clight ')
152  arad     = get_value('probe ', 'arad ')
153  alfa     = get_value('probe ', 'alfa ')
154  freq0    = get_value('probe ', 'freq0 ')
155  bunch    = get_value('probe ', 'kbunch ')
156  deltap   = get_value('probe ','deltap ')
157
158  !  ****************** Test print ********
159  print *, '  '
160  print *, 'TOUSCHEK MODULE '
161  print *, 'particle ', sequ_name(:lp)
162  print *, 'Charge ', charge
163  print *, 'gammas ', gammas
164  print *, 'gamma ', gamma
165  print *, 'Energy ', en0
166  print *, 'Mass ', amass
167  print *, 'Ex ', ex
168  print *, 'Ey ', ey
169  print *, 'Et ', et
170  print *, 'sigt ', sigt
171  print *, 'sige ', sige
172  print *, 'parnum ', parnum
173  print *, 'circ ', circ
174  print *, 'currnt ', currnt
175  print *, 'betas ', betas
176  print *, 'beta ', beta
177  print *, 'clight ', clight
178  print *, 'arad ', arad
179  print *, 'alfa ', alfa
180  print *, 'freq0 ', freq0
181  print *, 'kbunch ', bunch
182  print *, 'deltap ', deltap
183  ! ***************************************
184
185  rr = arad*arad
186  beta2 = beta*beta
187  gamma2 = gamma*gamma
188  ccost = rr*clight*parnum/(8d0*sqrt(pi)*gamma2*gamma2*beta2)
189  tol  = get_value('touschek ', 'tolerance ')
190  n    = get_option('touschek_table ')
191
192  ! ****** Start new Twiss Table reading *****************
193  !
194  call table_range('twiss ', '#s/#e ', range)
195  !      print *, 'Range for Table ', range(1), range(2)
196  flag = double_from_table_row('twiss ', 's ', range(1), s1)
197  if (flag .ne. 0)  goto 102
198  flag = double_from_table_row('twiss ', 'betx ', range(1), bx1)
199  if (flag .ne. 0)  goto 102
200  flag = double_from_table_row('twiss ', 'bety ', range(1), by1)
201  if (flag .ne. 0)  goto 102
202  flag = double_from_table_row('twiss ', 'alfx ', range(1), ax1)
203  if (flag .ne. 0)  goto 102
204  flag = double_from_table_row('twiss ', 'alfy ', range(1), ay1)
205  if (flag .ne. 0)  goto 102
206  flag = double_from_table_row('twiss ', 'dx ', range(1), dx1)
207  if (flag .ne. 0)  goto 102
208  flag = double_from_table_row('twiss ', 'dpx ', range(1), dpx1)
209  if (flag .ne. 0)  goto 102
210  flag = double_from_table_row('twiss ', 'dy ', range(1), dy1)
211  if (flag .ne. 0)  goto 102
212  flag = double_from_table_row('twiss ', 'dpy ', range(1), dpy1)
213  if (flag .ne. 0)  goto 102
214
215  ! ********** Start Do loop ***************
216  !
217  j = restart_sequ()
218  do i = range(1)+1, range(2)
219     j = advance_to_pos('twiss ', i)
220     flag = string_from_table_row('twiss ', 'name ', i, name)
221     flag = double_from_table_row('twiss ', 's ', i, s2)
222     if (flag .ne. 0)  goto 102
223     flag = double_from_table_row('twiss ', 'betx ', i, bx2)
224     if (flag .ne. 0)  goto 102
225     flag = double_from_table_row('twiss ', 'bety ', i, by2)
226     if (flag .ne. 0)  goto 102
227     flag = double_from_table_row('twiss ', 'alfx ', i, ax2)
228     if (flag .ne. 0)  goto 102
229     flag = double_from_table_row('twiss ', 'alfy ', i, ay2)
230     if (flag .ne. 0)  goto 102
231     flag = double_from_table_row('twiss ', 'dx ', i, dx2)
232     if (flag .ne. 0)  goto 102
233     flag = double_from_table_row('twiss ', 'dpx ', i, dpx2)
234     if (flag .ne. 0)  goto 102
235     flag = double_from_table_row('twiss ', 'dy ', i, dy2)
236     if (flag .ne. 0)  goto 102
237     flag = double_from_table_row('twiss ', 'dpy ', i, dpy2)
238     if (flag .ne. 0)  goto 102
239
240     dels = s2-s1
241     sdum = half * (s2 + s1)
242     bx  = half * (bx2 + bx1)
243     by  = half * (by2 + by1)
244     ax    = half * (ax2 + ax1)
245     ay    = half * (ay2 + ay1)
246     dx     = half * (dx2 + dx1)
247     dpx    = half * (dpx2 + dpx1)
248     dy     = half * (dy2 + dy1)
249     dpy    = half * (dpy2 + dpy1)
250
251     sigx2 = ex*bx
252     sigy2 = ey*by
253
254     ddx2 = (dpx*bx+dx*ax)**2
255     ddy2 = (dpy*by+dy*ay)**2
256
257     disigh2 = (1d0/sige**2)+((dx**2+ddx2)/sigx2)+((dy**2+ddy2)/sigy2)
258     sigh2 = (1d0/disigh2)
259     fact = sqrt(sigh2)*bx*by/(sigt*sige*sigx2*sigy2)
260     fa0 = 2d0*beta2*gamma2
261     fa1 = (bx**2/sigx2)*(1d0-sigh2*ddx2/sigx2)
262     fa2 = (by**2/sigy2)*(1d0-sigh2*ddy2/sigy2)
263     fb1 = (fa1 + fa2)/fa0
264     fb2 = sqrt(fb1**2-( bx**2*by**2*sigh2/(beta2*beta2*gamma2       &
265          *gamma2*sigx2*sigy2)*((1d0/sige**2)+(dx**2/sigx2)                 &
266          +(dy**2/sigy2))))
267
268     call cavtouschek(um,uloss,iflag)
269     um1 = um
270
271     if (um1.eq.0) then
272        call aawarn('TOUSCHEK ', '  rf voltage = 0, rest skipped ')
273        goto 101
274     endif
275
276     if (iflag .eq.1) then
277        call aawarn('TOUSCHEK ', ' uloss = 0 missing chrom in twiss ')
278     endif
279
280     !---- calculates the Piwinski integral
281
282     km = ATAN(sqrt(um1))
283
284     pi2 = pi/2.d0
285
286     piwint = DGAUSS(ftousch,km,pi2,tol)
287     litousch = ccost*fact*piwint
288     litouschw = litousch*dels/circ
289
290     !---- Accumulate contributions.
291     tlitouschek = tlitouschek + litousch*dels/circ
292
293     ! *************** Fill "touschek_table"  *********************
294
295     if(n.ne.0) then
296        call string_to_table_curr('touschek ', 'name ', name )
297        call double_to_table_curr('touschek ','s ', sdum)
298        call double_to_table_curr('touschek ','tli ', litousch)
299        call double_to_table_curr('touschek ','tliw ', litouschw)
300        call double_to_table_curr('touschek ','tlitot ', tlitouschek)
301        call augment_count('touschek ')
302     endif
303
304     s1   = s2
305     bx1  = bx2
306     by1  = by2
307     ax1  = ax2
308     ay1  = ay2
309     dx1  = dx2
310     dpx1 = dpx2
311
312  enddo
313  goto 101
314102 continue
315  call aawarn('TOUSCHEK ', 'table value not found, rest skipped ')
316101 continue
317
318  tltouschek=1d0/tlitouschek
319  tltouschekhr=1d0/tlitouschek/60d0/60d0
320  print *, '  '
321  print *, 'Energy radiated per turn [MeV]'
322  print *, uloss
323  print *, 'Touschek Lifetime [seconds/hours]'
324  print *, tltouschek, tltouschekhr
325
326  RETURN
327end subroutine touschek
328! ***************************************************************
329double precision function ftousch(k)
330
331  use physconsfi
332  use touschekfi
333  implicit none
334
335  double precision  k, pi, get_variable,z,aftoush
336  double precision ZR,ZI,BJOR,BJOI,BJIR,BJII,                       &
337       BYOR,BYOI,BYIR,BYII
338
339  integer  iflag
340
341  z = TAN(k)**2
342  pi=get_variable('pi ')
343
344  ZI = fb2*z
345  ZR = 0.
346
347  call CJYDBB(ZR,ZI,BJOR,BJOI,BJIR,BJII,                            &
348       BYOR,BYOI,BYIR,BYII,IFLAG)
349
350  aftoush = 2d0*( ((2d0*z+1d0)**2*(z/um1/(1d0+z)-1d0)/z)+z-         &
351       sqrt(z*um1*(1d0+z))-(2d0+1d0/(2d0*z))*log(z/um1/(1d0+z)) )        &
352       *sqrt(1d0+z)
353
354  if (iflag.eq.0) then
355     ftousch = aftoush*exp(-fb1*z)*BJOR
356  else
357     ftousch = aftoush*BJOR*(exp(-(fb1-fb2)*z)+                      &
358          exp(-(fb1+fb2)*z))/2.d0
359  end if
360
361  return
362end function ftousch
363
364! *********************************************************************
365!
366! $Id: touschek.f90,v 1.4 2010-03-13 02:25:54 frs Exp $
367!
368! $Log: not supported by cvs2svn $
369! Revision 1.3  2009/07/01 17:50:07  frs
370! Introducing elements:
371! placeholder 38 equivalent with instrument 24
372! tkicker 39 equivalent with kicker 15
373!
374! Revision 1.2  2009/04/06 23:30:07  frs
375! Fortran Clean-up: indenting, remove potentially uninitialized variable and
376! also remove unused variables
377!
378! Revision 1.1  2009/03/27 09:37:53  frs
379! New files needed for MAD-X Version 4
380!
381! Revision 1.8  2009/03/21 00:03:43  frs
382! Clean-up
383!
384! Revision 1.7  2008/03/07 10:09:53  frankz
385! add ion charge in rf bucket height calculation of Touschek module
386!
387! Revision 1.6  2005/08/24 13:42:19  frs
388! simple clean-up
389!
390! Revision 1.5  2005/02/05 00:54:43  frs
391! Fortran clean-up with mymod, f2c, lf95, f95
392!
393! Revision 1.4  2005/01/31 09:13:53  frs
394! Well tested touschek module
395!
396! Revision 1.1.1.1  1996/02/15 17:48:17  mclareni
397! Kernlib
398!
399!
400!#include "kernnum/pilot.h"
401DOUBLE PRECISION FUNCTION DGAUSS(F,A,B,EPS)
402  implicit none
403  integer i,LGFILE
404  DOUBLE PRECISION F,A,B,EPS
405  DOUBLE PRECISION W(12),X(12),AA,BB,C1,C2,U,S8,S16,CONST
406  LOGICAL MFLAG,RFLAG
407  EXTERNAL F
408  !
409  !     ******************************************************************
410  !
411  !     ADAPTIVE DOUBLE PRECISION GAUSSIAN QUADRATURE.
412  !
413  !     DGAUSS IS SET EQUAL TO THE APPROXIMATE VALUE OF THE INTEGRAL OF
414  !     THE FUNCTION F OVER THE INTERVAL (A,B), WITH ACCURACY PARAMETER
415  !     EPS.
416  !
417  !     ******************************************************************
418  !
419  !
420  DATA W / 0.1012285362903762591525313543D0,                        &
421       0.2223810344533744705443559944D0,                                 &
422       0.3137066458778872873379622020D0,                                 &
423       0.3626837833783619829651504493D0,                                 &
424       0.2715245941175409485178057246D-1,                                &
425       0.6225352393864789286284383699D-1,                                &
426       0.9515851168249278480992510760D-1,                                &
427       0.1246289712555338720524762822D0,                                 &
428       0.1495959888165767320815017305D0,                                 &
429       0.1691565193950025381893120790D0,                                 &
430       0.1826034150449235888667636680D0,                                 &
431       0.1894506104550684962853967232D0/
432  !
433  DATA X / 0.9602898564975362316835608686D0,                        &
434       0.7966664774136267395915539365D0,                                 &
435       0.5255324099163289858177390492D0,                                 &
436       0.1834346424956498049394761424D0,                                 &
437       0.9894009349916499325961541735D0,                                 &
438       0.9445750230732325760779884155D0,                                 &
439       0.8656312023878317438804678977D0,                                 &
440       0.7554044083550030338951011948D0,                                 &
441       0.6178762444026437484466717640D0,                                 &
442       0.4580167776572273863424194430D0,                                 &
443       0.2816035507792589132304605015D0,                                 &
444       0.9501250983763744018531933543D-1/
445  !
446  !     ******************************************************************
447  !
448  !  START.
449  DGAUSS=0.0D0
450  IF(B.EQ.A) RETURN
451  CONST=0.005D0/(B-A)
452  BB=A
453  !
454  !  COMPUTATIONAL LOOP.
4551 AA=BB
456  BB=B
4572 C1=0.5D0*(BB+AA)
458  C2=0.5D0*(BB-AA)
459  S8=0.0D0
460  DO I=1,4
461     U=C2*X(I)
462     S8=S8+W(I)*(F(C1+U)+F(C1-U))
463  enddo
464  S8=C2*S8
465  S16=0.0D0
466  DO I=5,12
467     U=C2*X(I)
468     S16=S16+W(I)*(F(C1+U)+F(C1-U))
469  enddo
470  S16=C2*S16
471  IF( ABS(S16-S8) .LE. EPS*(1.+ABS(S16)) ) GO TO 5
472  BB=C1
473  IF( 1.D0+ABS(CONST*C2) .NE. 1.D0) GO TO 2
474  DGAUSS=0.0D0
475  CALL KERMTR('D103.1',LGFILE,MFLAG,RFLAG)
476  IF(MFLAG) THEN
477     IF(LGFILE.EQ.0) THEN
478        WRITE(*,6)
479     ELSE
480        WRITE(LGFILE,6)
481     ENDIF
482  ENDIF
483  IF(.NOT. RFLAG) CALL ABEND
484  RETURN
4855 DGAUSS=DGAUSS+S16
486  IF(BB.NE.B) GO TO 1
487  RETURN
488  !
4896 FORMAT( 4X, 'FUNCTION DGAUSS ... TOO HIGH ACCURACY REQUIRED')
490END FUNCTION DGAUSS
491
492!********************************************************************
493!
494! $Id: touschek.f90,v 1.4 2010-03-13 02:25:54 frs Exp $
495!
496! $Log: not supported by cvs2svn $
497! Revision 1.3  2009/07/01 17:50:07  frs
498! Introducing elements:
499! placeholder 38 equivalent with instrument 24
500! tkicker 39 equivalent with kicker 15
501!
502! Revision 1.2  2009/04/06 23:30:07  frs
503! Fortran Clean-up: indenting, remove potentially uninitialized variable and
504! also remove unused variables
505!
506! Revision 1.1  2009/03/27 09:37:53  frs
507! New files needed for MAD-X Version 4
508!
509! Revision 1.8  2009/03/21 00:03:43  frs
510! Clean-up
511!
512! Revision 1.7  2008/03/07 10:09:53  frankz
513! add ion charge in rf bucket height calculation of Touschek module
514!
515! Revision 1.6  2005/08/24 13:42:19  frs
516! simple clean-up
517!
518! Revision 1.5  2005/02/05 00:54:43  frs
519! Fortran clean-up with mymod, f2c, lf95, f95
520!
521! Revision 1.4  2005/01/31 09:13:53  frs
522! Well tested touschek module
523!
524! Revision 1.1.1.1  1996/02/15 17:48:35  mclareni
525! Kernlib
526!
527!
528SUBROUTINE KERSET(ERCODE,LGFILE,LIMITM,LIMITR)
529  IMPLICIT NONE
530
531  INTEGER KOUNTE,L,I,LIMITM,LIMITR,LOG,LOGF,LGFILE
532  PARAMETER(KOUNTE  =  27)
533  CHARACTER(6)        ERCODE,   CODE(KOUNTE)
534  LOGICAL             MFLAG,    RFLAG
535  INTEGER             KNTM(KOUNTE),       KNTR(KOUNTE)
536  DATA      LOGF      /  0  /
537  DATA      CODE(1), KNTM(1), KNTR(1)  / 'C204.1', 255, 255 /
538  DATA      CODE(2), KNTM(2), KNTR(2)  / 'C204.2', 255, 255 /
539  DATA      CODE(3), KNTM(3), KNTR(3)  / 'C204.3', 255, 255 /
540  DATA      CODE(4), KNTM(4), KNTR(4)  / 'C205.1', 255, 255 /
541  DATA      CODE(5), KNTM(5), KNTR(5)  / 'C205.2', 255, 255 /
542  DATA      CODE(6), KNTM(6), KNTR(6)  / 'C305.1', 255, 255 /
543  DATA      CODE(7), KNTM(7), KNTR(7)  / 'C308.1', 255, 255 /
544  DATA      CODE(8), KNTM(8), KNTR(8)  / 'C312.1', 255, 255 /
545  DATA      CODE(9), KNTM(9), KNTR(9)  / 'C313.1', 255, 255 /
546  DATA      CODE(10),KNTM(10),KNTR(10) / 'C336.1', 255, 255 /
547  DATA      CODE(11),KNTM(11),KNTR(11) / 'C337.1', 255, 255 /
548  DATA      CODE(12),KNTM(12),KNTR(12) / 'C341.1', 255, 255 /
549  DATA      CODE(13),KNTM(13),KNTR(13) / 'D103.1', 255, 255 /
550  DATA      CODE(14),KNTM(14),KNTR(14) / 'D106.1', 255, 255 /
551  DATA      CODE(15),KNTM(15),KNTR(15) / 'D209.1', 255, 255 /
552  DATA      CODE(16),KNTM(16),KNTR(16) / 'D509.1', 255, 255 /
553  DATA      CODE(17),KNTM(17),KNTR(17) / 'E100.1', 255, 255 /
554  DATA      CODE(18),KNTM(18),KNTR(18) / 'E104.1', 255, 255 /
555  DATA      CODE(19),KNTM(19),KNTR(19) / 'E105.1', 255, 255 /
556  DATA      CODE(20),KNTM(20),KNTR(20) / 'E208.1', 255, 255 /
557  DATA      CODE(21),KNTM(21),KNTR(21) / 'E208.2', 255, 255 /
558  DATA      CODE(22),KNTM(22),KNTR(22) / 'F010.1', 255,   0 /
559  DATA      CODE(23),KNTM(23),KNTR(23) / 'F011.1', 255,   0 /
560  DATA      CODE(24),KNTM(24),KNTR(24) / 'F012.1', 255,   0 /
561  DATA      CODE(25),KNTM(25),KNTR(25) / 'F406.1', 255,   0 /
562  DATA      CODE(26),KNTM(26),KNTR(26) / 'G100.1', 255, 255 /
563  DATA      CODE(27),KNTM(27),KNTR(27) / 'G100.2', 255, 255 /
564  LOGF  =  LGFILE
565
566  L  =  0
567  IF(ERCODE .NE. ' ')  THEN
568     DO  L = 1, 6
569        IF(ERCODE(1:L) .EQ. ERCODE)  GOTO 12
570     enddo
57112   CONTINUE
572  ENDIF
573  DO     I  =  1, KOUNTE
574     IF(L .EQ. 0)  GOTO 13
575     IF(CODE(I)(1:L) .NE. ERCODE(1:L))  GOTO 14
57613   IF(LIMITM.GE.0) KNTM(I)  =  LIMITM
577     IF(LIMITR.GE.0) KNTR(I)  =  LIMITR
57814   CONTINUE
579  enddo
580  RETURN
581
582  ENTRY KERMTR(ERCODE,LOG,MFLAG,RFLAG)
583  LOG  =  LOGF
584  DO     I  =  1, KOUNTE
585     IF(ERCODE .EQ. CODE(I))  GOTO 21
586  enddo
587  WRITE(*,1000)  ERCODE
588  CALL ABEND
589  RETURN
59021 RFLAG  =  KNTR(I) .GE. 1
591  IF(RFLAG  .AND.  (KNTR(I) .LT. 255))  KNTR(I)  =  KNTR(I) - 1
592  MFLAG  =  KNTM(I) .GE. 1
593  IF(MFLAG  .AND.  (KNTM(I) .LT. 255))  KNTM(I)  =  KNTM(I) - 1
594  IF(.NOT. RFLAG)  THEN
595     IF(LOGF .LT. 1)  THEN
596        WRITE(*,1001)  CODE(I)
597     ELSE
598        WRITE(LOGF,1001)  CODE(I)
599     ENDIF
600  ENDIF
601  IF(MFLAG .AND. RFLAG)  THEN
602     IF(LOGF .LT. 1)  THEN
603        WRITE(*,1002)  CODE(I)
604     ELSE
605        WRITE(LOGF,1002)  CODE(I)
606     ENDIF
607  ENDIF
608  RETURN
6091000 FORMAT(' KERNLIB LIBRARY ERROR. ' /                               &
610       ' ERROR CODE ',A6,' NOT RECOGNIZED BY KERMTR',                    &
611       ' ERROR MONITOR. RUN ABORTED.')
6121001 FORMAT(/' ***** RUN TERMINATED BY CERN LIBRARY ERROR ',           &
613       'CONDITION ',A6)
6141002 FORMAT(/' ***** CERN LIBRARY ERROR CONDITION ',A6)
615END SUBROUTINE KERSET
616
617!_________________________________________________
618!
619!     THE CJYDBB  ROUTINE
620!         given the arguments:  (Re=ZR, Im=ZI),  ZI can be .GT. 150
621!         returns the modified Bessel functions:
622!                J0 (Re=BJOR,Im= BJOI), J1(Re=BJIR, Im=BJII)
623!       Y0 (Re=BYOR,Im= BYOI), Y1(Re=BYIR, Im=BYII)
624!         if ZI .GT. 150 then iflag = 1 and   BJOR = BJOR/cosh(ZI)
625!_________________________________________________
626
627SUBROUTINE CJYDBB(ZR,ZI,BJOR,BJOI,BJIR,BJII,                      &
628     BYOR,BYOI,BYIR,BYII,IFLAG)
629
630  IMPLICIT NONE
631
632  INTEGER iflag,k,j,l,n,m
633  DOUBLE PRECISION ZR,ZI,BJOR,BJOI,BJIR,BJII,                       &
634       BYOR,BYOI,BYIR,BYII,zmag,angz,cang,                               &
635       cangz,sang,sangz,aabs,cterm0,cterm2,                              &
636       ang,u,ui,angu,ur,szi,zzi,hcosx,hsinx,                             &
637       zimag,sxir,cxir,sxor,cxor,ti,t,zir,zor,                           &
638       z2mag,angt,psi,fk,rk,cterm,cterm1,                                &
639       bjiit,byort,byoit,byirt,byiit,pi2,gam,tyr,                        &
640       tyi,sum1,sum0,cmu0,cmu1,tk2,s,smag,                               &
641       spor,spir,spii,sqor,sqoi,sqir,sqii,                               &
642       por,poi,spoi,pii,pi4,htanx,byiti,byitr,                           &
643       bjiti,bjitr,byoti,byotr,bjoti,bjotr,                              &
644       sii,sir,cii,cir,soi,sor,coi,cor,z1,qii,                           &
645       bjirt,rl,angs,tr,pir,qor,qoi,qir,sink
646
647  DOUBLE PRECISION CJOR(18),CJOI(18),CJIR(18),CJII(18),             &
648       CYOR(18),CYOI(18),CYIR(18),CYII(18),CPO(12),                      &
649       CPI(12),CQO(12),CQI(12),CPOR(12),CPOI(12),                        &
650       CPIR(12),CPII(12),CQOR(12),CQOI(12),                              &
651       CQIR(12),CQII(12)
652
653  !---- Initialize.
654  call dzero(CJOR,18)
655  call dzero(CJOI,18)
656  call dzero(CJIR,18)
657  call dzero(CJII,18)
658  call dzero(CYOR,18)
659  call dzero(CYOI,18)
660  call dzero(CYIR,18)
661  call dzero(CYII,18)
662  call dzero(CPO,12)
663  call dzero(CPI,12)
664  call dzero(CQO,12)
665  call dzero(CQI,12)
666  call dzero(CPOR,12)
667  call dzero(CPOI,12)
668  call dzero(CPIR,12)
669  call dzero(CPII,12)
670  call dzero(CQOR,12)
671  call dzero(CQOI,12)
672  call dzero(CQIR,12)
673  call dzero(CQII,12)
674  iflag = 0
675
676  Z2MAG=ZR*ZR+ZI*ZI
677  ZMAG=SQRT(Z2MAG)
678  ANGZ=ATAN2(ZI,ZR)
679  IF(ZMAG - .1D+02.gt.0) goto 21
680  T=Z2MAG/.4D+01
681  ANGT=ANGZ*.2D+01
682  PSI=0.0D+00
683  FK=.1D+01
684  DO K=1,18
685     RK=DBLE(FLOAT(K))
686     ANG=ANGT*RK
687     SINK=(-.1D+01)**K
688     FK=FK*RK
689     SANG=SIN(ANG)
690     CANG=COS(ANG)
691     CTERM=SINK/(FK*FK)
692     CJOR(K)=CTERM*CANG
693     CJOI(K)=CTERM*SANG
694     CTERM1=CTERM/(RK+.1D+01)
695     CJIR(K)=CTERM1*CANG
696     CJII(K)=CTERM1*SANG
697     PSI=PSI+.1D+01/RK
698     CTERM0=CTERM*PSI
699     CYOR(K)=CTERM0*CANG
700     CYOI(K)=CTERM0*SANG
701     CTERM2=CTERM1*(PSI+.5D+00/(RK+.1D+01))
702     CYIR(K)=CTERM2*CANG
703     CYII(K)=CTERM2*SANG
704  enddo
705  BJOR=(((((((((((((((((CJOR(18)*T+CJOR(17))*T+                     &
706       CJOR(16))*T+CJOR(15))*T+CJOR(14))*T+CJOR(13))*T                   &
707       +CJOR(12))*T+CJOR(11))*T+CJOR(10))*T+CJOR(9))*T                   &
708       +CJOR(8))*T+CJOR(7))*T+CJOR(6))*T+CJOR(5))*T+                     &
709       CJOR(4))*T+CJOR(3))*T+CJOR(2))*T+CJOR(1))*T+                      &
710       .1D+01
711  BJOI=(((((((((((((((((CJOI(18)*T+CJOI(17))*T+                     &
712       CJOI(16))*T+CJOI(15))*T+CJOI(14))*T+CJOI(13))*T                   &
713       +CJOI(12))*T+CJOI(11))*T+CJOI(10))*T+CJOI(9))*T                   &
714       +CJOI(8))*T+CJOI(7))*T+CJOI(6))*T+CJOI(5))*T+                     &
715       CJOI(4))*T+CJOI(3))*T+CJOI(2))*T+CJOI(1))*T
716  BJIRT=(((((((((((((((((CJIR(18)*T+CJIR(17))*T+                    &
717       CJIR(16))*T+CJIR(15))*T+CJIR(14))*T+CJIR(13))*T                   &
718       +CJIR(12))*T+CJIR(11))*T+CJIR(10))*T+CJIR(9))*T                   &
719       +CJIR(8))*T+CJIR(7))*T+CJIR(6))*T+CJIR(5))*T+                     &
720       CJIR(4))*T+CJIR(3))*T+CJIR(2))*T+CJIR(1))*T+                      &
721       .1D+01
722  BJIIT=(((((((((((((((((CJII(18)*T+CJII(17))*T+CJII(16))           &
723       *T+CJII(15))*T+CJII(14))*T+CJII(13))*T+CJII(12))*T                &
724       +CJII(11))*T+CJII(10))*T+CJII(9))*T+CJII(8))*T+                   &
725       CJII(7))*T+CJII(6))*T+CJII(5))*T+CJII(4))*T+                      &
726       CJII(3))*T+CJII(2))*T+CJII(1))*T
727  BYORT=(((((((((((((((((CYOR(18)*T+CYOR(17))*T+CYOR(16))           &
728       *T+CYOR(15))*T+CYOR(14))*T+CYOR(13))*T+CYOR(12))*T+               &
729       CYOR(11))*T+CYOR(10))*T+CYOR(9))*T+CYOR(8))*T+                    &
730       CYOR(7))*T+CYOR(6))*T+CYOR(5))*T+CYOR(4))*T+CYOR(3))              &
731       *T+CYOR(2))*T+CYOR(1))*T
732  BYOIT=(((((((((((((((((CYOI(18)*T+CYOI(17))*T+CYOI(16))*T         &
733       +CYOI(15))*T+CYOI(14))*T+CYOI(13))*T+CYOI(12))*T+                 &
734       CYOI(11))*T+CYOI(10))*T+CYOI(9))*T+CYOI(8))*T+CYOI(7)             &
735       )*T+CYOI(6))*T+CYOI(5))*T+CYOI(4))*T+CYOI(3))*                    &
736       T+CYOI(2))*T+CYOI(1))*T
737  BYIRT=(((((((((((((((((CYIR(18)*T+CYIR(17))*T+CYIR(16))*T+        &
738       CYIR(15))*T+CYIR(14))*T+CYIR(13))*T+CYIR(12))*T+                  &
739       CYIR(11))*T+CYIR(10))*T+CYIR(9))*T+CYIR(8))*T+CYIR(7)             &
740       )*T+CYIR(6))*T+CYIR(5))*T+CYIR(4))*T+CYIR(3))*T+                  &
741       CYIR(2))*T+CYIR(1))*T+.5D+00
742  BYIIT=(((((((((((((((((CYII(18)*T+CYII(17))*T+CYII(16))*T         &
743       +CYII(15))*T+CYII(14))*T+CYII(13))*T+CYII(12))*T+                 &
744       CYII(11))*T+CYII(10))*T+CYII(9))*T+CYII(8))*T+CYII(7)             &
745       )*T+CYII(6))*T+CYII(5))*T+CYII(4))*T+CYII(3))*T+                  &
746       CYII(2))*T+CYII(1))*T
747  CANGZ=COS(ANGZ)
748  SANGZ=SIN(ANGZ)
749  BJIR=(ZR*BJIRT-ZI*BJIIT)/.2D+01
750  BJII=(ZI*BJIRT+ZR*BJIIT)/.2D+01
751  PI2=.2D+01*ATAN(.1D+01)
752  GAM=.577215664901533D+00
753  TYR=LOG(ZMAG/.2D+01)+GAM
754  TYI=ANGZ
755  BYOR=(TYR*BJOR-TYI*BJOI-BYORT)/PI2
756  BYOI=(TYI*BJOR+TYR*BJOI-BYOIT)/PI2
757  BYIR=(TYR*BJIR-TYI*BJII-CANGZ/ZMAG-(BYIRT*ZR-                     &
758       BYIIT*ZI)/.2D+01)/PI2
759  BYII=(TYI*BJIR+TYR*BJII+SANGZ/ZMAG-(BYIIT*ZR+BYIRT*ZI)/.2D+01)/   &
760       &PI2
761  RETURN
76221 FK=.1D+01
763  if (ZR .ne.0) then
764     AABS=ABS(ZI/ZR)
765     ANG=ATAN(AABS)
766  else
767     ang = ASIN(1d0)
768  end if
769  IF(ZI.lt.0) goto 50
770  ANGZ=ANG
771  GOTO 52
77250 ANGZ=-ANG
77352 CONTINUE
774  SUM0=.1D+01
775  SUM1=.1D+01
776  CMU0=0.0
777  CMU1=.4D+01
778  L=1
779  DO K=1,24
780     L=-L
781     RK=DBLE(FLOAT(K))
782     FK=FK
783     TK2=4.*RK-4.+1./RK
784     SUM0=SUM0*(CMU0-TK2)
785     SUM1=SUM1*(CMU1/RK-TK2)
786     IF(L.gt.0) goto 23
787     N=(K+1)/2
788     CQO(N)=SUM0/FK
789     CQI(N)=SUM1/FK
790     GOTO 30
79123   M=K/2
792     CPO(M)=SUM0/FK
793     CPI(M)=SUM1/FK
79430   CONTINUE
795  enddo
796  RL=-.1D+01
797  DO J=1,12
798     RL=-RL
799     CQO(J)=CQO(J)*RL
800     CQI(J)=CQI(J)*RL
801     CPO(J)=CPO(J)*RL
802     CPI(J)=CPI(J)*RL
803  enddo
804  ANGS=-ANGZ*.2D+01
805  SMAG=.1D+01/(.64D+02*Z2MAG)
806  S=SMAG
807  DO L=1,12
808     RL=DBLE(FLOAT(L))
809     ANG=ANGS*RL
810     CANG=COS(ANG)
811     SANG=SIN(ANG)
812     CPOR(L)=CPO(L)*CANG
813     CPOI(L)=CPO(L)*SANG
814     CPIR(L)=CPI(L)*CANG
815     CPII(L)=CPI(L)*SANG
816     CQOR(L)=CQO(L)*CANG
817     CQOI(L)=CQO(L)*SANG
818     CQIR(L)=CQI(L)*CANG
819     CQII(L)=CQI(L)*SANG
820  enddo
821  SPOR=(((((((((((CPOR(12)*S+CPOR(11))*S+CPOR(10))*S+               &
822       CPOR(9))*S+CPOR(8))*S+CPOR(7))*S+CPOR(6))*S+CPOR(5))              &
823       *S+CPOR(4))*S+CPOR(3))*S+CPOR(2))*S+CPOR(1))*S
824  SPOI=(((((((((((CPOI(12)*S+CPOI(11))*S+CPOI(10))*S+               &
825       CPOI(9))*S+CPOI(8))*S+CPOI(7))*S+CPOI(6))*S+CPOI(5))              &
826       *S+CPOI(4))*S+CPOI(3))*S+CPOI(2))*S+CPOI(1))*S
827  SPIR=(((((((((((CPIR(12)*S+CPIR(11))*S+CPIR(10))*S+               &
828       CPIR(9))*S+CPIR(8))*S+CPIR(7))*S+CPIR(6))*S+CPIR(5))              &
829       *S+CPIR(4))*S+CPIR(3))*S+CPIR(2))*S+CPIR(1))*S
830  SPII=(((((((((((CPII(12)*S+CPII(11))*S+CPII(10))*S+               &
831       CPII(9))*S+CPII(8))*S+CPII(7))*S+CPII(6))*S+CPII(5))              &
832       *S+CPII(4))*S+CPII(3))*S+CPII(2))*S+CPII(1))*S
833  SQOR=(((((((((((CQOR(12)*S+CQOR(11))*S+CQOR(10))*S+               &
834       CQOR(9))*S+CQOR(8))*S+CQOR(7))*S+CQOR(6))*S+CQOR(5))              &
835       *S+CQOR(4))*S+CQOR(3))*S+CQOR(2))*S+CQOR(1))*S
836  SQOI=(((((((((((CQOI(12)*S+CQOI(11))*S+CQOI(10))*S+               &
837       CQOI(9))*S+CQOI(8))*S+CQOI(7))*S+CQOI(6))*S+CQOI(5))              &
838       *S+CQOI(4))*S+CQOI(3))*S+CQOI(2))*S+CQOI(1))*S
839  SQIR=(((((((((((CQIR(12)*S+CQIR(11))*S+CQIR(10))*S+               &
840       CQIR(9))*S+CQIR(8))*S+CQIR(7))*S+CQIR(6))*S+CQIR(5))              &
841       *S+CQIR(4))*S+CQIR(3))*S+CQIR(2))*S+CQIR(1))*S
842  SQII=(((((((((((CQII(12)*S+CQII(11))*S+CQII(10))*S+               &
843       CQII(9))*S+CQII(8))*S+CQII(7))*S+CQII(6))*S+CQII(5))              &
844       *S+CQII(4))*S+CQII(3))*S+CQII(2))*S+CQII(1))*S
845  T=.8D+01*ZMAG
846  ANGT=ANGZ
847  TR=T*COS(ANGT)
848  TI=T*SIN(ANGT)
849  POR=.1D+01-SPOR
850  POI=-SPOI
851  PIR=.1D+01-SPIR
852  PII=-SPII
853  QOR=TR*SQOR-TI*SQOI
854  QOI=TI*SQOR+TR*SQOI
855  QIR=TR*SQIR-TI*SQII
856  QII=TI*SQIR+TR*SQII
857  PI4=ATAN(.1D+01)
858  PI2=.2D+01*PI4
859  Z1=ABS(ZR)
860  ZOR=Z1-PI4
861  ZIR=Z1-(PI2+PI4)
862  ZZI=ZI
863  CXOR=COS(ZOR)
864  SXOR=SIN(ZOR)
865  CXIR=COS(ZIR)
866  SXIR=SIN(ZIR)
867  ZIMAG=ABS(ZZI)
868  !***
869  IF(ZIMAG .GT. .15D+03) then
870     !         WRITE(*,100)
871     !         WRITE(*,*) iflag,zimag
872     iflag = 1
873  end if
874100 FORMAT(1X,'ZIMAG .GT. .15D+03 (SUBROUTINE CJYDBB)')
875  !***
876  IF(ZIMAG - .15D+03.gt.0) goto 61
877  HSINX=(EXP(ZZI)-EXP(-ZZI))/.2D+01
878  HCOSX=(EXP(ZZI)+EXP(-ZZI))/.2D+01
879  COR=CXOR*HCOSX
880  COI=-SXOR*HSINX
881  SOR=SXOR*HCOSX
882  SOI=CXOR*HSINX
883  CIR=CXIR*HCOSX
884  CII=-SXIR*HSINX
885  SIR=SXIR*HCOSX
886  SII=CXIR*HSINX
887  GOTO 62
88861 SZI=SNGL(ZZI)
889  HTANX=DBLE(TANH(SZI))
890  COR=CXOR
891  COI=-SXOR*HTANX
892  SOR=SXOR
893  SOI=CXOR*HTANX
894  CIR=CXIR
895  CII=-SXIR*HTANX
896  SIR=SXIR
897  SII=CXIR*HTANX
89862 CONTINUE
899  BJOTR=(POR*COR-POI*COI)-(QOR*SOR-QOI*SOI)
900  BJOTI=(POI*COR+POR*COI)-(QOI*SOR+QOR*SOI)
901  BYOTR=(POR*SOR-POI*SOI)+(QOR*COR-QOI*COI)
902  BYOTI=(POI*SOR+POR*SOI)+(QOI*COR+QOR*COI)
903  BJITR=(PIR*CIR-PII*CII)-(QIR*SIR-QII*SII)
904  BJITI=(PII*CIR+PIR*CII)-(QII*SIR+QIR*SII)
905  BYITR=(PIR*SIR-PII*SII)+(QIR*CIR-QII*CII)
906  BYITI=(PII*SIR+PIR*SII)+(QII*CIR+QIR*CII)
907  U=SQRT(.1D+01/(PI2*ZMAG))
908  ANGU=-ANGZ/.2D+01
909  UR=U*COS(ANGU)
910  UI=U*SIN(ANGU)
911  BJOR=(UR*BJOTR-UI*BJOTI)
912  BJOI=(UI*BJOTR+UR*BJOTI)
913  BYOR=(UR*BYOTR-UI*BYOTI)
914  BYOI=(UI*BYOTR+UR*BYOTI)
915  BJIR=(UR*BJITR-UI*BJITI)
916  BJII=(UI*BJITR+UR*BJITI)
917  BYIR=(UR*BYITR-UI*BYITI)
918  BYII=(UI*BYITR+UR*BYITI)
919  IF(ZR.lt.0) goto 40
920  RETURN
92140 BYOR=BYOR+.2D+01*BJOI
922  BYOI=-BYOI+.2D+01*BJOR
923  BYIR=-BYIR-.2D+01*BJII
924  BYII=BYII-.2D+01*BJIR
925  BJOI=-BJOI
926  BJIR=-BJIR
927
928  RETURN
929END SUBROUTINE CJYDBB
930!   ____________
931!
Note: See TracBrowser for help on using the repository browser.