Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
dpm25nulib.f
Go to the documentation of this file.
1  SUBROUTINE timdat
2  IMPLICIT DOUBLE PRECISION (a-h,o-z)
3  character*8 datum
4  character*11 zeit
5 * CALL GETDAT(IYEAR,IMONTH,IDAY)
6 * CALL GETTIM(IHOUR,IMINUT,ISECND,IHSCND)
7 * WRITE(6,'(A)') ' YEAR, MONTH,DAY '
8 * WRITE(6,'(1X,3I4)')IYEAR,IMONTH,IDAY
9 * WRITE(6,'(A)') ' HOUR, MINUTE,SECOND,...'
10 * WRITE(6,'(1X,4I4)')IHOUR,IMINUT,ISECND,IHSCND
11 C call date(datum)
12 C call time(zeit)
13 C WRITE(6,'(4A)') ' YEAR / MONTH / DAY : ',datum,' time: ',zeit
14  RETURN
15  END
16 *CMZ : 1.00/00 25/11/91 16.46.51 by H.-J. M¶hring
17 *-- Author :
18 C===========================================================
19 C#######################################################################
20 C THIS IS A PACKAGE CONTAINING A RANDOM NUMBER GENERATOR AND
21 C SERVICE ROUTINES.
22 C THE ALGORITHM IS FROM
23 C 'TOWARD A UNVERSAL RANDOM NUMBER GENERATOR'
24 C G.MARSAGLIA, A.ZAMAN ; FLORIDA STATE UNIV. PREPRINT FSU-SCRI-87-
25 C IMPLEMENTATION BY K. HAHN DEC. 88,
26 C THIS GENERATOR SHOULD NOT DEPEND ON THE HARD WARE ( IF A REAL HAS
27 C AT LEAST 24 SIGNIFICANT BITS IN INTERNAL REPRESENTATION ),
28 C THE PERIOD IS ABOUT 2**144,
29 C TIME FOR ONE CALL AT IBM-XT IS ABOUT 0.7 MILLISECONDS,
30 C THE PACKAGE CONTAINS
31 C FUNCTION RNDM(I) : GENERATOR
32 C SUBROUTINE RNDMST(NA1,NA2,NA3,NB4) : INITIALIZATION
33 C SUBROUTINE RNDMIN(U,C,CD,CM,I,J) : PUT SEED TO GENERATOR
34 C SUBROUTINE RNDMOU(U,C,CD,CM,I,J) : TAKE SEED FROM GENERATOR
35 C SUBROUTINE RNDMTE(IO) : TEST OF GENERATOR
36 C---
37 C FUNCTION RNDM(I)
38 C GIVES UNIFORMLY DISTRIBUTED RANDOM NUMBERS IN 0..1)
39 C I - DUMMY VARIABLE, NOT USED
40 C SUBROUTINE RNDMST(NA1,NA2,NA3,NB1)
41 C INITIALIZES THE GENERATOR, MUST BE CALLED BEFORE USING RNDM
42 C NA1,NA2,NA3,NB1 - VALUES FOR INITIALIZING THE GENERATOR
43 C NA? MUST BE IN 1..178 AND NOT ALL 1
44 C 12,34,56 ARE THE STANDARD VALUES
45 C NB1 MUST BE IN 1..168
46 C 78 IS THE STANDARD VALUE
47 C SUBROUTINE RNDMIN(U,C,CD,CM,I,J)
48 C PUTS SEED TO GENERATOR ( BRINGS GENERATOR IN THE SAME STATUS
49 C AS AFTER THE LAST RNDMOU CALL )
50 C U(97),C,CD,CM,I,J - SEED VALUES AS TAKEN FROM RNDMOU
51 C SUBROUTINE RNDMOU(U,C,CD,CM,I,J)
52 C TAKES SEED FROM GENERATOR
53 C U(97),C,CD,CM,I,J - SEED VALUES
54 C SUBROUTINE RNDMTE(IO)
55 C TEST OF THE GENERATOR
56 C IO - DEFINES OUTPUT
57 C = 0 OUTPUT ONLY IF AN ERROR IS DETECTED
58 C = 1 OUTPUT INDEPENDEND ON AN ERROR
59 C RNDMTE USES RNDMIN AND RNDMOU TO BRING GENERATOR TO SAME STATUS
60 C AS BEFORE CALL OF RNDMTE
61 C#######################################################################
62 C===========================================================
63 C DOUBLE PRECISION FUNCTION RNDM(VDUMMY)
64 C IMPLICIT DOUBLE PRECISION (A-H,O-Z)
65 C COMMON /RANDOO/ U(97),C,CD,CM,I,J
66 C RNDM = U(I)-U(J)
67 C IF ( RNDM.LT.0.0D0 ) RNDM = RNDM+1.0D0
68 C U(I) = RNDM
69 C I = I-1
70 C IF ( I.EQ.0 ) I = 97
71 C J = J-1
72 C IF ( J.EQ.0 ) J = 97
73 C C = C-CD
74 C IF ( C.LT.0.0D0 ) C = C+CM
75 C RNDM = RNDM-C
76 C IF ( RNDM.LT.0.0D0 ) RNDM = RNDM+1.0D0
77 C RETURN
78 C END
79 *-- Author :
80 C===========================================================
81  SUBROUTINE rndmst(NA1,NA2,NA3,NB1)
82  IMPLICIT DOUBLE PRECISION (a-h,o-z)
83  COMMON /randoo/ u(97),c,cd,cm,i,j
84  ma1 = na1
85  ma2 = na2
86  ma3 = na3
87  mb1 = nb1
88  i = 97
89  j = 33
90  DO 20 ii2 = 1,97
91  s = 0
92  t = 0.5d0
93  DO 10 ii1 = 1,24
94  mat = mod(mod(ma1*ma2,179)*ma3,179)
95  ma1 = ma2
96  ma2 = ma3
97  ma3 = mat
98  mb1 = mod(53*mb1+1,169)
99  IF ( mod(mb1*mat,64).GE.32 ) s = s+t
100  10 t = 0.5*t
101  20 u(ii2) = s
102  c = 362436.d0/16777216.d0
103  cd = 7654321.d0/16777216.d0
104  cm = 16777213.d0/16777216.d0
105  RETURN
106  END
107 *-- Author :
108 C==========================================================
109  SUBROUTINE rndmin(UIN,CIN,CDIN,CMIN,IIN,JIN)
110  IMPLICIT DOUBLE PRECISION (a-h,o-z)
111  dimension uin(97)
112  COMMON /randoo/ u(97),c,cd,cm,i,j
113  DO 10 kkk = 1,97
114  10 u(kkk) = uin(kkk)
115  c = cin
116  cd = cdin
117  cm = cmin
118  i = iin
119  j = jin
120  RETURN
121  END
122 *-- Author :
123 C==========================================================
124  SUBROUTINE rndmou(UOUT,COUT,CDOUT,CMOUT,IOUT,JOUT)
125  IMPLICIT DOUBLE PRECISION (a-h,o-z)
126  dimension uout(97)
127  COMMON /randoo/ u(97),c,cd,cm,i,j
128  DO 10 kkk = 1,97
129  10 uout(kkk) = u(kkk)
130  cout = c
131  cdout = cd
132  cmout = cm
133  iout = i
134  jout = j
135  RETURN
136  END
137 *-- Author :
138 C==========================================================
139  SUBROUTINE rndmte(IO)
140  IMPLICIT DOUBLE PRECISION (a-h,o-z)
141  dimension uu(97),u(6),x(6),d(6)
142  DATA u / 6533892.d0, 14220222.d0, 7275067.d0, 6172232.d0,
143  +8354498.d0, 10633180.d0/
144  CALL rndmou(uu,cc,ccd,ccm,ii,jj)
145  CALL rndmst(12,34,56,78)
146  DO 10 ii1 = 1,20000
147  10 xx = rndm(v)
148  sd = 0
149  DO 20 ii2 = 1,6
150  x(ii2) = 4096.d0*(4096.d0*rndm(v))
151  d(ii2) = x(ii2)-u(ii2)
152  20 sd = sd+d(ii2)
153  CALL rndmin(uu,cc,ccd,ccm,ii,jj)
154 C IF ( IO.EQ. 1.OR. SD.NE.0. 0) WRITE(6,500) (U(I),X(I),D(I),I=1,6)
155  RETURN
156 C==END OF RANDOM GENERATOR PACKAGE==========================
157  500 FORMAT(' === TEST OF THE RANDOM-GENERATOR ===',/,
158  +' EXPECTED VALUE CALCULATED VALUE DIFFERENCE',/, 6(f17.
159  +1,f20.1,f15.3,/), ' === END OF TEST ;',
160  +' GENERATOR HAS THE SAME STATUS AS BEFORE CALLING RNDMTE')
161  END
162 *-- Author :
163 C 03/10/89 910141255 MEMBER NAME KKEVT (KK89.S) F77
164  SUBROUTINE rannor(X,Y)
165  IMPLICIT DOUBLE PRECISION (a-h,o-z)
166 C X UND Y SIND GAUSSVERTEILTE ZUFALLSZAHLEN
167  CALL dsfecf(sfe,cfe)
168  a=sqrt(-2.d0*log(rndm(v)))
169  x=a*sfe
170  y=a*cfe
171  RETURN
172  END
173 *-- Author :
174 ****************************************************** PLOT**********
175 *
176  SUBROUTINE plot (X,Y,N,M,MM,XO,DX,YO,DY)
177  IMPLICIT DOUBLE PRECISION (a-h,o-z)
178 *
179 * initial version
180 * J. Ranft, (FORTRAN-Programmierung,J.R.,Teubner, Leipzig, 72)
181 * This is a subroutine of fluka to plot Y across the page
182 * as a function of X down the page. Up to 37 curves can be
183 * plotted in the same picture with different plotting characters.
184 * Output of first 10 overprinted characters addad by FB 88
185 * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
186 *
187 * Input Variables:
188 * X = array containing the values of X
189 * Y = array containing the values of Y
190 * N = number of values in X and in Y
191 * can exceed the fixed number of lines
192 * M = number of different curves X,Y are containing
193 * MM = number of points in each curve i.e. N=M*MM
194 * XO = smallest value of X to be plotted
195 * DX = increment of X between subsequent lines
196 * YO = smallest value of Y to be plotted
197 * DY = increment of Y between subsequent character spaces
198 *
199 * other variables used inside:
200 * XX = numbers along the X-coordinate axis
201 * YY = numbers along the Y-coordinate axis
202 * LL = ten lines temporary storage for the plot
203 * L = character set used to plot different curves
204 * LOV = memorizes overprinted symbols
205 * the first 10 overprinted symbols are printed on
206 * the end of the line to avoid ambiguities
207 * (added by FB as considered quite helpful)
208 *
209 *********************************************************************
210 *
211  dimension xx(61),yy(61),ll(101,10)
212  dimension x(n),y(n),l(40),lov(40,10)
213  DATA l/
214  11h*,1h2,1h3,1h4,1h5,1h6,1h7,1h8,1h9,1hz,
215  21h+,1ha,1ho,1hb,1hc,1hd,1he,1hf,1hg,1hh,
216  31hi,1hj,1hk,1hl,1hm,1hn,1ho,1hp,1hq,1hr,
217  41hs,1ht,1hu,1hv,1hw,1hx,1hy,1h1,1h-,1h /
218 *
219 *
220  mn=51
221  amn=mn
222  DO 10 i=1,mn
223  ai=i-1
224  10 xx(i)=xo+ai*dx
225  DO 20 i=1,11
226  ai=i-1
227  20 yy(i)=yo+10.*ai*dy
228  WRITE(6, 500) (yy(i),i=1,11)
229  mmn=mn-1
230 *
231 *
232  DO 90 jj=1,mmn,10
233  jjj=jj-1
234  DO 30 i=1,101
235  DO 30 j=1,10
236  30 ll(i,j)=l(40)
237  DO 40 i=1,101
238  40 ll(i,1)=l(39)
239  DO 50 i=1,101,10
240  DO 50 j=1,10
241  50 ll(i,j)=l(38)
242  DO 60 i=1,40
243  DO 60 j=1,10
244  60 lov(i,j)=l(40)
245 *
246 *
247  DO 70 i=1,m
248  DO 70 j=1,mm
249  ii=j+(i-1)*mm
250  aix=(x(ii)-(xo-dx/2.))/dx+1.
251  aiy=(y(ii)-(yo-dy/2.))/dy+1.
252  aix=aix-float(jjj)
253 * changed Sept.88 by FB to avoid INTEGER OVERFLOW
254  IF( aix .GT. 1.d0.AND. aix .LT. 11.d0.AND. aiy .GT. 1.d0.and
255  + . aiy .LT. 102.d0) THEN
256  ix=aix
257  iy=aiy
258  IF( ix.GT. 0.AND. ix.LE. 10.AND. iy.GT. 0.AND. iy.LE. 101)
259  + THEN
260  IF(ll(iy,ix).NE.l(38).AND.ll(iy,ix).NE.l(39)) lov(i,ix)
261  + =ll(iy,ix)
262  ll(iy,ix)=l(i)
263  ENDIF
264  ENDIF
265  70 CONTINUE
266 *
267 *
268  DO 80 i=1,10
269  ii=i+jjj
270  iii=ii+1
271  WRITE(6,510) xx(ii),xx(iii) , (ll(j,i),j=1,101) , (lov(j,i),j
272  + =1,10)
273  80 CONTINUE
274  90 CONTINUE
275 *
276 *
277  WRITE(6, 520)
278  WRITE(6, 500) (yy(i),i=1,11)
279  RETURN
280 *
281  500 FORMAT(11x,11(1pe10.2),11hoverprinted)
282  510 FORMAT(1x,2(1pe10.2),101a1,1h ,10a1)
283  520 FORMAT(20x,10('1---------'),'1')
284  END
285 *-- Author :
286 C
287 C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
288 C
289  DOUBLE PRECISION FUNCTION dbetar(GAM,ETA)
290  IMPLICIT DOUBLE PRECISION (a-h,o-z)
291 C
292 C********************************************************************
293 C
294 C RANDOM NUMBER GENERATION FROM BETA
295 C DISTRIBUTION IN REGION 0.LE.X.LE.1.
296 C F(X) = X**(GAM-1.)*(1.-X)**(ETA-1)*GAMM(ETA+GAM) / (GAMM(GAM
297 C *GAMM(ETA))
298 C
299 C********************************************************************
300 C
301  y = dgamrn(1.d0,gam)
302  z = dgamrn(1.d0,eta)
303  dbetar = y/(y+z)
304  RETURN
305  END
306 *-- Author :
307  DOUBLE PRECISION FUNCTION dgamrn(ALAM,ETA)
308  IMPLICIT DOUBLE PRECISION (a-h,o-z)
309 C
310 C********************************************************************
311 C
312 C RANDOM NUMBER SELECTION FROM GAMMA DISTRIBUTION
313 C F(X) = ALAM**ETA*X**(ETA-1)*EXP(-ALAM*X) / GAM(ETA)
314 C
315 C********************************************************************
316 C
317  ncou=0
318  n = eta
319  f = eta - n
320  IF(f.EQ.0.d0) go to 20
321  10 r = rndm(r)
322  ncou=ncou+1
323  IF (ncou.GE.11) go to 20
324  IF(r.LT.f/(f+2.71828d0)) go to 30
325  yyy=log(rndm(y)+1.0d-9)/f
326  IF(abs(yyy).GT.50.0d0) go to 20
327  y = exp(yyy)
328  IF(log(rndm(r)+1.0d-9).GT.-y) go to 10
329  go to 40
330  20 y = 0
331  go to 50
332  30 y = 1.-log(rndm(y)+1.0d-9)
333  IF(rndm(r).GT.y**(f-1.d0)) go to 10
334  40 IF(n.EQ.0) go to 70
335  50 z = 1
336  DO 60 i = 1,n
337  60 z = z*rndm(z)
338  y = y-log(z+1.0d-9)
339  70 dgamrn = y/alam
340  RETURN
341  END
342 *-- Author :
343 C*****************************************************************
344  DOUBLE PRECISION FUNCTION betrej(GAM,ETA,XMIN,XMAX)
345  IMPLICIT DOUBLE PRECISION (a-h,o-z)
346  IF (xmin.GE.xmax)THEN
347  WRITE (6,500)xmin,xmax
348  stop
349  ENDIF
350  10 CONTINUE
351  xx=xmin+(xmax-xmin)*rndm(v)
352  betmax=xmin**(gam-1.)*(1.-xmin)**(eta-1.)
353  yy=betmax*rndm(w)
354  betxx=xx**(gam-1.)*(1.-xx)**(eta-1.)
355  IF(yy.GT.betxx) go to 10
356  betrej=xx
357  RETURN
358  500 FORMAT(' XMIN<XMAX IN BETREJ STOP ',2f10.5)
359  END
360 *-- Author :
361 C
362 C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
363 C
364  FUNCTION mcihad(MCIND)
365  IMPLICIT DOUBLE PRECISION (a-h,o-z)
366 C*** CALCULATION OF THE PARTICLE INDEX ACCORDING TO THE PDG PROPOSAL.
367 C MCIHAD PARTICLE NUMBER AS IN DECAY, BAMJET, HADEVT, FLUKA ETC.
368 C MCIND PDG PARTICLE NUMBER
369  INTEGER hamcin
370  COMMON /hamcin/ iamcin(410)
371  ih=0
372  mcihad=0
373  IF((mcind.EQ.0).OR.(mcind.GT.70000))RETURN
374  DO 10 i=1,410
375  ih=i
376  IF (iamcin(i).EQ.mcind) go to 20
377  10 CONTINUE
378  20 CONTINUE
379  mcihad=ih
380  RETURN
381  END
382 *-- Author :
383 C
384 C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
385 C
386  FUNCTION mpdgha(MCIND)
387  IMPLICIT DOUBLE PRECISION (a-h,o-z)
388 C*** CALCULATION OF THE PARTICLE INDEX ACCORDING TO THE PDG PROPOSAL.
389 C MCIND PARTICLE NUMBER AS IN DECAY, BAMJET, HADEVT, FLUKA ETC.
390 C MPDGHA PDG PARTICLE NUMBER
391  INTEGER hamcin
392  COMMON /hamcin/ iamcin(410)
393  mpdgha=iamcin(mcind)
394  RETURN
395  END
396 *-- Author :
397 C---------------------------------
398  BLOCK DATA radini
399  IMPLICIT DOUBLE PRECISION (a-h,o-z)
400  INTEGER hamcin
401  COMMON /rhamcin/ hamcin(200)
402  DATA hamcin/
403  *2212,-2212,-11,11,12, -12,22,2112,-2112,-13,
404  *13,130,211,-211,321, -321,3122,-3122,310,3112,
405  *3222,3212,111,311,-311, 0,0,0,0,0,
406  *221,213,113,-213,223, 323,313,-323,-313,99999,
407  *99999,99999,99999,30323,30313, -30323,-30313,3224,3214,3114,
408  *3216,3218,2224,2214,2114, 1114,12224,12214,12114,11114,
409  *12212,12112,22212,22112,99999,99999,-2224,-2214,-2114,-1114,
410  *10*99999,
411  *10*99999,
412  *4*99999,331,333,3322,3312,-3222,-3212,
413  *-3112,-3322,-3312,3224,3214,3114,3324,3314,3334,-3114,
414  *-3214,-3224,-3324,-3314,-3334,421,411,-411,-421,431,
415  *-431,441,8*99999,
416  *6*99999,4122,4232,4132,4222,
417  *4212,4112,6*99999,-4122,-4232,
418  *-4132,-4222,-4212,-4112,6*99999,
419  *40*99999/
420  END
421 
422 
423 C**********************************************************************
424  BLOCK DATA hadini
425 C**********************************************************************
426 C
427 C conversion table BAMJET --> PDG
428 C
429 C**********************************************************************
430  IMPLICIT DOUBLE PRECISION (a-h,o-z)
431  SAVE
432 C translation table version filled up by r.e. 25.01.94
433  COMMON /hamcin/ iamcin(410)
434 C
435  DATA iamcin /
436  &2212,-2212,11,-11,12, -12,22,2112,-2112,-13,
437  &13,130,211,-211,321, -321,3122,-3122,310,3112,
438  &3222,3212,111,311,-311, 0,0,0,0,0,
439  &221,213,113,-213,223, 323,313,-323,-313,10323,
440  &10313,-10323,-10313,30323,30313, -30323,-30313,3224,3214,3114,
441  &3216,3218,2224,2214,2114, 1114,12224,12214,12114,11114,
442  &12212,12112,22212,22112,32124, 31214,-2224,-2214,-2114,-1114,
443  &-12224,-12214,-12114,-11114,-2124, -1214,4*99999,
444  &5*99999, 5*99999,
445  &4*99999,331, 333,3322,3312,-3222,-3212,
446  &-3112,-3322,-3312,3224,3214, 3114,3324,3314,3334,-3224,
447  &-3214,-3114,-3324,-3314,-3334, 421,411,-411,-421,431,
448  &-431,441,423,413,-413, -423,433,-433,20443,443,
449  &-15,15,16,-16,14, -14,4122,4232,4132,4222,
450  & 4212,4112,4322,4312,4332, 4422,4412,4432,-4122,-4232,
451  & -4132,-4222,-4212,-4112,-4322, -4312,-4332,-4422,-4412,-4432,
452  & 4224,4214,4114,4324,4314, 4334,4424,4414,4434,4444,
453  & -4224,-4214,-4114,-4324,-4314, -4334,-4424,-4414,-4434,-4444,
454  &5*99999 , 20211,20111,-20211,99999,20321,
455  &-20321,20311,-20311,7*99999 ,
456  &7*99999,12212,12112,99999,
457  &115 ,215 ,-215 ,225 ,315 ,-315 ,325 ,-325 ,335 ,415 ,
458  &-415 ,425 ,-425 ,435 ,-435 ,445 ,511 ,-511 ,513 ,-513 ,
459  &515 ,-515 ,521 ,-521 ,523 ,-523 ,525 ,-525 ,531 ,-531 ,
460  &533 ,-533 ,535 ,-535 ,551 ,553 ,555 ,661 ,663 ,665 ,
461  &5112,-5112,5114,-5114,5122,-5122,5132,-5132,5212,-5212,
462  &5214 ,-5214,5222,-5222,5224,-5224,5232,-5232,5312,-5312,
463  &5314 ,-5314,5322,-5322,5324,-5324,5332,-5332,5334,-5334,
464  &10111 ,10113,10211,-10211,10213,-10213,10221,10223,
465  &10311,-10311,
466  & 10321,-10321,10331,10333,10411,-10411,10413,-10413,
467  &10421,-10421,
468  & 10423,-10423,10431,-10431,10433,10433,10411,10443,
469  &10511,-10511,
470  & 10513,-10513,10521,-10521,10523,-10523,10531,-10531,
471  &10533,-10533,
472  & 10551,10553,10661,10663,20113,20213,-20213,20223,
473  &20313,-20313,
474  & 20323,-20323,20333,20413,-20413,20423,-20423,20433,
475  &-20433,20443,
476  & 20513,-20513,20523,-20523,20533,-20533,20553,20663,
477  *99999,99999,
478  &40*99999 ,
479  &1,-1,2,-2,3,-3,4,-4,5,-5,
480  &6,-6,21,23,24,-24,25,99999,99999,99999/
481 C
482  END
483 C
484 *-- Author :
485 C
486 C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
487 C
488  FUNCTION mchad(ITDTU)
489  IMPLICIT DOUBLE PRECISION (a-h,o-z)
490 C*** CALCULATION OF THE PARTICLE INDEX TO BE USED FOR HADRIN
491 C IN SECONDARY COLLISIONS
492 C ITDTU PARTICLE NUMBER AS IN DECAY, BAMJET, HADEVT, FLUKA ETC.
493 C MCHAD PARTICLE NUMBER TO BE USED IN HADRIN
494  dimension itrans(210)
495  DATA itrans / 1, 2, -1, -1, -1, -1, -1, 8, 9, -1, -1, 24, 13, 14,
496  +15, 16, 8, 9, 25, 8, 1, 8, 23, 24, 25, -1, -1, -1, -1, -1, 23, 13,
497  +23, 14, 23, 15, 24, 16, 25, 15, 24, 16, 25, 15, 24, 16, 25, 1, 8,
498  +8, 8, 1, 1, 1, 8, 8, 1, 1, 8, 8, 1, 8, 1, 8, 1, 8, 2, 2, 9, 9, 2,
499  +2, 9, 9, 2, 9, 1, 13, 23, 14, 1, 1, 8, 8, 1, 1, 23, 14, 1, 8, 1,
500  +8, 1, 8, 23, 23, 8, 8, 2, 9, 9, 9, 9, 1, 8, 8, 8, 8, 8, 2, 9, 9,
501  +9, 9, 9, 85*- 1,7*-1,1,8,-1/
502 C
503  mchad=itrans(itdtu)
504  RETURN
505  END
506 *-- Author :
507 C
508 C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
509 C
510  SUBROUTINE dtrans(XO,YO,ZO,CDE,SDE,SFE,CFE,X,Y,Z)
511  IMPLICIT DOUBLE PRECISION (a-h,o-z)
512 C NEW TRANS VERSION OCTOBER 1987 J.RANFT
513 C ROTATION OF COORDINATE FRAME (1) DE RORATION AROUND Y AXIS
514 C (2) FE ROTATION AROUND Z AXIS
515  x= cde*cfe*xo-sfe*yo+sde*cfe*zo
516  y= cde*sfe*xo+cfe*yo+sde*sfe*zo
517  z=-sde *xo +cde *zo
518  RETURN
519  END
520 *-- Author :
521 C
522 C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
523 C
524 C SUBROUTINE DALTRA(GA,BGX,BGY,BGZ,PCX,PCY,PCZ,EC,P,PX,PY,PZ,E)
525 C IMPLICIT DOUBLE PRECISION (A-H,O-Z)
526 C*** ARBITRARY LORENTZ T\RANSFORM
527 C HUGE=1.D15
528 C EP=PCX*BGX+PCY*BGY+PCZ*BGZ
529 C PE=EP/(GA+1.)+EC
530 C PX=PC X+BG X*PE
531 C PY=PC Y+BG Y*PE
532 C PZ=PC Z+BG Z*PE
533 C PX=MIN(HUGE,MAX(-HUGE,P X))
534 C PY=MIN(HUGE,MAX(-HUGE,P Y))
535 C PZ=MIN(HUGE,MAX(-HUGE,P Z))
536 C P=SQRT(PX*PX+PY*PY+PZ*PZ)
537 C E=GA*EC+EP
538 C RETURN
539 C END
540 *===daltra=============================================================*
541 *
542  SUBROUTINE daltra(GA,BGX,BGY,BGZ,PCX,PCY,PCZ,EC,P,PX,PY,PZ,E)
543 
544 ************************************************************************
545 * Arbitrary Lorentz-transformation. *
546 * Adopted from the original by S. Roesler. This version dated 15.01.95 *
547 ************************************************************************
548 
549  IMPLICIT DOUBLE PRECISION (a-h,o-z)
550  parameter(huge=1.0d50,one=1.0d0)
551 
552  ep = pcx*bgx+pcy*bgy+pcz*bgz
553  pe = ep/(ga+one)+ec
554  px = pcx+bgx*pe
555  py = pcy+bgy*pe
556  pz = pcz+bgz*pe
557  px = min(huge,max(-huge,px))
558  py = min(huge,max(-huge,py))
559  pz = min(huge,max(-huge,pz))
560  p = sqrt(px*px+py*py+pz*pz)
561  e = ga*ec+ep
562 
563  RETURN
564  END
565 *-- Author :
566 C
567 C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
568 C
569  SUBROUTINE faltra(GA,BGA,CX,CY,CZ,COD,COF,SIF,PC,EC,P,PX,PY,PZ,E)
570  IMPLICIT DOUBLE PRECISION (a-h,o-z)
571  bgx=bga*cx
572  bgy=bga*cy
573  bgz=bga*cz
574  cod2=cod*cod
575  IF (cod2.GT.0.999999d0) cod2=0.999999d0
576  sid=sqrt(1.-cod2)*pc
577  pcx=sid*cof
578  pcy=sid*sif
579  pcz=cod*pc
580  ep=pcx*bgx+pcy*bgy+pcz*bgz
581  pe=ep/(ga+1.)+ec
582  px=pcx+bgx*pe
583  py=pcy+bgy*pe
584  pz=pcz+bgz*pe
585  p=sqrt(px*px+py*py+pz*pz)
586  pm=1./p
587  px=px*pm
588  py=py*pm
589  pz=pz*pm
590  e=ga*ec+ep
591  RETURN
592  END
593 *-- Author :
594 C
595 C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
596 C
597  SUBROUTINE dpoli(CS,SI)
598  IMPLICIT DOUBLE PRECISION (a-h,o-z)
599  u=rndm(v)
600  cs=rndm(vv)
601  IF (u.LT.0.5d0) cs=-cs
602  si=sqrt(1.-cs*cs+1.d-10)
603  RETURN
604  END
605 *-- Author :
606 C
607 C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
608 C
609  SUBROUTINE dfermi(GPART)
610  IMPLICIT DOUBLE PRECISION (a-h,o-z)
611  dimension g(3)
612  DO 10 i=1,3
613  g(i)=rndm(g(i))
614  10 CONTINUE
615 C
616 C FIND LARGEST OF 3 RANDOM NUMBERS
617 C
618  IF (g(3).LT.g(2)) go to 40
619  IF (g(3).LT.g(1)) go to 30
620  gpart=g(3)
621  20 RETURN
622  30 gpart=g(1)
623  go to 20
624  40 IF (g(2).LT.g(1)) go to 30
625  gpart=g(2)
626  go to 20
627  END
628 *-- Author :
629 C
630 C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
631 C--------------------------------------------------------------
632 C
633 C GQUAD.FOR
634 C
635 C----------------------------------------------------------------
636  DOUBLE PRECISION FUNCTION gquad(F,AX,BX,NX)
637  IMPLICIT DOUBLE PRECISION (a-h,o-z)
638 C N-POINT GAUSSIAN QUADRATURE OF FUNCTION F OVER INTERVAL (AX,BX).
639 C
640  EXTERNAL f
641  COMMON /gqcom/a(273),x(273),ktab(96)
642 C
643  CALL d106bd
644 C-----TEST N
645  n=nx
646  alpha=0.5*(bx+ax)
647  beta=0.5*(bx-ax)
648  IF(n.LT.1) go to 50
649  IF(n.NE.1) go to 10
650  gquad=(bx-ax)*f(alpha)
651  RETURN
652 C
653  10 IF(n.LE.16) go to 20
654  IF(n.EQ.20) go to 20
655  IF(n.EQ.24) go to 20
656  IF(n.EQ.32) go to 20
657  IF(n.EQ.40) go to 20
658  IF(n.EQ.48) go to 20
659  IF(n.EQ.64) go to 20
660  IF(n.EQ.80) go to 20
661  IF(n.EQ.96) go to 20
662  go to 50
663 C
664 C----- SET K EQUAL TO INITIAL SUBSCRIPT AND INTEGRATE
665  20 k=ktab(n)
666  m=n/2
667  sum=0.0
668  jmax=k-1+m
669 C
670  DO 30 j=k,jmax
671  delta=beta*x(j)
672  sum=sum+a(j)*(f(alpha+delta)+f(alpha-delta))
673  30 CONTINUE
674 C
675  IF(n-m-m.EQ.0) go to 40
676  jmid=k+m
677  sum=sum+a(jmid)*f(alpha)
678  40 gquad=beta*sum
679  RETURN
680 C
681  50 gquad=0.0
682  zn=n
683  WRITE(6, 500)zn
684  RETURN
685 C
686  500 FORMAT( 42h gquad ... n has the non-permissible value ,e11. 3)
687  END
688 *-- Author :
689  SUBROUTINE gset(AX,BX,NX,Z,W)
690  IMPLICIT DOUBLE PRECISION (a-h,o-z)
691 C N-POINT GAUSS ZEROS AND WEIGHTS FOR THE INTERVAL (AX,BX) ARE
692 C STORED IN ARRAYS Z AND W RESPECTIVELY.
693 C
694  COMMON /gqcom/a(273),x(273),ktab(96)
695  dimension z(192),w(192)
696 C
697  CALL d106bd
698 C-----TEST N
699  n=nx
700  alpha=0.5*(bx+ax)
701  beta=0.5*(bx-ax)
702  IF(n.LT.1) go to 40
703  IF(n.NE.1) go to 10
704  z(1)=alpha
705  w(1)=bx-ax
706  RETURN
707 C
708  10 IF(n.LE.16) go to 20
709  IF(n.EQ.20) go to 20
710  IF(n.EQ.24) go to 20
711  IF(n.EQ.32) go to 20
712  IF(n.EQ.40) go to 20
713  IF(n.EQ.48) go to 20
714  IF(n.EQ.64) go to 20
715  IF(n.EQ.80) go to 20
716  IF(n.EQ.96) go to 20
717  go to 40
718 C
719 C----- SET K EQUAL TO INITIAL SUBSCRIPT AND STORE RESULTS
720  20 k=ktab(n)
721  m=n/2
722 C
723  DO 30 j=1,m
724  jtab=k-1+j
725  wtemp=beta*a(jtab)
726  delta=beta*x(jtab)
727  z(j)=alpha-delta
728  w(j)=wtemp
729  jp=n+1-j
730  z(jp)=alpha+delta
731  w(jp)=wtemp
732  30 CONTINUE
733 C
734  IF((n-m-m).EQ.0) RETURN
735  z(m+1)=alpha
736  jmid=k+m
737  w(m+1)=beta*a(jmid)
738  RETURN
739 C
740  40 zn=n
741  WRITE(6, 500)zn
742  RETURN
743 C
744  500 FORMAT( 41h gset ... n has the non-permissible value ,e11. 3)
745  END
746 *-- Author :
747  SUBROUTINE d106bd
748  IMPLICIT DOUBLE PRECISION (a-h,o-z)
749  COMMON /gqcom/ b(273),y(273),ltab(96)
750  dimension a(273),x(273),ktab(96)
751 C
752 C-----TABLE OF INITIAL SUBSCRIPTS FOR N=2(1)16(4)96
753  DATA ktab(2)/1/
754  DATA ktab(3)/2/
755  DATA ktab(4)/4/
756  DATA ktab(5)/6/
757  DATA ktab(6)/9/
758  DATA ktab(7)/12/
759  DATA ktab(8)/16/
760  DATA ktab(9)/20/
761  DATA ktab(10)/25/
762  DATA ktab(11)/30/
763  DATA ktab(12)/36/
764  DATA ktab(13)/42/
765  DATA ktab(14)/49/
766  DATA ktab(15)/56/
767  DATA ktab(16)/64/
768  DATA ktab(20)/72/
769  DATA ktab(24)/82/
770  DATA ktab(28)/82/
771  DATA ktab(32)/94/
772  DATA ktab(36)/94/
773  DATA ktab(40)/110/
774  DATA ktab(44)/110/
775  DATA ktab(48)/130/
776  DATA ktab(52)/130/
777  DATA ktab(56)/130/
778  DATA ktab(60)/130/
779  DATA ktab(64)/154/
780  DATA ktab(68)/154/
781  DATA ktab(72)/154/
782  DATA ktab(76)/154/
783  DATA ktab(80)/186/
784  DATA ktab(84)/186/
785  DATA ktab(88)/186/
786  DATA ktab(92)/186/
787  DATA ktab(96)/226/
788 C
789 C-----TABLE OF ABSCISSAE (X) AND WEIGHTS (A) FOR INTERVAL (-1,+1).
790 C
791 C-----N=2
792  DATA x(1)/0.577350269189626 /, a(1)/1.000000000000000 /
793 C-----N=3
794  DATA x(2)/0.774596669241483 /, a(2)/0.555555555555556 /
795  DATA x(3)/0.000000000000000 /, a(3)/0.888888888888889 /
796 C-----N=4
797  DATA x(4)/0.861136311594053 /, a(4)/0.347854845137454 /
798  DATA x(5)/0.339981043584856 /, a(5)/0.652145154862546 /
799 C-----N=5
800  DATA x(6)/0.906179845938664 /, a(6)/0.236926885056189 /
801  DATA x(7)/0.538469310105683 /, a(7)/0.478628670499366 /
802  DATA x(8)/0.000000000000000 /, a(8)/0.568888888888889 /
803 C-----N=6
804  DATA x(9)/0.932469514203152 /, a(9)/0.171324492379170 /
805  DATA x(10)/0.661209386466265 /, a(10)/0.360761573048139 /
806  DATA x(11)/0.238619186083197 /, a(11)/0.467913934572691 /
807 C-----N=7
808  DATA x(12)/0.949107912342759 /, a(12)/0.129484966168870 /
809  DATA x(13)/0.741531185599394 /, a(13)/0.279705391489277 /
810  DATA x(14)/0.405845151377397 /, a(14)/0.381830050505119 /
811  DATA x(15)/0.000000000000000 /, a(15)/0.417959183673469 /
812 C-----N=8
813  DATA x(16)/0.960289856497536 /, a(16)/0.101228536290376 /
814  DATA x(17)/0.796666477413627 /, a(17)/0.222381034453374 /
815  DATA x(18)/0.525532409916329 /, a(18)/0.313706645877887 /
816  DATA x(19)/0.183434642495650 /, a(19)/0.362683783378362 /
817 C-----N=9
818  DATA x(20)/0.968160239507626 /, a(20)/0.081274388361574 /
819  DATA x(21)/0.836031107326636 /, a(21)/0.180648160694857 /
820  DATA x(22)/0.613371432700590 /, a(22)/0.260610696402935 /
821  DATA x(23)/0.324253423403809 /, a(23)/0.312347077040003 /
822  DATA x(24)/0.000000000000000 /, a(24)/0.330239355001260 /
823 C-----N=10
824  DATA x(25)/0.973906528517172 /, a(25)/0.066671344308688 /
825  DATA x(26)/0.865063366688985 /, a(26)/0.149451349150581 /
826  DATA x(27)/0.679409568299024 /, a(27)/0.219086362515982 /
827  DATA x(28)/0.433395394129247 /, a(28)/0.269266719309996 /
828  DATA x(29)/0.148874338981631 /, a(29)/0.295524224714753 /
829 C-----N=11
830  DATA x(30)/0.978228658146057 /, a(30)/0.055668567116174 /
831  DATA x(31)/0.887062599768095 /, a(31)/0.125580369464905 /
832  DATA x(32)/0.730152005574049 /, a(32)/0.186290210927734 /
833  DATA x(33)/0.519096129206812 /, a(33)/0.233193764591990 /
834  DATA x(34)/0.269543155952345 /, a(34)/0.262804544510247 /
835  DATA x(35)/0.000000000000000 /, a(35)/0.272925086777901 /
836 C-----N=12
837  DATA x(36)/0.981560634246719 /, a(36)/0.047175336386512 /
838  DATA x(37)/0.904117256370475 /, a(37)/0.106939325995318 /
839  DATA x(38)/0.769902674194305 /, a(38)/0.160078328543346 /
840  DATA x(39)/0.587317954286617 /, a(39)/0.203167426723066 /
841  DATA x(40)/0.367831498998180 /, a(40)/0.233492536538355 /
842  DATA x(41)/0.125233408511469 /, a(41)/0.249147045813403 /
843 C-----N=13
844  DATA x(42)/0.984183054718588 /, a(42)/0.040484004765316 /
845  DATA x(43)/0.917598399222978 /, a(43)/0.092121499837728 /
846  DATA x(44)/0.801578090733310 /, a(44)/0.138873510219787 /
847  DATA x(45)/0.642349339440340 /, a(45)/0.178145980761946 /
848  DATA x(46)/0.448492751036447 /, a(46)/0.207816047536889 /
849  DATA x(47)/0.230458315955135 /, a(47)/0.226283180262897 /
850  DATA x(48)/0.000000000000000 /, a(48)/0.232551553230874 /
851 C-----N=14
852  DATA x(49)/0.986283808696812 /, a(49)/0.035119460331752 /
853  DATA x(50)/0.928434883663574 /, a(50)/0.080158087159760 /
854  DATA x(51)/0.827201315069765 /, a(51)/0.121518570687903 /
855  DATA x(52)/0.687292904811685 /, a(52)/0.157203167158194 /
856  DATA x(53)/0.515248636358154 /, a(53)/0.185538397477938 /
857  DATA x(54)/0.319112368927890 /, a(54)/0.205198463721296 /
858  DATA x(55)/0.108054948707344 /, a(55)/0.215263853463158 /
859 C-----N=15
860  DATA x(56)/0.987992518020485 /, a(56)/0.030753241996117 /
861  DATA x(57)/0.937273392400706 /, a(57)/0.070366047488108 /
862  DATA x(58)/0.848206583410427 /, a(58)/0.107159220467172 /
863  DATA x(59)/0.724417731360170 /, a(59)/0.139570677926154 /
864  DATA x(60)/0.570972172608539 /, a(60)/0.166269205816994 /
865  DATA x(61)/0.394151347077563 /, a(61)/0.186161000015562 /
866  DATA x(62)/0.201194093997435 /, a(62)/0.198431485327111 /
867  DATA x(63)/0.000000000000000 /, a(63)/0.202578241925561 /
868 C-----N=16
869  DATA x(64)/0.989400934991650 /, a(64)/0.027152459411754 /
870  DATA x(65)/0.944575023073233 /, a(65)/0.062253523938648 /
871  DATA x(66)/0.865631202387832 /, a(66)/0.095158511682493 /
872  DATA x(67)/0.755404408355003 /, a(67)/0.124628971255534 /
873  DATA x(68)/0.617876244402644 /, a(68)/0.149595988816577 /
874  DATA x(69)/0.458016777657227 /, a(69)/0.169156519395003 /
875  DATA x(70)/0.281603550779259 /, a(70)/0.182603415044924 /
876  DATA x(71)/0.095012509837637 /, a(71)/0.189450610455069 /
877 C-----N=20
878  DATA x(72)/0.993128599185094 /, a(72)/0.017614007139152 /
879  DATA x(73)/0.963971927277913 /, a(73)/0.040601429800386 /
880  DATA x(74)/0.912234428251325 /, a(74)/0.062672048334109 /
881  DATA x(75)/0.839116971822218 /, a(75)/0.083276741576704 /
882  DATA x(76)/0.746331906460150 /, a(76)/0.101930119817240 /
883  DATA x(77)/0.636053680726515 /, a(77)/0.118194531961518 /
884  DATA x(78)/0.510867001950827 /, a(78)/0.131688638449176 /
885  DATA x(79)/0.373706088715419 /, a(79)/0.142096109318382 /
886  DATA x(80)/0.227785851141645 /, a(80)/0.149172986472603 /
887  DATA x(81)/0.076526521133497 /, a(81)/0.152753387130725 /
888 C-----N=24
889  DATA x(82)/0.995187219997021 /, a(82)/0.012341229799987 /
890  DATA x(83)/0.974728555971309 /, a(83)/0.028531388628933 /
891  DATA x(84)/0.938274552002732 /, a(84)/0.044277438817419 /
892  DATA x(85)/0.886415527004401 /, a(85)/0.059298584915436 /
893  DATA x(86)/0.820001985973902 /, a(86)/0.073346481411080 /
894  DATA x(87)/0.740124191578554 /, a(87)/0.086190161531953 /
895  DATA x(88)/0.648093651936975 /, a(88)/0.097618652104113 /
896  DATA x(89)/0.545421471388839 /, a(89)/0.107444270115965 /
897  DATA x(90)/0.433793507626045 /, a(90)/0.115505668053725 /
898  DATA x(91)/0.315042679696163 /, a(91)/0.121670472927803 /
899  DATA x(92)/0.191118867473616 /, a(92)/0.125837456346828 /
900  DATA x(93)/0.064056892862605 /, a(93)/0.127938195346752 /
901 C-----N=32
902  DATA x(94)/0.997263861849481 /, a(94)/0.007018610009470 /
903  DATA x(95)/0.985611511545268 /, a(95)/0.016274394730905 /
904  DATA x(96)/0.964762255587506 /, a(96)/0.025392065309262 /
905  DATA x(97)/0.934906075937739 /, a(97)/0.034273862913021 /
906  DATA x(98)/0.896321155766052 /, a(98)/0.042835898022226 /
907  DATA x(99)/0.849367613732569 /, a(99)/0.050998059262376 /
908  DATA x(100)/0.794483795967942/, a(100)/0.058684093478535/
909  DATA x(101)/0.732182118740289/, a(101)/0.065822222776361/
910  DATA x(102)/0.663044266930215/, a(102)/0.072345794108848/
911  DATA x(103)/0.587715757240762/, a(103)/0.078193895787070/
912  DATA x(104)/0.506899908932229/, a(104)/0.083311924226946/
913  DATA x(105)/0.421351276130635/, a(105)/0.087652093004403/
914  DATA x(106)/0.331868602282127/, a(106)/0.091173878695763/
915  DATA x(107)/0.239287362252137/, a(107)/0.093844399080804/
916  DATA x(108)/0.144471961582796/, a(108)/0.095638720079274/
917  DATA x(109)/0.048307665687738/, a(109)/0.096540088514727/
918 C-----N=40
919  DATA x(110)/0.998237709710559/, a(110)/0.004521277098533/
920  DATA x(111)/0.990726238699457/, a(111)/0.010498284531152/
921  DATA x(112)/0.977259949983774/, a(112)/0.016421058381907/
922  DATA x(113)/0.957916819213791/, a(113)/0.022245849194166/
923  DATA x(114)/0.932812808278676/, a(114)/0.027937006980023/
924  DATA x(115)/0.902098806968874/, a(115)/0.033460195282547/
925  DATA x(116)/0.865959503212259/, a(116)/0.038782167974472/
926  DATA x(117)/0.824612230833311/, a(117)/0.043870908185673/
927  DATA x(118)/0.778305651426519/, a(118)/0.048695807635072/
928  DATA x(119)/0.727318255189927/, a(119)/0.053227846983936/
929  DATA x(120)/0.671956684614179/, a(120)/0.057439769099391/
930  DATA x(121)/0.612553889667980/, a(121)/0.061306242492928/
931  DATA x(122)/0.549467125095128/, a(122)/0.064804013456601/
932  DATA x(123)/0.483075801686178/, a(123)/0.067912045815233/
933  DATA x(124)/0.413779204371605/, a(124)/0.070611647391286/
934  DATA x(125)/0.341994090825758/, a(125)/0.072886582395804/
935  DATA x(126)/0.268152185007253/, a(126)/0.074723169057968/
936  DATA x(127)/0.192697580701371/, a(127)/0.076110361900626/
937  DATA x(128)/0.116084070675255/, a(128)/0.077039818164247/
938  DATA x(129)/0.038772417506050/, a(129)/0.077505947978424/
939 C-----N=48
940  DATA x(130)/0.998771007252426/, a(130)/0.003153346052305/
941  DATA x(131)/0.993530172266350/, a(131)/0.007327553901276/
942  DATA x(132)/0.984124583722826/, a(132)/0.011477234579234/
943  DATA x(133)/0.970591592546247/, a(133)/0.015579315722943/
944  DATA x(134)/0.952987703160430/, a(134)/0.019616160457355/
945  DATA x(135)/0.931386690706554/, a(135)/0.023570760839324/
946  DATA x(136)/0.905879136715569/, a(136)/0.027426509708356/
947  DATA x(137)/0.876572020274247/, a(137)/0.031167227832798/
948  DATA x(138)/0.843588261624393/, a(138)/0.034777222564770/
949  DATA x(139)/0.807066204029442/, a(139)/0.038241351065830/
950  DATA x(140)/0.767159032515740/, a(140)/0.041545082943464/
951  DATA x(141)/0.724034130923814/, a(141)/0.044674560856694/
952  DATA x(142)/0.677872379632663/, a(142)/0.047616658492490/
953  DATA x(143)/0.628867396776513/, a(143)/0.050359035553854/
954  DATA x(144)/0.577224726083972/, a(144)/0.052890189485193/
955  DATA x(145)/0.523160974722233/, a(145)/0.055199503699984/
956  DATA x(146)/0.466902904750958/, a(146)/0.057277292100403/
957  DATA x(147)/0.408686481990716/, a(147)/0.059114839698395/
958  DATA x(148)/0.348755886292160/, a(148)/0.060704439165893/
959  DATA x(149)/0.287362487355455/, a(149)/0.062039423159892/
960  DATA x(150)/0.224763790394689/, a(150)/0.063114192286254/
961  DATA x(151)/0.161222356068891/, a(151)/0.063924238584648/
962  DATA x(152)/0.097004699209462/, a(152)/0.064466164435950/
963  DATA x(153)/0.032380170962869/, a(153)/0.064737696812683/
964 C-----N=64
965  DATA x(154)/0.999305041735772/, a(154)/0.001783280721696/
966  DATA x(155)/0.996340116771955/, a(155)/0.004147033260562/
967  DATA x(156)/0.991013371476744/, a(156)/0.006504457968978/
968  DATA x(157)/0.983336253884625/, a(157)/0.008846759826363/
969  DATA x(158)/0.973326827789910/, a(158)/0.011168139460131/
970  DATA x(159)/0.961008799652053/, a(159)/0.013463047896718/
971  DATA x(160)/0.946411374858402/, a(160)/0.015726030476024/
972  DATA x(161)/0.929569172131939/, a(161)/0.017951715775697/
973  DATA x(162)/0.910522137078502/, a(162)/0.020134823153530/
974  DATA x(163)/0.889315445995114/, a(163)/0.022270173808383/
975  DATA x(164)/0.865999398154092/, a(164)/0.024352702568710/
976  DATA x(165)/0.840629296252580/, a(165)/0.026377469715054/
977  DATA x(166)/0.813265315122797/, a(166)/0.028339672614259/
978  DATA x(167)/0.783972358943341/, a(167)/0.030234657072402/
979  DATA x(168)/0.752819907260531/, a(168)/0.032057928354851/
980  DATA x(169)/0.719881850171610/, a(169)/0.033805161837141/
981  DATA x(170)/0.685236313054233/, a(170)/0.035472213256882/
982  DATA x(171)/0.648965471254657/, a(171)/0.037055128540240/
983  DATA x(172)/0.611155355172393/, a(172)/0.038550153178615/
984  DATA x(173)/0.571895646202634/, a(173)/0.039953741132720/
985  DATA x(174)/0.531279464019894/, a(174)/0.041262563242623/
986  DATA x(175)/0.489403145707052/, a(175)/0.042473515123653/
987  DATA x(176)/0.446366017253464/, a(176)/0.043583724529323/
988  DATA x(177)/0.402270157963991/, a(177)/0.044590558163756/
989  DATA x(178)/0.357220158337668/, a(178)/0.045491627927418/
990  DATA x(179)/0.311322871990210/, a(179)/0.046284796581314/
991  DATA x(180)/0.264687162208767/, a(180)/0.046968182816210/
992  DATA x(181)/0.217423643740007/, a(181)/0.047540165714830/
993  DATA x(182)/0.169644420423992/, a(182)/0.047999388596458/
994  DATA x(183)/0.121462819296120/, a(183)/0.048344762234802/
995  DATA x(184)/0.072993121787799/, a(184)/0.048575467441503/
996  DATA x(185)/0.024350292663424/, a(185)/0.048690957009139/
997 C-----N=80
998  DATA x(186)/0.999553822651630/, a(186)/0.001144950003186/
999  DATA x(187)/0.997649864398237/, a(187)/0.002663533589512/
1000  DATA x(188)/0.994227540965688/, a(188)/0.004180313124694/
1001  DATA x(189)/0.989291302499755/, a(189)/0.005690922451403/
1002  DATA x(190)/0.982848572738629/, a(190)/0.007192904768117/
1003  DATA x(191)/0.974909140585727/, a(191)/0.008683945269260/
1004  DATA x(192)/0.965485089043799/, a(192)/0.010161766041103/
1005  DATA x(193)/0.954590766343634/, a(193)/0.011624114120797/
1006  DATA x(194)/0.942242761309872/, a(194)/0.013068761592401/
1007  DATA x(195)/0.928459877172445/, a(195)/0.014493508040509/
1008  DATA x(196)/0.913263102571757/, a(196)/0.015896183583725/
1009  DATA x(197)/0.896675579438770/, a(197)/0.017274652056269/
1010  DATA x(198)/0.878722567678213/, a(198)/0.018626814208299/
1011  DATA x(199)/0.859431406663111/, a(199)/0.019950610878141/
1012  DATA x(200)/0.838831473580255/, a(200)/0.021244026115782/
1013  DATA x(201)/0.816954138681463/, a(201)/0.022505090246332/
1014  DATA x(202)/0.793832717504605/, a(202)/0.023731882865930/
1015  DATA x(203)/0.769502420135041/, a(203)/0.024922535764115/
1016  DATA x(204)/0.744000297583597/, a(204)/0.026075235767565/
1017  DATA x(205)/0.717365185362099/, a(205)/0.027188227500486/
1018  DATA x(206)/0.689637644342027/, a(206)/0.028259816057276/
1019  DATA x(207)/0.660859898986119/, a(207)/0.029288369583267/
1020  DATA x(208)/0.631075773046871/, a(208)/0.030272321759557/
1021  DATA x(209)/0.600330622829751/, a(209)/0.031210174188114/
1022  DATA x(210)/0.568671268122709/, a(210)/0.032100498673487/
1023  DATA x(211)/0.536145920897131/, a(211)/0.032941939397645/
1024  DATA x(212)/0.502804111888784/, a(212)/0.033733214984611/
1025  DATA x(213)/0.468696615170544/, a(213)/0.034473120451753/
1026  DATA x(214)/0.433875370831756/, a(214)/0.035160529044747/
1027  DATA x(215)/0.398393405881969/, a(215)/0.035794393953416/
1028  DATA x(216)/0.362304753499487/, a(216)/0.036373749905835/
1029  DATA x(217)/0.325664370747701/, a(217)/0.036897714638276/
1030  DATA x(218)/0.288528054884511/, a(218)/0.037365490238730/
1031  DATA x(219)/0.250952358392272/, a(219)/0.037776364362001/
1032  DATA x(220)/0.212994502857666/, a(220)/0.038129711314477/
1033  DATA x(221)/0.174712291832646/, a(221)/0.038424993006959/
1034  DATA x(222)/0.136164022809143/, a(222)/0.038661759774076/
1035  DATA x(223)/0.097408398441584/, a(223)/0.038839651059051/
1036  DATA x(224)/0.058504437152420/, a(224)/0.038958395962769/
1037  DATA x(225)/0.019511383256793/, a(225)/0.039017813656306/
1038 C-----N=96
1039  DATA x(226)/0.999689503883230/, a(226)/0.000796792065552/
1040  DATA x(227)/0.998364375863181/, a(227)/0.001853960788946/
1041  DATA x(228)/0.995981842987209/, a(228)/0.002910731817934/
1042  DATA x(229)/0.992543900323762/, a(229)/0.003964554338444/
1043  DATA x(230)/0.988054126329623/, a(230)/0.005014202742927/
1044  DATA x(231)/0.982517263563014/, a(231)/0.006058545504235/
1045  DATA x(232)/0.975939174585136/, a(232)/0.007096470791153/
1046  DATA x(233)/0.968326828463264/, a(233)/0.008126876925698/
1047  DATA x(234)/0.959688291448742/, a(234)/0.009148671230783/
1048  DATA x(235)/0.950032717784437/, a(235)/0.010160770535008/
1049  DATA x(236)/0.939370339752755/, a(236)/0.011162102099838/
1050  DATA x(237)/0.927712456722308/, a(237)/0.012151604671088/
1051  DATA x(238)/0.915071423120898/, a(238)/0.013128229566961/
1052  DATA x(239)/0.901460635315852/, a(239)/0.014090941772314/
1053  DATA x(240)/0.886894517402420/, a(240)/0.015038721026994/
1054  DATA x(241)/0.871388505909296/, a(241)/0.015970562902562/
1055  DATA x(242)/0.854959033434601/, a(242)/0.016885479864245/
1056  DATA x(243)/0.837623511228187/, a(243)/0.017782502316045/
1057  DATA x(244)/0.819400310737931/, a(244)/0.018660679627411/
1058  DATA x(245)/0.800308744139140/, a(245)/0.019519081140145/
1059  DATA x(246)/0.780369043867433/, a(246)/0.020356797154333/
1060  DATA x(247)/0.759602341176647/, a(247)/0.021172939892191/
1061  DATA x(248)/0.738030643744400/, a(248)/0.021966644438744/
1062  DATA x(249)/0.715676812348967/, a(249)/0.022737069658329/
1063  DATA x(250)/0.692564536642171/, a(250)/0.023483399085926/
1064  DATA x(251)/0.668718310043916/, a(251)/0.024204841792364/
1065  DATA x(252)/0.644163403784967/, a(252)/0.024900633222483/
1066  DATA x(253)/0.618925840125468/, a(253)/0.025570036005349/
1067  DATA x(254)/0.593032364777572/, a(254)/0.026212340735672/
1068  DATA x(255)/0.566510418561397/, a(255)/0.026826866725591/
1069  DATA x(256)/0.539388108324357/, a(256)/0.027412962726029/
1070  DATA x(257)/0.511694177154667/, a(257)/0.027970007616848/
1071  DATA x(258)/0.483457973920596/, a(258)/0.028497411065085/
1072  DATA x(259)/0.454709422167743/, a(259)/0.028994614150555/
1073  DATA x(260)/0.425478988407300/, a(260)/0.029461089958167/
1074  DATA x(261)/0.395797649828908/, a(261)/0.029896344136328/
1075  DATA x(262)/0.365696861472313/, a(262)/0.030299915420827/
1076  DATA x(263)/0.335208522892625/, a(263)/0.030671376123669/
1077  DATA x(264)/0.304364944354496/, a(264)/0.031010332586313/
1078  DATA x(265)/0.273198812591049/, a(265)/0.031316425596861/
1079  DATA x(266)/0.241743156163840/, a(266)/0.031589330770727/
1080  DATA x(267)/0.210031310460567/, a(267)/0.031828758894411/
1081  DATA x(268)/0.178096882367618/, a(268)/0.032034456231992/
1082  DATA x(269)/0.145973714654896/, a(269)/0.032206204794030/
1083  DATA x(270)/0.113695850110665/, a(270)/0.032343822568575/
1084  DATA x(271)/0.081297495464425/, a(271)/0.032447163714064/
1085  DATA x(272)/0.048812985136049/, a(272)/0.032516118713868/
1086  DATA x(273)/0.016276744849602/, a(273)/0.032550614492363/
1087  DATA ibd/0/
1088  IF(ibd.NE.0) RETURN
1089  ibd=1
1090  DO 10 i=1,273
1091  b(i) = a(i)
1092  10 y(i) = x(i)
1093  DO 20 i=1,96
1094  20 ltab(i) = ktab(i)
1095  RETURN
1096  END
1097 *-- Author :
1098 C*********************************************************************
1099  DOUBLE PRECISION FUNCTION samsqx(X1,X2)
1100  IMPLICIT DOUBLE PRECISION (a-h,o-z)
1101 C SAMPLING FROM F(X)=1/X**0.5 BETWEEN X1 AND X2
1102  r=rndm(v)
1103  samsqx=(r*sqrt(x2)+(1.-r)*sqrt(x1))**2
1104  RETURN
1105  END
1106 C*********************************************************************
1107  DOUBLE PRECISION FUNCTION sampey(X1,X2)
1108  IMPLICIT DOUBLE PRECISION (a-h,o-z)
1109  DATA isampe/1/
1110 C SAMPLING FROM F(X)=1/X BETWEEN X1 AND X2
1111  r=rndm(v)
1112  al1=log(x1)
1113  al2=log(x2)
1114  sampey=exp((1.-r)*al1+r*al2)
1115  RETURN
1116  END
1117 C*********************************************************************
1118  DOUBLE PRECISION FUNCTION sampex(X1,X2)
1119  IMPLICIT DOUBLE PRECISION (a-h,o-z)
1120  DATA isampe/1/
1121 C SAMPLING FROM F(X)=1/X BETWEEN X1 AND X2
1122  IF(isampe.EQ.0)THEN
1123  r=rndm(v)
1124  al1=log(x1)
1125  al2=log(x2)
1126  sampex=exp((1.-r)*al1+r*al2)
1127  ELSEIF(isampe.EQ.1)THEN
1128  sampex=samsqx(x1,x2)
1129  ENDIF
1130  RETURN
1131  END
1132 C*********************************************************************
1133  DOUBLE PRECISION FUNCTION sampxb(X1,X2,B)
1134  IMPLICIT DOUBLE PRECISION (a-h,o-z)
1135 C Sampling x values between x1 and x2 from
1136 C f(x)=1./SQRT(X**2+B**2)
1137 C
1138  a1=log(x1+sqrt(x1**2+b**2))
1139  a2=log(x2+sqrt(x2**2+b**2))
1140  an=a2-a1
1141  a=an*rndm(v)+a1
1142  bb=exp(a)
1143  sampxb=(bb**2-b**2)/(2.*bb)
1144  RETURN
1145  END
1146 C
1147  SUBROUTINE sorti(A,N)
1148  IMPLICIT DOUBLE PRECISION (a-h,o-z)
1149  dimension a(n)
1150  m=n
1151  10 CONTINUE
1152  m=n-1
1153  IF(m.LE.0) RETURN
1154  l=0
1155  DO 20 i=1,m
1156  j=i+1
1157  IF (a(i).LE.a(j)) go to 20
1158  b=a(i)
1159  a(i)=a(j)
1160  a(j)=b
1161  l=1
1162  20 CONTINUE
1163  IF(l.EQ.1) go to 10
1164  RETURN
1165  END
1166 *$ CREATE RM48.FOR
1167 *COPY RM48
1168  SUBROUTINE rm48(RVEC,LENV)
1169 C Double-precision version of
1170 C Universal random number generator proposed by Marsaglia and Zaman
1171 C in report FSU-SCRI-87-50
1172 C based on RANMAR, modified by F. James, to generate vectors
1173 C of pseudorandom numbers RVEC of length LENV, where the numbers
1174 C in RVEC are numbers with at least 48-bit mantissas.
1175 C Input and output entry points: RM48IN, RM48UT.
1176 C!!! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1177 C!!! Calling sequences for RM48: ++
1178 C!!! CALL RM48 (RVEC, LEN) returns a vector RVEC of LEN ++
1179 C!!! 64-bit random floating point numbers between ++
1180 C!!! zero and one. ++
1181 C!!! CALL RM48IN(I1,N1,N2) initializes the generator from one ++
1182 C!!! 64-bit integer I1, and number counts N1,N2 ++
1183 C!!! (for initializing, set N1=N2=0, but to restart ++
1184 C!!! a previously generated sequence, use values ++
1185 C!!! output by RM48UT) ++
1186 C!!! CALL RM48UT(I1,N1,N2) outputs the value of the original ++
1187 C!!! seed and the two number counts, to be used ++
1188 C!!! for restarting by initializing to I1 and ++
1189 C!!! skipping N2*100000000+N1 numbers. ++
1190 C!!! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1191 C for 32-bit machines, use IMPLICIT DOUBLE PRECISION
1192 C INCLUDE '(DBLPRC)'
1193 *$ CREATE DBLPRC.ADD
1194  IMPLICIT DOUBLE PRECISION (a-h,o-z)
1195  parameter( kalgnm = 2 )
1196  parameter( anglgb = 5.0d-16 )
1197  parameter( anglsq = 2.5d-31 )
1198  parameter( axcssv = 0.2d+16 )
1199  parameter( andrfl = 1.0d-38 )
1200  parameter( avrflw = 1.0d+38 )
1201  parameter( ainfnt = 1.0d+30 )
1202  parameter( azrzrz = 1.0d-30 )
1203  parameter( einfnt = +69.07755278982137 d+00 )
1204  parameter( ezrzrz = -69.07755278982137 d+00 )
1205  parameter( onemns = 0.999999999999999 d+00 )
1206  parameter( onepls = 1.000000000000001 d+00 )
1207  parameter( csnnrm = 2.0d-15 )
1208  parameter( dmxtrn = 1.0d+08 )
1209  parameter( zerzer = 0.d+00 )
1210  parameter( oneone = 1.d+00 )
1211  parameter( twotwo = 2.d+00 )
1212  parameter( thrthr = 3.d+00 )
1213  parameter( foufou = 4.d+00 )
1214  parameter( fivfiv = 5.d+00 )
1215  parameter( sixsix = 6.d+00 )
1216  parameter( sevsev = 7.d+00 )
1217  parameter( eigeig = 8.d+00 )
1218  parameter( aninen = 9.d+00 )
1219  parameter( tenten = 10.d+00 )
1220  parameter( hlfhlf = 0.5d+00 )
1221  parameter( onethi = oneone / thrthr )
1222  parameter( twothi = twotwo / thrthr )
1223  parameter( onefou = oneone / foufou )
1224  parameter( thrtwo = thrthr / twotwo )
1225  parameter( pipipi = 3.141592653589793238462643383279d+00 )
1226  parameter( twopip = 6.283185307179586476925286766559d+00 )
1227  parameter( pip5o2 = 7.853981633974483096156608458199d+00 )
1228  parameter( pipisq = 9.869604401089358618834490999876d+00 )
1229  parameter( pihalf = 1.570796326794896619231321691640d+00 )
1230  parameter( erfa00 = 0.886226925452758013649083741671d+00 )
1231  parameter( eneper = 2.718281828459045235360287471353d+00 )
1232  parameter( sqrent = 1.648721270700128146848650787814d+00 )
1233  parameter( sqrsix = 2.449489742783178098197284074706d+00 )
1234  parameter( sqrsev = 2.645751311064590590501615753639d+00 )
1235  parameter( sqrt12 = 3.464101615137754587054892683012d+00 )
1236  parameter( clight = 2.99792458 d+10 )
1237  parameter( avogad = 6.0221367 d+23 )
1238  parameter( boltzm = 1.380658 d-23 )
1239  parameter( amelgr = 9.1093897 d-28 )
1240  parameter( plckbr = 1.05457266 d-27 )
1241  parameter( elccgs = 4.8032068 d-10 )
1242  parameter( elcmks = 1.60217733 d-19 )
1243  parameter( amugrm = 1.6605402 d-24 )
1244  parameter( ammumu = 0.113428913 d+00 )
1245  parameter( amprmu = 1.007276470 d+00 )
1246  parameter( amnemu = 1.008664904 d+00 )
1247  parameter( alpfsc = 7.2973530791728595 d-03 )
1248  parameter( fscto2 = 5.3251361962113614 d-05 )
1249  parameter( fscto3 = 3.8859399018437826 d-07 )
1250  parameter( fscto4 = 2.8357075508200407 d-09 )
1251  parameter( plabrc = 0.197327053 d+00 )
1252  parameter( amelct = 0.51099906 d-03 )
1253  parameter( amugev = 0.93149432 d+00 )
1254  parameter( ammuon = 0.105658389 d+00 )
1255  parameter( amprtn = 0.93827231 d+00 )
1256  parameter( amntrn = 0.93956563 d+00 )
1257  parameter( amdeut = 1.87561339 d+00 )
1258  parameter( cougfm = elccgs * elccgs / elcmks * 1.d-07 * 1.d+13
1259  & * 1.d-09 )
1260  parameter( rclsel = 2.8179409183694872 d-13 )
1261  parameter( bltzmn = 8.617385 d-14 )
1262  parameter( gevmev = 1.0 d+03 )
1263  parameter( emvgev = 1.0 d-03 )
1264  parameter( algvmv = 6.90775527898214 d+00 )
1265  parameter( raddeg = 180.d+00 / pipipi )
1266  parameter( degrad = pipipi / 180.d+00 )
1267  LOGICAL lgbias, lgbana
1268  COMMON / global / lgbias, lgbana
1269 C INCLUDE '(DIMPAR)'
1270 *$ CREATE DIMPAR.ADD
1271  parameter( mxxrgn = 5000 )
1272  parameter( mxxmdf = 56 )
1273  parameter( mxxmde = 50 )
1274  parameter( mfstck = 1000 )
1275  parameter( mestck = 100 )
1276  parameter( nallwp = 39 )
1277  parameter( mpdpdx = 8 )
1278  parameter( icomax = 180 )
1279  parameter( nstbis = 304 )
1280  parameter( idmaxp = 210 )
1281  parameter( idmxdc = 620 )
1282  parameter( mkbmx1 = 1 )
1283  parameter( mkbmx2 = 1 )
1284 C INCLUDE '(IOUNIT)'
1285 *$ CREATE IOUNIT.ADD
1286  parameter( lunin = 5 )
1287  parameter( lunout = 6 )
1288  parameter( lunerr = 15 )
1289  parameter( lunber = 14 )
1290  parameter( lunech = 8 )
1291  parameter( lunflu = 13 )
1292  parameter( lungeo = 16 )
1293  parameter( lunpgs = 12 )
1294  parameter( lunran = 2 )
1295  parameter( lunxsc = 9 )
1296  parameter( lundet = 17 )
1297  parameter( lunray = 10 )
1298  parameter( lunrdb = 1 )
1299 *
1300 * IMPLICIT DOUBLE PRECISION (A-H,O-Z)
1301 *
1302  dimension rvec(*)
1303  common/r48st1/u(97),c,i97,j97
1304  parameter(modcns=1000000000)
1305  SAVE cd, cm, twom24, zero, one, ntot, ntot2, ijkl
1306  DATA ntot,ntot2,ijkl/-1,0,0/
1307 C
1308  IF (ntot .GE. 0) go to 50
1309 C
1310 C Default initialization. User has called RM48 without RM48IN.
1311  ijkl = 54217137
1312  ntot = 0
1313  ntot2 = 0
1314  kalled = 0
1315  go to 1
1316 C
1317  entry rm48in(ijklin, ntotin,ntot2n)
1318 C Initializing routine for RM48, may be called before
1319 C generating pseudorandom numbers with RM48. The input
1320 C values should be in the ranges: 0<=IJKLIN<=900 OOO OOO
1321 C 0<=NTOTIN<=999 999 999
1322 C 0<=NTOT2N<<999 999 999!
1323 C To get the standard values in Marsaglia's paper, IJKLIN=54217137
1324 C NTOTIN,NTOT2N=0
1325  ijkl = ijklin
1326  ntot = max(ntotin,0)
1327  ntot2= max(ntot2n,0)
1328  kalled = 1
1329 C always come here to initialize
1330  1 CONTINUE
1331  ij = ijkl/30082
1332  kl = ijkl - 30082*ij
1333  i = mod(ij/177, 177) + 2
1334  j = mod(ij, 177) + 2
1335  k = mod(kl/169, 178) + 1
1336  l = mod(kl, 169)
1337  WRITE(lunout,'(A,I10,2X,2I10)')
1338  & ' RM48 INITIALIZED:',ijkl,ntot,ntot2
1339 CCC PRINT '(A,4I10)', ' I,J,K,L= ',I,J,K,L
1340  one = 1.d+00
1341  half = 0.5d+00
1342  zero = 0.d+00
1343  DO 2 ii= 1, 97
1344  s = 0.d+00
1345  t = half
1346  DO 3 jj= 1, 48
1347  m = mod(mod(i*j,179)*k, 179)
1348  i = j
1349  j = k
1350  k = m
1351  l = mod(53*l+1, 169)
1352  IF (mod(l*m,64) .GE. 32) s = s+t
1353  3 t = half*t
1354  2 u(ii) = s
1355  twom24 = one
1356  DO 4 i24= 1, 24
1357  4 twom24 = half*twom24
1358  c = 362436.d+00*twom24
1359  cd = 7654321.d+00*twom24
1360  cm = 16777213.d+00*twom24
1361  i97 = 97
1362  j97 = 33
1363 C Complete initialization by skipping
1364 C (NTOT2*MODCNS + NTOT) random numbers
1365  DO 45 loop2= 1, ntot2+1
1366  now = modcns
1367  IF (loop2 .EQ. ntot2+1) now=ntot
1368  IF (now .GT. 0) THEN
1369  WRITE(lunout,'(A,I15)') ' RM48IN SKIPPING OVER ',now
1370  DO 40 idum = 1, ntot
1371  uni = u(i97)-u(j97)
1372  IF (uni .LT. zero) uni=uni+one
1373  u(i97) = uni
1374  i97 = i97-1
1375  IF (i97 .EQ. 0) i97=97
1376  j97 = j97-1
1377  IF (j97 .EQ. 0) j97=97
1378  c = c - cd
1379  IF (c .LT. zero) c=c+cm
1380  40 CONTINUE
1381  ENDIF
1382  45 CONTINUE
1383  IF (kalled .EQ. 1) RETURN
1384 C
1385 C Normal entry to generate LENV random numbers
1386  50 CONTINUE
1387  DO 100 ivec= 1, lenv
1388  uni = u(i97)-u(j97)
1389  IF (uni .LT. zero) uni=uni+one
1390  u(i97) = uni
1391  i97 = i97-1
1392  IF (i97 .EQ. 0) i97=97
1393  j97 = j97-1
1394  IF (j97 .EQ. 0) j97=97
1395  c = c - cd
1396  IF (c .LT. zero) c=c+cm
1397  uni = uni-c
1398  IF (uni .LT. zero) uni=uni+one
1399  rvec(ivec) = uni
1400 CC Replace exact zeros by uniform distr. *2**-24
1401 C IF (UNI .EQ. 0.) THEN
1402 C ZUNI = TWOM24*U(2)
1403 CC An exact zero here is very unlikely, but let's be safe.
1404 C IF (ZUNI .EQ. 0.) ZUNI= TWOM24*TWOM24
1405 C RVEC(IVEC) = ZUNI
1406 C ENDIF
1407  100 CONTINUE
1408  ntot = ntot + lenv
1409  IF (ntot .GE. modcns) THEN
1410  ntot2 = ntot2 + 1
1411  ntot = ntot - modcns
1412  ENDIF
1413  RETURN
1414 C Entry to output current status
1415  entry rm48ut(ijklut,ntotut,ntot2t)
1416  ijklut = ijkl
1417  ntotut = ntot
1418  ntot2t = ntot2
1419  RETURN
1420 C
1421  entry rm48wr(ioseed)
1422 C Output routine for RM48, without skipping numbers
1423  WRITE (ioseed,'(2Z8)') ntot,ntot2
1424  WRITE (ioseed,'(2Z8,Z16)') i97,j97,c
1425  WRITE (ioseed,'(24(4Z16,/),Z16)') u
1426  RETURN
1427 C
1428  entry rm48rd(ioseed)
1429 C Initializing routine for RM48, without skipping numbers
1430  READ (ioseed,'(2Z8)') ntot,ntot2
1431  READ (ioseed,'(2Z8,Z16)') i97,j97,c
1432  READ (ioseed,'(24(4Z16,/),Z16)') u
1433  CLOSE (unit=ioseed)
1434  ijkl = 54217137
1435  ij = ijkl/30082
1436  kl = ijkl - 30082*ij
1437  i = mod(ij/177, 177) + 2
1438  j = mod(ij, 177) + 2
1439  k = mod(kl/169, 178) + 1
1440  l = mod(kl, 169)
1441  WRITE (lunout,'(A,I10,2X,2I10)')
1442  & ' RM48 INITIALIZED:',ijkl,ntot,ntot2
1443 CCC PRINT '(A,4I10)', ' I,J,K,L= ',I,J,K,L
1444  one = 1.d+00
1445  half = 0.5d+00
1446  zero = 0.d+00
1447  twom24 = one
1448  DO 400 i24= 1, 24
1449  400 twom24 = half*twom24
1450  cd = 7654321.d+00*twom24
1451  cm = 16777213.d+00*twom24
1452  RETURN
1453  END
1454 
1455 *$ CREATE RNDM.FOR
1456 *COPY RNDM
1457 *
1458 *=== rndm ===========================================================*
1459 *
1460  DOUBLE PRECISION FUNCTION rndm (RDUMMY)
1461 
1462 C INCLUDE '(DBLPRC)'
1463 *$ CREATE DBLPRC.ADD
1464  IMPLICIT DOUBLE PRECISION (a-h,o-z)
1465  parameter( kalgnm = 2 )
1466  parameter( anglgb = 5.0d-16 )
1467  parameter( anglsq = 2.5d-31 )
1468  parameter( axcssv = 0.2d+16 )
1469  parameter( andrfl = 1.0d-38 )
1470  parameter( avrflw = 1.0d+38 )
1471  parameter( ainfnt = 1.0d+30 )
1472  parameter( azrzrz = 1.0d-30 )
1473  parameter( einfnt = +69.07755278982137 d+00 )
1474  parameter( ezrzrz = -69.07755278982137 d+00 )
1475  parameter( onemns = 0.999999999999999 d+00 )
1476  parameter( onepls = 1.000000000000001 d+00 )
1477  parameter( csnnrm = 2.0d-15 )
1478  parameter( dmxtrn = 1.0d+08 )
1479  parameter( zerzer = 0.d+00 )
1480  parameter( oneone = 1.d+00 )
1481  parameter( twotwo = 2.d+00 )
1482  parameter( thrthr = 3.d+00 )
1483  parameter( foufou = 4.d+00 )
1484  parameter( fivfiv = 5.d+00 )
1485  parameter( sixsix = 6.d+00 )
1486  parameter( sevsev = 7.d+00 )
1487  parameter( eigeig = 8.d+00 )
1488  parameter( aninen = 9.d+00 )
1489  parameter( tenten = 10.d+00 )
1490  parameter( hlfhlf = 0.5d+00 )
1491  parameter( onethi = oneone / thrthr )
1492  parameter( twothi = twotwo / thrthr )
1493  parameter( onefou = oneone / foufou )
1494  parameter( thrtwo = thrthr / twotwo )
1495  parameter( pipipi = 3.141592653589793238462643383279d+00 )
1496  parameter( twopip = 6.283185307179586476925286766559d+00 )
1497  parameter( pip5o2 = 7.853981633974483096156608458199d+00 )
1498  parameter( pipisq = 9.869604401089358618834490999876d+00 )
1499  parameter( pihalf = 1.570796326794896619231321691640d+00 )
1500  parameter( erfa00 = 0.886226925452758013649083741671d+00 )
1501  parameter( eneper = 2.718281828459045235360287471353d+00 )
1502  parameter( sqrent = 1.648721270700128146848650787814d+00 )
1503  parameter( sqrsix = 2.449489742783178098197284074706d+00 )
1504  parameter( sqrsev = 2.645751311064590590501615753639d+00 )
1505  parameter( sqrt12 = 3.464101615137754587054892683012d+00 )
1506  parameter( clight = 2.99792458 d+10 )
1507  parameter( avogad = 6.0221367 d+23 )
1508  parameter( boltzm = 1.380658 d-23 )
1509  parameter( amelgr = 9.1093897 d-28 )
1510  parameter( plckbr = 1.05457266 d-27 )
1511  parameter( elccgs = 4.8032068 d-10 )
1512  parameter( elcmks = 1.60217733 d-19 )
1513  parameter( amugrm = 1.6605402 d-24 )
1514  parameter( ammumu = 0.113428913 d+00 )
1515  parameter( amprmu = 1.007276470 d+00 )
1516  parameter( amnemu = 1.008664904 d+00 )
1517  parameter( alpfsc = 7.2973530791728595 d-03 )
1518  parameter( fscto2 = 5.3251361962113614 d-05 )
1519  parameter( fscto3 = 3.8859399018437826 d-07 )
1520  parameter( fscto4 = 2.8357075508200407 d-09 )
1521  parameter( plabrc = 0.197327053 d+00 )
1522  parameter( amelct = 0.51099906 d-03 )
1523  parameter( amugev = 0.93149432 d+00 )
1524  parameter( ammuon = 0.105658389 d+00 )
1525  parameter( amprtn = 0.93827231 d+00 )
1526  parameter( amntrn = 0.93956563 d+00 )
1527  parameter( amdeut = 1.87561339 d+00 )
1528  parameter( cougfm = elccgs * elccgs / elcmks * 1.d-07 * 1.d+13
1529  & * 1.d-09 )
1530  parameter( rclsel = 2.8179409183694872 d-13 )
1531  parameter( bltzmn = 8.617385 d-14 )
1532  parameter( gevmev = 1.0 d+03 )
1533  parameter( emvgev = 1.0 d-03 )
1534  parameter( algvmv = 6.90775527898214 d+00 )
1535  parameter( raddeg = 180.d+00 / pipipi )
1536  parameter( degrad = pipipi / 180.d+00 )
1537  LOGICAL lgbias, lgbana
1538  COMMON / global / lgbias, lgbana
1539 C INCLUDE '(DIMPAR)'
1540 *$ CREATE DIMPAR.ADD
1541  parameter( mxxrgn = 5000 )
1542  parameter( mxxmdf = 56 )
1543  parameter( mxxmde = 50 )
1544  parameter( mfstck = 1000 )
1545  parameter( mestck = 100 )
1546  parameter( nallwp = 39 )
1547  parameter( mpdpdx = 8 )
1548  parameter( icomax = 180 )
1549  parameter( nstbis = 304 )
1550  parameter( idmaxp = 210 )
1551  parameter( idmxdc = 620 )
1552  parameter( mkbmx1 = 1 )
1553  parameter( mkbmx2 = 1 )
1554 C INCLUDE '(IOUNIT)'
1555 *$ CREATE IOUNIT.ADD
1556  parameter( lunin = 5 )
1557  parameter( lunout = 6 )
1558  parameter( lunerr = 15 )
1559  parameter( lunber = 14 )
1560  parameter( lunech = 8 )
1561  parameter( lunflu = 13 )
1562  parameter( lungeo = 16 )
1563  parameter( lunpgs = 12 )
1564  parameter( lunran = 2 )
1565  parameter( lunxsc = 9 )
1566  parameter( lundet = 17 )
1567  parameter( lunray = 10 )
1568  parameter( lunrdb = 1 )
1569 *
1570 *----------------------------------------------------------------------*
1571 * *
1572 * This routine merely acts as an interface to F.James RM48 gen. *
1573 * *
1574 * Created on 03 april 1992 by Alfredo Ferrari & Paola Sala *
1575 * Infn - Milan *
1576 * *
1577 * Last change on 16-sep-93 by Alfredo Ferrari *
1578 * *
1579 * *
1580 *----------------------------------------------------------------------*
1581 *
1582  dimension rndnum(2)
1583  CALL rm48( rndnum, 1 )
1584  rndm = rndnum(1)
1585  RETURN
1586  entry rd2in(iseed1,iseed2)
1587 * The following card just to avoid warning messages on the HP compiler
1588  rd2in = pipipi
1589  CALL rm48in(54217137,iseed1,iseed2)
1590  RETURN
1591  entry rd2out(iseed1,iseed2)
1592 * The following card just to avoid warning messages on the HP compiler
1593  rd2out = pipipi
1594  CALL rm48ut(idummy,iseed1,iseed2)
1595 *=== End of function rndm =============================================*
1596  RETURN
1597  END
1598 
Generated on Sat May 25 2013 14:32:37 for Geant4 by   doxygen 1.8.4