← Back to branch summary

~maddevelopers/mg5amcnlo/2.5.4_run_py8_at_evtgen

« back to all changes in this revision

Viewing changes to tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%V0_udx_wp%check_sa_born_splitOrders.f

  • Committer: olivier Mattelaer
  • Date: 2016-05-12 11:00:18 UTC
  • mfrom: (262.1.150 2.3.4)
  • Revision ID: olivier.mattelaer@uclouvain.be-20160512110018-sevb79f0wm4g8mpp
pass to 2.4.0

Show diffs side-by-side

added added

removed removed

Lines of Context:
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%V0_udx_wp%check_sa_born_splitOrders.f
 
1
      PROGRAM DRIVER
 
2
C     *****************************************************************
 
3
C     *********
 
4
C     THIS IS THE DRIVER FOR CHECKING THE STANDALONE MATRIX ELEMENT.
 
5
C     IT USES A SIMPLE PHASE SPACE GENERATOR
 
6
C     Fabio Maltoni - 3rd Febraury 2007
 
7
C     *****************************************************************
 
8
C     *********
 
9
      IMPLICIT NONE
 
10
C     
 
11
C     CONSTANTS  
 
12
C     
 
13
      REAL*8 ZERO
 
14
      PARAMETER (ZERO=0D0)
 
15
      INTEGER NSPLITORDERS
 
16
      PARAMETER (NSPLITORDERS=1)
 
17
C     
 
18
C     INCLUDE FILES
 
19
C     
 
20
C     ---  the include file with the values of the parameters and
 
21
C      masses   
 
22
      INCLUDE 'coupl.inc'
 
23
C     integer nexternal and number particles (incoming+outgoing) in
 
24
C      the me 
 
25
      INTEGER NEXTERNAL, NINCOMING
 
26
      PARAMETER (NEXTERNAL=3,NINCOMING=2)
 
27
C     ---  particle masses
 
28
      REAL*8 PMASS(NEXTERNAL)
 
29
      REAL*8 TOTALMASS
 
30
C     ---  integer    n_max_cg
 
31
      INCLUDE 'ngraphs.inc'  !how many diagrams (could be useful to know...)
 
32
 
 
33
C     
 
34
C     LOCAL
 
35
C     
 
36
      INTEGER I,J,K
 
37
      REAL*8 P(0:3,NEXTERNAL)  ! four momenta. Energy is the zeroth component.
 
38
      REAL*8 SQRTS  ! sqrt(s)= center of mass energy
 
39
      REAL*8 MATELEM, MATELEMS(0:NSPLITORDERS)
 
40
      REAL*8 PIN(0:3), POUT(0:3)
 
41
      CHARACTER*120 BUFF(NEXTERNAL)
 
42
 
 
43
      INTEGER NCHOSEN
 
44
      CHARACTER*20 CHOSEN_SO_INDICES(NSPLITORDERS)
 
45
 
 
46
C     
 
47
C     EXTERNAL
 
48
C     
 
49
      REAL*8 DOT
 
50
      EXTERNAL DOT
 
51
 
 
52
      LOGICAL CHOSEN_SO_CONFIGS(NSPLITORDERS)
 
53
      COMMON/CHOSEN_BORN_SQSO/CHOSEN_SO_CONFIGS
 
54
 
 
55
C     -----
 
56
C     BEGIN CODE
 
57
C     -----
 
58
C     
 
59
C     Start by initializing what is the squared split orders indices
 
60
C      chosen
 
61
      NCHOSEN=0
 
62
      DO I=1,NSPLITORDERS
 
63
        IF (CHOSEN_SO_CONFIGS(I)) THEN
 
64
          NCHOSEN=NCHOSEN+1
 
65
          WRITE(CHOSEN_SO_INDICES(NCHOSEN),'(I3,A1)') I,')'
 
66
        ENDIF
 
67
      ENDDO
 
68
 
 
69
C     ---  INITIALIZATION CALLS
 
70
C     
 
71
C     ---  Call to initialize the values of the couplings, masses and
 
72
C      widths 
 
73
C     used in the evaluation of the matrix element. The primary
 
74
C      parameters of the
 
75
C     models are read from Cards/param_card.dat. The secondary
 
76
C      parameters are calculated
 
77
C     in Source/MODEL/couplings.f. The values are stored in common
 
78
C      blocks that are listed
 
79
C     in coupl.inc .
 
80
 
 
81
      CALL SETPARA('param_card.dat')  !first call to setup the paramaters
 
82
      INCLUDE 'pmass.inc'  !set up masses
 
83
 
 
84
      TOTALMASS = 0.0D0
 
85
      DO I=1,NEXTERNAL
 
86
        TOTALMASS = TOTALMASS + PMASS(I)
 
87
      ENDDO
 
88
 
 
89
C     ---  Now use a simple multipurpose PS generator (RAMBO) just to
 
90
C      get a 
 
91
C     RANDOM set of four momenta of given masses pmass(i) to be used
 
92
C      to evaluate 
 
93
C     the madgraph matrix-element.       
 
94
C     Alternatevely, here the user can call or set the four momenta at
 
95
C      his will, see below.
 
96
C     
 
97
      IF(NINCOMING.EQ.1) THEN
 
98
        SQRTS=PMASS(1)
 
99
      ELSE
 
100
        SQRTS=1000D0  !CMS energy in GEV
 
101
        IF (SQRTS.LE.2.0D0*TOTALMASS) THEN
 
102
          SQRTS = 2.0D0*TOTALMASS
 
103
        ENDIF
 
104
      ENDIF
 
105
 
 
106
      CALL PRINTOUT()
 
107
 
 
108
      CALL GET_MOMENTA(SQRTS,PMASS,P)
 
109
C     
 
110
C     write the information on the four momenta 
 
111
C     
 
112
      WRITE (*,*)
 
113
      WRITE (*,*) ' Phase space point:'
 
114
      WRITE (*,*)
 
115
      WRITE (*,*) '-----------------------------'
 
116
      WRITE (*,*)  'n   E   px   py   pz   m '
 
117
      DO I=1,NEXTERNAL
 
118
        WRITE (*,'(i2,1x,5e15.7)') I, P(0,I),P(1,I),P(2,I),P(3,I)
 
119
     $   ,DSQRT(DABS(DOT(P(0,I),P(0,I))))
 
120
      ENDDO
 
121
      WRITE (*,*) '-----------------------------'
 
122
 
 
123
C     
 
124
C     Now we can call the matrix element!
 
125
C     
 
126
      CALL SMATRIX_SPLITORDERS(P,MATELEMS)
 
127
      MATELEM=MATELEMS(0)
 
128
      WRITE(*,*) '1) Matrix element for (QCD=0) = ',MATELEMS(1)
 
129
C     
 
130
      IF (NCHOSEN.NE.NSPLITORDERS) THEN
 
131
        WRITE (*,*) 'Selected squared coupling orders combination for'
 
132
     $   //' the sum below:'
 
133
        WRITE (*,*) (CHOSEN_SO_INDICES(I),I=1,NCHOSEN)
 
134
      ENDIF
 
135
      WRITE (*,*) 'Total Matrix element = ', MATELEM, ' GeV^',-(2
 
136
     $ *NEXTERNAL-8)
 
137
      WRITE (*,*) '-----------------------------'
 
138
 
 
139
C     c
 
140
C     c      Copy down here (or read in) the four momenta as a string. 
 
141
C     c      
 
142
C     c
 
143
C     buff(1)=" 1   0.5630480E+04  0.0000000E+00  0.0000000E+00 
 
144
C      0.5630480E+04"
 
145
C     buff(2)=" 2   0.5630480E+04  0.0000000E+00  0.0000000E+00
 
146
C      -0.5630480E+04"
 
147
C     buff(3)=" 3   0.5466073E+04  0.4443190E+03  0.2446331E+04
 
148
C      -0.4864732E+04"
 
149
C     buff(4)=" 4   0.8785819E+03 -0.2533886E+03  0.2741971E+03 
 
150
C      0.7759741E+03"
 
151
C     buff(5)=" 5   0.4916306E+04 -0.1909305E+03 -0.2720528E+04 
 
152
C      0.4088757E+04"
 
153
C     c
 
154
C     c      Here the k,E,px,py,pz are read from the string into the
 
155
C      momenta array.
 
156
C     c      k=1,2          : incoming
 
157
C     c      k=3,nexternal  : outgoing
 
158
C     c
 
159
C     do i=1,nexternal
 
160
C     read (buff(i),*) k, P(0,i),P(1,i),P(2,i),P(3,i)
 
161
C     enddo
 
162
C     
 
163
C     c- print the momenta out
 
164
C     
 
165
C     do i=1,nexternal
 
166
C     write (*,'(i2,1x,5e15.7)') i, P(0,i),P(1,i),P(2,i),P(3,i), 
 
167
C     .dsqrt(dabs(DOT(p(0,i),p(0,i))))
 
168
C     enddo
 
169
C     
 
170
C     CALL SMATRIX(P,MATELEM)
 
171
C     
 
172
C     write (*,*) "-------------------------------------------------"
 
173
C     write (*,*) "Matrix element = ", MATELEM, " GeV^",-(2*nexternal-8)
 
174
C       
 
175
C     write (*,*) "-------------------------------------------------"
 
176
 
 
177
      END
 
178
 
 
179
 
 
180
 
 
181
 
 
182
      DOUBLE PRECISION FUNCTION DOT(P1,P2)
 
183
C     *****************************************************************
 
184
C     ***********
 
185
C     4-Vector Dot product
 
186
C     *****************************************************************
 
187
C     ***********
 
188
      IMPLICIT NONE
 
189
      DOUBLE PRECISION P1(0:3),P2(0:3)
 
190
      DOT=P1(0)*P2(0)-P1(1)*P2(1)-P1(2)*P2(2)-P1(3)*P2(3)
 
191
      END
 
192
 
 
193
 
 
194
      SUBROUTINE GET_MOMENTA(ENERGY,PMASS,P)
 
195
C     ---- auxiliary function to change convention between madgraph
 
196
C      and rambo
 
197
C     ---- four momenta.          
 
198
      IMPLICIT NONE
 
199
C     integer nexternal and number particles (incoming+outgoing) in
 
200
C      the me 
 
201
      INTEGER NEXTERNAL, NINCOMING
 
202
      PARAMETER (NEXTERNAL=3,NINCOMING=2)
 
203
C     ARGUMENTS
 
204
      REAL*8 ENERGY,PMASS(NEXTERNAL),P(0:3,NEXTERNAL),PRAMBO(4,10),WGT
 
205
C     LOCAL
 
206
      INTEGER I
 
207
      REAL*8 ETOT2,MOM,M1,M2,E1,E2
 
208
 
 
209
      ETOT2=ENERGY**2
 
210
      M1=PMASS(1)
 
211
      M2=PMASS(2)
 
212
      MOM=(ETOT2**2 - 2*ETOT2*M1**2 + M1**4 - 2*ETOT2*M2**2 - 2*M1**2
 
213
     $ *M2**2 + M2**4)/(4.*ETOT2)
 
214
      MOM=DSQRT(MOM)
 
215
      E1=DSQRT(MOM**2+M1**2)
 
216
      E2=DSQRT(MOM**2+M2**2)
 
217
      WRITE (*,*) E1+E2,MOM
 
218
 
 
219
      IF(NINCOMING.EQ.2) THEN
 
220
 
 
221
        P(0,1)=E1
 
222
        P(1,1)=0D0
 
223
        P(2,1)=0D0
 
224
        P(3,1)=MOM
 
225
 
 
226
        P(0,2)=E2
 
227
        P(1,2)=0D0
 
228
        P(2,2)=0D0
 
229
        P(3,2)=-MOM
 
230
 
 
231
        CALL RAMBO(NEXTERNAL-2,ENERGY,PMASS(3),PRAMBO,WGT)
 
232
        DO I=3, NEXTERNAL
 
233
          P(0,I)=PRAMBO(4,I-2)
 
234
          P(1,I)=PRAMBO(1,I-2)
 
235
          P(2,I)=PRAMBO(2,I-2)
 
236
          P(3,I)=PRAMBO(3,I-2)
 
237
        ENDDO
 
238
 
 
239
      ELSEIF(NINCOMING.EQ.1) THEN
 
240
 
 
241
        P(0,1)=ENERGY
 
242
        P(1,1)=0D0
 
243
        P(2,1)=0D0
 
244
        P(3,1)=0D0
 
245
 
 
246
        CALL RAMBO(NEXTERNAL-1,ENERGY,PMASS(2),PRAMBO,WGT)
 
247
        DO I=2, NEXTERNAL
 
248
          P(0,I)=PRAMBO(4,I-1)
 
249
          P(1,I)=PRAMBO(1,I-1)
 
250
          P(2,I)=PRAMBO(2,I-1)
 
251
          P(3,I)=PRAMBO(3,I-1)
 
252
        ENDDO
 
253
      ENDIF
 
254
 
 
255
      RETURN
 
256
      END
 
257
 
 
258
 
 
259
      SUBROUTINE RAMBO(N,ET,XM,P,WT)
 
260
C     *****************************************************************
 
261
C     ******
 
262
C     *                       RAMBO                                   
 
263
C           *
 
264
C     *    RA(NDOM)  M(OMENTA)  B(EAUTIFULLY)  O(RGANIZED)            
 
265
C           *
 
266
C     *                                                               
 
267
C           *
 
268
C     *    A DEMOCRATIC MULTI-PARTICLE PHASE SPACE GENERATOR          
 
269
C           *
 
270
C     *    AUTHORS:  S.D. ELLIS,  R. KLEISS,  W.J. STIRLING           
 
271
C           *
 
272
C     *    THIS IS VERSION 1.0 -  WRITTEN BY R. KLEISS                
 
273
C           *
 
274
C     *    -- ADJUSTED BY HANS KUIJF, WEIGHTS ARE LOGARITHMIC
 
275
C      (20-08-90)    *
 
276
C     *                                                               
 
277
C           *
 
278
C     *    N  = NUMBER OF PARTICLES                                   
 
279
C           *
 
280
C     *    ET = TOTAL CENTRE-OF-MASS ENERGY                           
 
281
C           *
 
282
C     *    XM = PARTICLE MASSES ( DIM=NEXTERNAL-nincoming )           
 
283
C           *
 
284
C     *    P  = PARTICLE MOMENTA ( DIM=(4,NEXTERNAL-nincoming) )      
 
285
C           *
 
286
C     *    WT = WEIGHT OF THE EVENT                                   
 
287
C           *
 
288
C     *****************************************************************
 
289
C     ******
 
290
      IMPLICIT REAL*8(A-H,O-Z)
 
291
C     integer nexternal and number particles (incoming+outgoing) in
 
292
C      the me 
 
293
      INTEGER NEXTERNAL, NINCOMING
 
294
      PARAMETER (NEXTERNAL=3,NINCOMING=2)
 
295
      DIMENSION XM(NEXTERNAL-NINCOMING),P(4,NEXTERNAL-NINCOMING)
 
296
      DIMENSION Q(4,NEXTERNAL-NINCOMING),Z(NEXTERNAL-NINCOMING),R(4),
 
297
     $  B(3),P2(NEXTERNAL-NINCOMING),XM2(NEXTERNAL-NINCOMING),
 
298
     $  E(NEXTERNAL-NINCOMING),V(NEXTERNAL-NINCOMING),IWARN(5)
 
299
      SAVE ACC,ITMAX,IBEGIN,IWARN
 
300
      DATA ACC/1.D-14/,ITMAX/6/,IBEGIN/0/,IWARN/5*0/
 
301
C     *
 
302
C     * INITIALIZATION STEP: FACTORIALS FOR THE PHASE SPACE WEIGHT
 
303
      IF(IBEGIN.NE.0) GOTO 103
 
304
      IBEGIN=1
 
305
      TWOPI=8.*DATAN(1.D0)
 
306
      PO2LOG=LOG(TWOPI/4.)
 
307
      Z(2)=PO2LOG
 
308
      DO 101 K=3,NEXTERNAL-NINCOMING-1
 
309
 101  Z(K)=Z(K-1)+PO2LOG-2.*LOG(DFLOAT(K-2))
 
310
      DO 102 K=3,NEXTERNAL-NINCOMING-1
 
311
 102  Z(K)=(Z(K)-LOG(DFLOAT(K-1)))
 
312
C     *
 
313
C     * CHECK ON THE NUMBER OF PARTICLES
 
314
 103  IF(N.GT.1.AND.N.LT.101) GOTO 104
 
315
      PRINT 1001,N
 
316
      STOP
 
317
C     *
 
318
C     * CHECK WHETHER TOTAL ENERGY IS SUFFICIENT; COUNT NONZERO MASSES
 
319
 104  XMT=0.
 
320
      NM=0
 
321
      DO 105 I=1,N
 
322
      IF(XM(I).NE.0.D0) NM=NM+1
 
323
 105  XMT=XMT+ABS(XM(I))
 
324
      IF(XMT.LE.ET) GOTO 201
 
325
      PRINT 1002,XMT,ET
 
326
      STOP
 
327
C     *
 
328
C     * THE PARAMETER VALUES ARE NOW ACCEPTED
 
329
C     *
 
330
C     * GENERATE N MASSLESS MOMENTA IN INFINITE PHASE SPACE
 
331
 201  DO 202 I=1,N
 
332
      R1=RN(1)
 
333
      C=2.*R1-1.
 
334
      S=SQRT(1.-C*C)
 
335
      F=TWOPI*RN(2)
 
336
      R1=RN(3)
 
337
      R2=RN(4)
 
338
      Q(4,I)=-LOG(R1*R2)
 
339
      Q(3,I)=Q(4,I)*C
 
340
      Q(2,I)=Q(4,I)*S*COS(F)
 
341
 202  Q(1,I)=Q(4,I)*S*SIN(F)
 
342
C     *
 
343
C     * CALCULATE THE PARAMETERS OF THE CONFORMAL TRANSFORMATION
 
344
      DO 203 I=1,4
 
345
 203  R(I)=0.
 
346
      DO 204 I=1,N
 
347
      DO 204 K=1,4
 
348
 204  R(K)=R(K)+Q(K,I)
 
349
      RMAS=SQRT(R(4)**2-R(3)**2-R(2)**2-R(1)**2)
 
350
      DO 205 K=1,3
 
351
 205  B(K)=-R(K)/RMAS
 
352
      G=R(4)/RMAS
 
353
      A=1./(1.+G)
 
354
      X=ET/RMAS
 
355
C     *
 
356
C     * TRANSFORM THE Q'S CONFORMALLY INTO THE P'S
 
357
      DO 207 I=1,N
 
358
      BQ=B(1)*Q(1,I)+B(2)*Q(2,I)+B(3)*Q(3,I)
 
359
      DO 206 K=1,3
 
360
 206  P(K,I)=X*(Q(K,I)+B(K)*(Q(4,I)+A*BQ))
 
361
 207  P(4,I)=X*(G*Q(4,I)+BQ)
 
362
C     *
 
363
C     * CALCULATE WEIGHT AND POSSIBLE WARNINGS
 
364
      WT=PO2LOG
 
365
      IF(N.NE.2) WT=(2.*N-4.)*LOG(ET)+Z(N)
 
366
      IF(WT.GE.-180.D0) GOTO 208
 
367
      IF(IWARN(1).LE.5) PRINT 1004,WT
 
368
      IWARN(1)=IWARN(1)+1
 
369
 208  IF(WT.LE. 174.D0) GOTO 209
 
370
      IF(IWARN(2).LE.5) PRINT 1005,WT
 
371
      IWARN(2)=IWARN(2)+1
 
372
C     *
 
373
C     * RETURN FOR WEIGHTED MASSLESS MOMENTA
 
374
 209  IF(NM.NE.0) GOTO 210
 
375
C     * RETURN LOG OF WEIGHT
 
376
      WT=WT
 
377
      RETURN
 
378
C     *
 
379
C     * MASSIVE PARTICLES: RESCALE THE MOMENTA BY A FACTOR X
 
380
 210  XMAX=SQRT(1.-(XMT/ET)**2)
 
381
      DO 301 I=1,N
 
382
      XM2(I)=XM(I)**2
 
383
 301  P2(I)=P(4,I)**2
 
384
      ITER=0
 
385
      X=XMAX
 
386
      ACCU=ET*ACC
 
387
 302  F0=-ET
 
388
      G0=0.
 
389
      X2=X*X
 
390
      DO 303 I=1,N
 
391
      E(I)=SQRT(XM2(I)+X2*P2(I))
 
392
      F0=F0+E(I)
 
393
 303  G0=G0+P2(I)/E(I)
 
394
      IF(ABS(F0).LE.ACCU) GOTO 305
 
395
      ITER=ITER+1
 
396
      IF(ITER.LE.ITMAX) GOTO 304
 
397
      PRINT 1006,ITMAX
 
398
      GOTO 305
 
399
 304  X=X-F0/(X*G0)
 
400
      GOTO 302
 
401
 305  DO 307 I=1,N
 
402
      V(I)=X*P(4,I)
 
403
      DO 306 K=1,3
 
404
 306  P(K,I)=X*P(K,I)
 
405
 307  P(4,I)=E(I)
 
406
C     *
 
407
C     * CALCULATE THE MASS-EFFECT WEIGHT FACTOR
 
408
      WT2=1.
 
409
      WT3=0.
 
410
      DO 308 I=1,N
 
411
      WT2=WT2*V(I)/E(I)
 
412
 308  WT3=WT3+V(I)**2/E(I)
 
413
      WTM=(2.*N-3.)*LOG(X)+LOG(WT2/WT3*ET)
 
414
C     *
 
415
C     * RETURN FOR  WEIGHTED MASSIVE MOMENTA
 
416
      WT=WT+WTM
 
417
      IF(WT.GE.-180.D0) GOTO 309
 
418
      IF(IWARN(3).LE.5) PRINT 1004,WT
 
419
      IWARN(3)=IWARN(3)+1
 
420
 309  IF(WT.LE. 174.D0) GOTO 310
 
421
      IF(IWARN(4).LE.5) PRINT 1005,WT
 
422
      IWARN(4)=IWARN(4)+1
 
423
C     * RETURN LOG OF WEIGHT
 
424
 310  WT=WT
 
425
      RETURN
 
426
C     *
 
427
 1001 FORMAT(' RAMBO FAILS: # OF PARTICLES =',I5,' IS NOT ALLOWED')
 
428
 1002 FORMAT(' RAMBO FAILS: TOTAL MASS =',D15.6,' IS NOT SMALLER THAN'
 
429
     $ //' TOTAL ENERGY =',D15.6)
 
430
 1004 FORMAT(' RAMBO WARNS: WEIGHT = EXP(',F20.9,') MAY UNDERFLOW')
 
431
 1005 FORMAT(' RAMBO WARNS: WEIGHT = EXP(',F20.9,') MAY  OVERFLOW')
 
432
 1006 FORMAT(' RAMBO WARNS:',I3,' ITERATIONS DID NOT GIVE THE DESIRED'
 
433
     $ //' ACCURACY =',D15.6)
 
434
      END
 
435
 
 
436
      FUNCTION RN(IDUMMY)
 
437
      REAL*8 RN,RAN
 
438
      SAVE INIT
 
439
      DATA INIT /1/
 
440
      IF (INIT.EQ.1) THEN
 
441
        INIT=0
 
442
        CALL RMARIN(1802,9373)
 
443
        END IF
 
444
C       *
 
445
 10     CALL RANMAR(RAN)
 
446
        IF (RAN.LT.1D-16) GOTO 10
 
447
        RN=RAN
 
448
C       *
 
449
        END
 
450
 
 
451
 
 
452
 
 
453
        SUBROUTINE RANMAR(RVEC)
 
454
C       *     -----------------
 
455
C       * Universal random number generator proposed by Marsaglia and
 
456
C        Zaman
 
457
C       * in report FSU-SCRI-87-50
 
458
C       * In this version RVEC is a double precision variable.
 
459
        IMPLICIT REAL*8(A-H,O-Z)
 
460
        COMMON/ RASET1 / RANU(97),RANC,RANCD,RANCM
 
461
        COMMON/ RASET2 / IRANMR,JRANMR
 
462
        SAVE /RASET1/,/RASET2/
 
463
        UNI = RANU(IRANMR) - RANU(JRANMR)
 
464
        IF(UNI .LT. 0D0) UNI = UNI + 1D0
 
465
        RANU(IRANMR) = UNI
 
466
        IRANMR = IRANMR - 1
 
467
        JRANMR = JRANMR - 1
 
468
        IF(IRANMR .EQ. 0) IRANMR = 97
 
469
        IF(JRANMR .EQ. 0) JRANMR = 97
 
470
        RANC = RANC - RANCD
 
471
        IF(RANC .LT. 0D0) RANC = RANC + RANCM
 
472
        UNI = UNI - RANC
 
473
        IF(UNI .LT. 0D0) UNI = UNI + 1D0
 
474
        RVEC = UNI
 
475
        END
 
476
 
 
477
        SUBROUTINE RMARIN(IJ,KL)
 
478
C       *     -----------------
 
479
C       * Initializing routine for RANMAR, must be called before
 
480
C        generating
 
481
C       * any pseudorandom numbers with RANMAR. The input values
 
482
C        should be in
 
483
C       * the ranges 0<=ij<=31328 ; 0<=kl<=30081
 
484
        IMPLICIT REAL*8(A-H,O-Z)
 
485
        COMMON/ RASET1 / RANU(97),RANC,RANCD,RANCM
 
486
        COMMON/ RASET2 / IRANMR,JRANMR
 
487
        SAVE /RASET1/,/RASET2/
 
488
C       * This shows correspondence between the simplified input seeds
 
489
C        IJ, KL
 
490
C       * and the original Marsaglia-Zaman seeds I,J,K,L.
 
491
C       * To get the standard values in the Marsaglia-Zaman paper
 
492
C        (i=12,j=34
 
493
C       * k=56,l=78) put ij=1802, kl=9373
 
494
        I = MOD( IJ/177 , 177 ) + 2
 
495
        J = MOD( IJ     , 177 ) + 2
 
496
        K = MOD( KL/169 , 178 ) + 1
 
497
        L = MOD( KL     , 169 )
 
498
        DO 300 II = 1 , 97
 
499
        S =  0D0
 
500
        T = .5D0
 
501
        DO 200 JJ = 1 , 24
 
502
        M = MOD( MOD(I*J,179)*K , 179 )
 
503
        I = J
 
504
        J = K
 
505
        K = M
 
506
        L = MOD( 53*L+1 , 169 )
 
507
        IF(MOD(L*M,64) .GE. 32) S = S + T
 
508
        T = .5D0*T
 
509
 200    CONTINUE
 
510
        RANU(II) = S
 
511
 300    CONTINUE
 
512
        RANC  =   362436D0 / 16777216D0
 
513
        RANCD =  7654321D0 / 16777216D0
 
514
        RANCM = 16777213D0 / 16777216D0
 
515
        IRANMR = 97
 
516
        JRANMR = 33
 
517
        END
 
518
 
 
519
 
 
520
 
 
521
 
 
522
 
 
523
 
 
524