C23456789012345678901234567890123456789012345678901234567890123456789012
C
C                  ----------------------------
C                  beowulfsubs.f  Version 1.1.1
C                  ----------------------------
C
C23456789012345678901234567890123456789012345678901234567890123456789012
      SUBROUTINE VERSION
      INTEGER NIN,NOUT
      COMMON /IOUNIT/ NIN,NOUT
      WRITE(NOUT,*)'beowulf 1.1.1 subroutines 16 August 2001'
      RETURN
      END
C23456789012345678901234567890123456789012345678901234567890123456789012
C
C Subroutines for numerical integration of jet cross sections in
C electron-positron annihilation. -- D. E. Soper
C
C First code: 29 November 1992.
C Latest revision: see Version subroutine above.
C Special note: modified 8 December 1999 to change tabs to spaces
C and to correct the header in line 4 above.
C
C The main program and subroutines that a user might want to modify
C are contained in the companion package, beowulf.f. In particular, a
C user may very well want to modify parameter settings in the main
C program and to change the observables calculated in the subroutine
C HROTHGAR and in the functions CALStype(NCUT,KCUT,index). Subroutines
C that can be modified only at extreme peril to the reliability of
C the results are in this package, beowulfsubs.f.
C
C There are two parallel calculations. Program beowulf calculates a
C sample integral, which by default is the average value of 
C (1 - thrust)^2. These are summed in the variable INTEGRAL and
C reported upon completion of the program. The program also computes
C a simple check integral in order to check on the jacobians etc.
C In the meantime, for each point in loop space and each final
C state cut, the program reports the corresponding point in the space
C of final state momenta along with the corresponding weight (Feynman
C diagram times jacobian factors) to the subroutine HROTHGAR, which
C multiplies by the measurement functions CALS corresponding to the
C measurements desired and accumulates the results.
C 
C In order to control roundoff errors, a point in loop space is rejected
C if the point is too near a singularity or if there is too much
C cancellation in the contribution from that point to INTEGRAL.
C 
C
C23456789012345678901234567890123456789012345678901234567890123456789012
C
C                        PROGRAM STRUCTURE
C
C * denotes routines found in beowulf.f. 
C Other routines are in beowulfsubs.f
C
C PROGRAM BEOWULF (*) <setup and final reporting>
C    SUBROUTINE MAKECUTINFO
C       SUBROUTINE NEWGRAPH
C          SUBROUTINE NEWCHOICE
C          SUBROUTINE NEXTCHOICE
C          FUNCTION ONEPI
C          SUBROUTINE CHECK
C             SUBROUTINE CHECKOUT
C             SUBROUTINE EXCHANGE
C       SUBROUTINE NEWCUT
C       SUBROUTINE FINDTYPES
C          SUBROUTINE NEWCUT
C    SUBROUTINE GETCOUNTS (*)
C    SUBROUTINE HROTHGAR (*) <see below>
C    SUBROUTINE RANDOMINIT
C       SUBROUTINE NEWRAN
C    SUBROUTINE DAYTIME (*)
C    SUBROUTINE VERSION
C    SUBROUTINE RENO <see below>
C    SUBROUTINE DIAGNOSTIC
C       SUBROUTINE NEWGRAPH ...
C       SUBROUTINE FINDTYPES ...
C       SUBROUTINE CHECKPOINT
C       SUBROUTINE CALCULATE <see below>
C end
C
C SUBROUTINE RENO <choses graphs and points>
C    SUBROUTINE TIMING (*)
C    SUBROUTINE HROTHGAR (*) <see below>
C    SUBROUTINE NEWGRAPH ...
C    SUBROUTINE FINDTYPES ...
C    SUBROUTINE FINDA <see below>
C       FUNCTION PROPSIGN
C    SUBROUTINE NEWPOINT
C       FUNCTION RANDOM
C          SUBROUTINE NEWRAN
C       SUBROUTINE CHOOSE3
C       SUBROUTINE CHOOSE2TO3D
C       SUBROUTINE CHOOSE2TO3E
C       SUBROUTINE CHOOSE2TO2T
C       SUBROUTINE CHOOSE2TO2S
C       SUBROUTINE CHOOSE2TO1
C    SUBROUTINE CHECKPOINT
C    SUBROUTINE CALCULATE <see below>
C return
C
C SUBROUTINE CALCULATE <the main calculation routines.
C    FUNCTION DENSITY
C       FUNCTION RHO3
C       FUNCTION RHO2TO3D
C       FUNCTION RHO2TO3E
C       FUNCTION RHO2TO2T
C       FUNCTION RHO2TO2S
C       FUNCTION RHO2TO1
C    FUNCTION CALS0 (*)
C       FUNCTION THRUST (*)
C    SUBROUTINE GETCUTINFO
C    SUBROUTINE DEFORM
C    SUBROUTINE FINDA <see below>
C    SUBROUTINE CHECKDEFORM
C    FUNCTION RNUMERATOR
C       SUBROUTINE QPROPR
C       SUBROUTINE GPROPR
C       SUBROUTINE QPROP
C          SUBROUTINE CHECKDEFORM2
C       FUNCTION RQQP3AQ
C       FUNCTION RQQP3BQ
C       FUNCTION RQQP3CQ
C       FUNCTION RQQG3AG
C       FUNCTION RQQG3BG
C       FUNCTION RQQG3CG
C       FUNCTION RQQG3AQ
C       FUNCTION RQQG3BQ
C       FUNCTION RQQG3CQ
C    FUNCTION NUMERATOR
C       SUBROUTINE QPROP
C          SUBROUTINE CHECKDEFORM2
C       SUBROUTINE GPROP
C    FUNCTION SMEAR
C    SUBROUTINE CHECKCALC
C    SUBROUTINE HROTHGAR (*) <see below>
C return
C 
C SUBROUTINE FINDA <finds matrix to get all k's from the loop l's>
C    FUNCTION PROPSIGN
C return
C 
C SUBROUTINE HROTHGAR (*) 
C  Beowulf serves Hothgar, who accepts the points in the space of final
C  state momenta with the corresponding weights, multiplies by
C  desired measurement functions CALS, and accumulates results.
C    FUNCTION CALSTHRUST (*)
C       FUNCTION THRUST (*)
C       FUNCTION KN (*)
C    FUNCTION CALS3JET (*)
C       SUBROUTINE COMBINEJETS (*)
C       FUNCTION BETHKE (*)
C    FUNCTION KN (*)
C    FUNCTION BETHKE (*)
C return
C
C Simple functions called from routines above, with calls
C not listed above:
C 
C SUBROUTINE AXES
C FUNCTION XXREAL
C FUNCTION XXIMAG
C FUNCTION COMPLEXSQRT
C FUNCTION FACTORIAL
C FUNCTION SINHINV
C FUNCTION DELTA
C 
C23456789012345678901234567890123456789012345678901234567890123456789012
C
C A brief introduction to the variables used:
C
C Size of the calculation:
C NLOOPS = number of loops (in cut photon self energy graph).
C NPROPS = number of propagators in graph, = 3 * NLOOPS - 1.
C NVERTS = number of vertices in graph, = 2 * NLOOPS.
C CUTMAX = NLOOPS + 1
C        = maximum number of cut propagators;
C        = number of independent loop momenta needed to determine the
C          propagator momenta, counting the virtual photon momentum.
C The current program is restricted to 0 and 1 virtual loops.
C
C Labels:
C L    = index of loop momenta, L = 0,1,...,NLOOPS.
C        L = 0 normally denontes the virtual photon momentum.
C P    = index of propagator, P = 0,1,...,NPROPS.
C        P = 0 denotes the virtual photon momentum.
C Q(L) = index P of propagator carrying the Lth loop momentum.
C V    = index of vertices, V = 1,...,NVERTS
C
C Momentum variables (MU = 0,1,2,3):
C K(P,MU)       = Momentum of Pth propagator.
C                 For P = 0, this is the virtual photon momentum:
C                 K(0,MU) = 0 for MU = 1, 2, 3 while K(0,0) = RTS.
C ABSK(P)       = Square of the three momentum of Pth propagator.
C KINLOOP(J,MU) = K(LOOPINDEX(J),MU) = momenta of loop propagators.
C KCUT(I,MU)    = K(CUTINDEX(I),MU) = momenta of cut propagators.
C K(Q(L),MU)    = Lth loop momentum, L = 0,...,NLOOPS;
C KC(P,MU)      = complex propagator momenta.
C A(P,L)        = Matrix relating propagator momenta to loop momenta.
C                 K(P,MU) = SUM_{L=0}^{NLOOPS} A(P,L) K(Q(L),MU)
C
C Variables from NEWGRAPH:
C VRTX(P,I) = Index of vertex at beginning (i= 1) and end (I = 2) of
C              of propagator P. Specifies the supergraph.
C PROP(V,I) = Index of Ith propagator attached to vertex V, I = 1,2,3.
C              Also specifies the supergraph.
C SELFPROP(P) = True if propagator P is part of a one loop self-energy
C               subgraph or attaches to a such a subgraph.
C
C Variables associated with NEWPOINT and FINDTYPES:
C NMAPS = Number of different maps from random x's to momenta.
C MAPNMUMBER = Number labelling a certain map.
C QS(MAPNUMBER,II) = Label of the IIth propagator that is special
C                   in map number MAPNUMBER.
C QSIGNS(MAPNUMBER,II) = sign needed to relate the conventional 
C       direction of the propagator to that in an elementary scattering
C MAPTYPES(MAPNUMBER) = T2TO3, T2TO2T, T2TO2S, T2TO1.
C
C JACNEWPOINT =1/DENSITY(GRAPHNUMBER,K,QS,QSIGNS,MAPTYPES,NMAPS)
C             = Jacobian for loop momenta L.
C
C Variables from NEWCUT:
C NEWCUTINIT: .TRUE. tells NEWCUT to initialize itself.
C NCUT = Number of cut propagators.
C ISIGN(P) = +1 if propagator P is left of cut, -1 if right, 0 if cut.
C CUTINDEX(I) = Index P of cut propagator I, I = 1,...,CUTMAX.
C CUTSIGN(I) = Sign of cut propagator I I = 1,...,CUTMAX.
C              (+1 if K(P,0) >0 for cut propagator.)
C LEFTLOOP = True iff there is a virtual loop to the left of the cut.
C RIGHTLOOP = True iff there is a virtual loop to the right of the cut.
C NINLOOP = Number of propagators in loop.
C LOOPINDEX(NP) = Index P of NPth propagator around the loop.
C LOOPSIGN(NP) = 1 if propagator direction is same as loop direction.
C               -1 if  direction is opposite to loop direction.
C NP = JCUT: Propagator cut by loopcut.
C CUTFOUND: .TRUE. if NEWCUT found a new cut.
C
C In RENO we use CUTINDEX to define CUT(P) = True if propagator
C    P is cut.
C
C Solving for the propagator energies:
C For NCUT = CUTMAX, cut propagators are P = CUTINDEX(I).
C   with direction of positive energy given by CUTSIGN(I).
C For NCUT = CUTMAX - 1, we define a "loopcut" on the propagator
C   numbered JCUT in order around the loop, 1.LE.JCUT.LE.NINLOOP:
C         CUTINDEX(CUTMAX) = LOOPINDEX(JCUT) and
C         CUTSIGN(CUTMAX) = LOOPSIGN(JCUT).
C Energies of cut propagators are
C   E(I-1) = K(CUTINDEX(I),0) for I = 1,...,CUTMAX.
C and are determined from
C   E(I-1) = CUTSIGN(I) * SQRT( Sum_J [  K(CUTINDEX(I),J)**2 ] ).
C This gives energies E(L) for L = 0,...,NLOOPS. We consider the
C propagators designated by QE(L) = CUTINDEX(L+1) as independent
C and generate the matrix AE(P,L) that gives the propagator energies
C in terms of these independent momenta. This gives the propagator
C energies.
C
C Contour deformation:
C   NEWKINLOOP(MU) =  addition to the momentum going around the loop
C                     caused by deforming the contour. We have
C   Im[ KC(LOOPINDEX(J,MU)) ] = LOOPSIGN(LOOPINDEX(J)) 
C                              * Im[ NEWKINLOOP(J,MU) ] for MU = 1,2,3.
C
C23456789012345678901234567890123456789012345678901234567890123456789012
C23456789012345678901234567890123456789012345678901234567890123456789012
C
      SUBROUTINE RENO(
     >          SUMR,ERRORR,SUMI,ERRORI,
     >          SUMBIS,ERRORBIS,
     >          SUMCHKR,ERRORCHKR,SUMCHKI,ERRORCHKI,FLUCT,
     >          INCLUDED,EXTRACOUNT,OMITTED,
     >          NVALPT,VALPTMAX,KBAD,BADGRAPHNUMBER,BADMAPNUMBER,
     >          NRENO,CPUTIME)
C Array sizes:
      INTEGER SIZE,MAXGRAPHS,MAXMAPS
      PARAMETER (SIZE = 3)
      PARAMETER (MAXGRAPHS = 10)
      PARAMETER (MAXMAPS = 64)
C Out:
      REAL*8 SUMR,ERRORR,SUMI,ERRORI
      REAL*8 SUMBIS,ERRORBIS
      REAL*8 SUMCHKR,ERRORCHKR,SUMCHKI,ERRORCHKI
      REAL*8 FLUCT(MAXGRAPHS,MAXMAPS)
      INTEGER*8 INCLUDED,EXTRACOUNT,OMITTED
      INTEGER NVALPT(-9:6)
      REAL*8 VALPTMAX
      REAL*8 KBAD(0:3*SIZE-1,0:3)
      INTEGER BADGRAPHNUMBER,BADMAPNUMBER
      INTEGER NRENO
      REAL*8 CPUTIME
C
C Computes the cross section integral by Monte Carlo integration.
C
C Latest revision 11 August 1999
C
      INTEGER NIN,NOUT
      COMMON /IOUNIT/ NIN,NOUT
      INTEGER NLOOPS,NPROPS,NVERTS,CUTMAX
      COMMON /SIZES/ NLOOPS,NPROPS,NVERTS,CUTMAX
      REAL*8 ENERGYSCALE
      COMMON /MSCALE/ ENERGYSCALE
      REAL*8 BADNESSLIMIT,CANCELLIMIT,THRUSTCUT
      COMMON /LIMITS/ BADNESSLIMIT,CANCELLIMIT,THRUSTCUT
      REAL*8 TIMELIMIT
      COMMON /MAXTIME/ TIMELIMIT
C Graphs to include
      LOGICAL USEGRAPH(MAXGRAPHS)
      COMMON /WHICHGRAPHS/ USEGRAPH
C How many graphs and how many cuts and maps for each:
      INTEGER NUMBEROFGRAPHS
      INTEGER NUMBEROFCUTS(MAXGRAPHS)
      INTEGER NUMBEROFMAPS(MAXGRAPHS)
      COMMON /GRAPHCOUNTS/ NUMBEROFGRAPHS,NUMBEROFCUTS,NUMBEROFMAPS
C Momenta:
      REAL*8 K(0:3*SIZE-1,0:3),ABSK(0:3*SIZE-1)
C Matrices:
      INTEGER A(0:3*SIZE-1,0:SIZE)
C NEWGRAPH variables:
      INTEGER VRTX(0:3*SIZE-1,2),PROP(2*SIZE,3)
      LOGICAL SELFPROP(3*SIZE-1)
      LOGICAL GRAPHFOUND
      INTEGER GRAPHNUMBER
C FINDA variable:
      LOGICAL QOK
C MAP variables:
      INTEGER NMAPS,MAPNUMBER
      INTEGER QS(MAXMAPS,0:SIZE),QSIGNS(MAXMAPS,0:SIZE)
      INTEGER Q(0:SIZE),QSIGN(0:SIZE)
      CHARACTER*6 MAPTYPES(MAXMAPS)
      CHARACTER*6 MAPTYPE
C Variable from CHECKPOINT:
      REAL*8 BADNESS
C Problem report from NEWPOINT
      LOGICAL BADNEWPOINT
C Logical variables to tell how to treat point:
      LOGICAL XTRAPOINTQ, BADPOINTQ
C Functions:
      REAL*8 XXREAL,XXIMAG
C Index variables:
      INTEGER L,P,MU
C Hrothgar dummy variables:
      REAL*8 KCUT0(SIZE+1,0:3)
C Reno size and counting variables:
      INTEGER GROUPSIZE(MAXGRAPHS,MAXMAPS)
      INTEGER GROUPSIZEGRAPH(MAXGRAPHS)
      INTEGER GROUPSIZETOTAL
      COMMON /MONTECARLO/GROUPSIZE,GROUPSIZEGRAPH,GROUPSIZETOTAL
      INTEGER POINT
C Reno results variables:
      REAL*8 SQRSUMR,SQRSUMCHKR
      REAL*8 SQRSUMI,SQRSUMCHKI
      REAL*8 SQRSUMBIS
      COMPLEX*16 INTEGRAL,INTEGRALCHK
      REAL*8 INTEGRALBIS
C Calculate variables:
      COMPLEX*16 VALUE,VALUECHK
      REAL*8 MAXPART
      REAL*8 VALPT,LOGVALPT
      LOGICAL REPORT,DETAILS
      COMMON /CALCULOOK/ REPORT,DETAILS
C Timing variables
      REAL*8 DELTATIME
C
C------------------------------ Begin ----------------------------------
C
C Dummy variables for Hrothgar.
C
      DO L = 1,SIZE+1
        DO MU = 0,3
         KCUT0(L,MU) = 1.0D0
        ENDDO
      ENDDO
C
C Initialize CPUTIME and NRENO. Call to TIMING starts the clock.
C
      CPUTIME = 0.0
      NRENO = 0
      CALL TIMING(DELTATIME)
C
C Initialize sums for loop over groups of Reno points. The sums
C will be updated for each group. Within a group, the quantities
C corresponding to SUMxxR + i SUMxxI are complex variables called
C INTEGRALxx.
C
      SUMR = 0.0D0
      SUMI = 0.0D0
      SUMBIS = 0.0D0
      SUMCHKR = 0.0D0
      SUMCHKI = 0.0D0
C
      SQRSUMR = 0.0D0
      SQRSUMI = 0.0D0
      SQRSUMBIS = 0.0D0
      SQRSUMCHKR = 0.0D0
      SQRSUMCHKI = 0.0D0
C
      DO GRAPHNUMBER = 1,NUMBEROFGRAPHS
      DO MAPNUMBER = 1,NUMBEROFMAPS(GRAPHNUMBER)
        FLUCT(GRAPHNUMBER,MAPNUMBER) = 0.0D0
      ENDDO
      ENDDO
C
      DO L = -9,6
        NVALPT(L) = 0
      ENDDO
      VALPTMAX = 0.0D0
      INCLUDED = 0
      EXTRACOUNT = 0
      OMITTED = 0
C
C Tell CALCULATE not to report its findings for each calculation
C
      REPORT = .FALSE.
C
C Initialize integrals for first group.
C
      INTEGRAL = (0.0D0,0.0D0)
      INTEGRALBIS = 0.0D0
      INTEGRALCHK = (0.0D0,0.0D0)
C
C Loop over groups of points.
C
      DO WHILE (CPUTIME.LT.TIMELIMIT)
      NRENO = NRENO + 1
C
C Call Hrothgar to tell him to that we are starting a new group.
C
      CALL HROTHGAR(1,KCUT0,1.0D0,1,'STARTGROUP')
C
C Get a new graph.
C
      GRAPHFOUND = .TRUE.
      GRAPHNUMBER = 0
C
      DO WHILE (GRAPHFOUND)
      CALL NEWGRAPH(VRTX,PROP,SELFPROP,GRAPHFOUND)
      IF (GRAPHFOUND) THEN
      GRAPHNUMBER = GRAPHNUMBER + 1
C
C Calculate number of maps NMAPS, index arrays QS,
C types MAPTYPES, and signs QSIGNS associated with the maps.
C
      CALL FINDTYPES(VRTX,PROP,NMAPS,QS,QSIGNS,MAPTYPES)
C
C Check if we were supposed to use this graph (USEGRAPH is
C set in the main program.)
C
      IF (USEGRAPH(GRAPHNUMBER)) THEN
C
C Loop over choices of maps from x's to loop momenta.
C
      DO MAPNUMBER = 1,NMAPS
C
      MAPTYPE = MAPTYPES(MAPNUMBER)
      DO L = 0,NLOOPS
        Q(L) = QS(MAPNUMBER,L)
        QSIGN(L) = QSIGNS(MAPNUMBER,L)
      ENDDO
C
      CALL FINDA(VRTX,Q,NLOOPS,A,QOK)
C
C Loop over Reno points within a group.
C
      DO POINT = 1,GROUPSIZE(GRAPHNUMBER,MAPNUMBER)
C
C Call Hrothgar to tell him that we are starting a new point.
C
      CALL HROTHGAR(1,KCUT0,1.0D0,1,'STARTPOINT')
C
C Get a new point. Check on its badness. If it is too bad,
C or if NEWPOINT reported a problem, we omit the point after
C notifying Hrothgar.
C
      BADPOINTQ = .FALSE.
      XTRAPOINTQ = .FALSE.
      CALL NEWPOINT(A,QSIGN,MAPTYPE,K,ABSK,BADNEWPOINT)
      IF (BADNEWPOINT) THEN
        CALL HROTHGAR(1,KCUT0,1.0D0,1,'BADPOINT  ')
        BADPOINTQ = .TRUE.
      ENDIF
      CALL CHECKPOINT(K,ABSK,PROP,BADNESS)
      IF (BADNESS.GT.100*BADNESSLIMIT) THEN
        CALL HROTHGAR(1,KCUT0,1.0D0,1,'BADPOINT  ')
        BADPOINTQ = .TRUE.
      ELSE IF (BADNESS.GT.BADNESSLIMIT) THEN
        CALL HROTHGAR(1,KCUT0,1.0D0,1,'XTRAPOINT ')
        XTRAPOINTQ = .TRUE.
      ENDIF
C
C If the point is not too bad, we can call CALCULATE.
C The final state momenta found, KCUT, along with the corresponding
C weights, are reported to Hrothgar by CACULATE.
C Then call Hrothgar to tell him that we are done with this point.
C
      IF (.NOT.BADPOINTQ) THEN
      CALL CALCULATE(VRTX,SELFPROP,GRAPHNUMBER,K,ABSK,
     >               QS,QSIGNS,MAPTYPES,NMAPS,VALUE,MAXPART,VALUECHK)
      ENDIF
C
C Add contribution from this point to integral.
C We count the point if Maxvalue/|Value| < Cancellimit.
C
      IF (.NOT.BADPOINTQ) THEN
        IF ( MAXPART.GT. 100*CANCELLIMIT*ABS(XXREAL(VALUE)) ) THEN
          CALL HROTHGAR(1,KCUT0,1.0D0,1,'BADPOINT  ')
          BADPOINTQ = .TRUE.
        ELSE IF ( MAXPART.GT. CANCELLIMIT*ABS(XXREAL(VALUE)) ) THEN
          CALL HROTHGAR(1,KCUT0,1.0D0,1,'XTRAPOINT ')
          XTRAPOINTQ = .TRUE.
        ENDIF
      ENDIF
C      
      IF ( (.NOT.BADPOINTQ).AND.(.NOT.XTRAPOINTQ) ) THEN
        INTEGRAL = INTEGRAL + VALUE
        FLUCT(GRAPHNUMBER,MAPNUMBER) = FLUCT(GRAPHNUMBER,MAPNUMBER)
     >    + XXREAL(VALUE)**2/GROUPSIZE(GRAPHNUMBER,MAPNUMBER)
        INTEGRALCHK = INTEGRALCHK + VALUECHK
        INCLUDED = INCLUDED + 1
C
C For diagnostic purposes, we need VALPT, the contribution to
C the integral being calculated from this point, normalized such
C that the integral is the sum over all points chosen of VALPT
C divided by the total number of points, NRENO * GROUPSIZETOTAL.
C
        VALPT = ABS(XXREAL(VALUE))*GROUPSIZETOTAL
        LOGVALPT = LOG10(VALPT)
        DO L = -9,6
          IF((LOGVALPT.GE.L).AND.(LOGVALPT.LT.(L+1))) THEN
            NVALPT(L) = NVALPT(L) + 1
          ENDIF
        ENDDO
        IF (VALPT.GT.VALPTMAX) THEN
          VALPTMAX = VALPT
          DO P = 1,NPROPS
          DO MU = 1,3
            KBAD(P,MU) = K(P,MU)
          ENDDO
          ENDDO
          BADGRAPHNUMBER = GRAPHNUMBER
          BADMAPNUMBER = MAPNUMBER
        ENDIF
      ELSE IF ((.NOT.BADPOINTQ).AND.(XTRAPOINTQ) ) THEN
C
C For points that are 'extra', we include the value of
C the integrand in the INTEGRALBIS, which will provide an estimate
C or the effect of the cutoffs.
C
        INTEGRALBIS = INTEGRALBIS + XXREAL(VALUE)
        EXTRACOUNT = EXTRACOUNT + 1
C
      ELSE
        OMITTED = OMITTED + 1
      ENDIF
C
C End of loop over POINT.
C
      CALL HROTHGAR(1,KCUT0,1.0D0,1,'POINTDONE ')
      ENDDO
C
C End of loop over MAPNUMBER.
C
      ENDDO
C
C End for IF (USEGRAPH(GRAPHNUMBER)) THEN
C
      ENDIF
C
C End of loop DO WHILE (GRAPHFOUND)/ IF (GRAPHFOUND).
C
      ENDIF
      ENDDO
C
C Call Hrothgar to tell him that we are done with this group.
C
      CALL HROTHGAR(1,KCUT0,1.0D0,1,'GROUPDONE ')
C
C Add results from this group to the SUM variables. 
C
      SUMR = SUMR + XXREAL(INTEGRAL)
      SUMI = SUMI + XXIMAG(INTEGRAL)
      SUMBIS = SUMBIS + INTEGRALBIS
      SUMCHKR = SUMCHKR + XXREAL(INTEGRALCHK)
      SUMCHKI = SUMCHKI + XXIMAG(INTEGRALCHK)
C
      SQRSUMR = SQRSUMR + XXREAL(INTEGRAL)**2
      SQRSUMI = SQRSUMI + XXIMAG(INTEGRAL)**2
      SQRSUMBIS = SQRSUMBIS + INTEGRALBIS**2
      SQRSUMCHKR = SQRSUMCHKR + XXREAL(INTEGRALCHK)**2
      SQRSUMCHKI = SQRSUMCHKI + XXIMAG(INTEGRALCHK)**2
C
C Reset the INTEGRAL variables for the next group.
C
      INTEGRAL = (0.0D0,0.0D0)
      INTEGRALBIS = 0.0D0
      INTEGRALCHK = (0.0D0,0.0D0)
C
C End of loop DO WHILE (CPUTIME.LT.TIMELIMIT)
C
      CALL TIMING(DELTATIME)
      CPUTIME = CPUTIME + DELTATIME
      ENDDO
C
C Calculate the SUM results.
C
      SUMR = SUMR/NRENO
      SUMI = SUMI/NRENO
      SUMBIS = SUMBIS/NRENO
      SUMCHKR = SUMCHKR/NRENO
      SUMCHKI = SUMCHKI/NRENO
C
      SQRSUMR = SQRSUMR/NRENO
      SQRSUMI = SQRSUMI/NRENO
      SQRSUMBIS = SQRSUMBIS/NRENO
      SQRSUMCHKR = SQRSUMCHKR/NRENO
      SQRSUMCHKI = SQRSUMCHKI/NRENO
C
      IF (NRENO.EQ.1) THEN
        WRITE(NOUT,*)'NRENO = 1  changed to 2 to avoid 1/0'
        WRITE(NOUT,*)' '
        NRENO = 2
      ENDIF
      ERRORR = SQRT((SQRSUMR - SUMR**2)/(NRENO - 1))
      ERRORI = SQRT((SQRSUMI - SUMI**2)/(NRENO - 1))
      ERRORBIS = SQRT((SQRSUMBIS - SUMBIS**2)/(NRENO - 1))
      ERRORCHKR = SQRT((SQRSUMCHKR - SUMCHKR**2)/(NRENO-1))
      ERRORCHKI = SQRT((SQRSUMCHKI - SUMCHKI**2)/(NRENO-1))
C
      DO GRAPHNUMBER = 1,NUMBEROFGRAPHS
      DO MAPNUMBER = 1,NUMBEROFMAPS(GRAPHNUMBER)
        FLUCT(GRAPHNUMBER,MAPNUMBER) =
     >         FLUCT(GRAPHNUMBER,MAPNUMBER)/NRENO
      ENDDO
      ENDDO
C
      RETURN
      END
C
C23456789012345678901234567890123456789012345678901234567890123456789012
C                 Subroutines associated with NEWGRAPH                 C
C23456789012345678901234567890123456789012345678901234567890123456789012
C
      SUBROUTINE NEWGRAPH(VRTX,PROP,SELFPROP,GRAPHFOUND)
C
      INTEGER SIZE
      PARAMETER (SIZE = 3)
C Out:
      INTEGER VRTX(0:3*SIZE-1,2),PROP(2*SIZE,3)
      LOGICAL SELFPROP(3*SIZE-1)
      LOGICAL GRAPHFOUND
C
C 8 November 1992 Home fixup of bugs.
C 28 November 1992 Add check that we get each graph only once.
C 13 July 1994
C 13 April 1996
C  1 January 1998 Add output variable SELFPROP. Omit NPERMS as output.
C----------
C Varibles:
C VRTX(P,I) = Index of vertex at beginning (i= 1) and end (I = 2) of
C             of propagator P. Specifies the supergraph for output.
C PROP(V,I) = Index of Ith propagator attached to vertex V, I = 1,2,3.
C              Also specifies the supergraph for output.
C SELFPROP(P) = True if propagator P is part of a one loop self-energy
C               subgraph or attaches to a such a subgraph.
C C(V,I) = Index of Ith vertex connected to vertex V.
C        V = 1,...,NVERTS; I =1,2,3; C(V,I) = 1,...,NVERTS and -1,-2.
C        Here C(V,1).LE.C(V,2).LE.C(V,3).
C        This is the fundamental specification of the supergraph.
C N = Number of permutations of the vertices that give same graph.
C GRAPHFOUND = True when the subroutine finds a new graph.
C COUNT(V) = Number of vertices connected to vertex V.
C Vertex 1 is automatically connected to the photon "-1":C(1,1) = -1.
C Vertex 2 is automatically connected to the photon "-2":C(2,1) = -2.
C The freedom to renumber the vertices 3,...,NVERTS is used to choose
C a standard numbering:
C  We choose the numbering with the smallest value of C(1,1);
C  For numberings with equal values of C(1,1) we choose the numbering
C    with the smallest value of C(1,2);
C  For numberings with equal values of C(1,2) we choose the numbering
C    with the smallest value of C(1,3);
C  For numberings with equal values of C(1,3) we choose the numbering
C    with the smallest value of C(2,1); et cetera.
C
C The connections are generated starting with vertex 1.  We make
C a choice of connections for vertex V, then move on to make a choice
C for connections to vertex V + 1.  When we are out of choices for
C connections to vertex V, we step back and try the next choice for
C vertex V - 1.
C
C Connections to the external boson:
C In C(V,I) we assign the first connection of vertex 1 to be vertex "-1"
C while the first connection of vertex 2 is vertex "-2." This numbering
C is convenient for working out C(V,I). In reporting the results, 
C however, we label the external boson with propagator 0, so that
C PROP(1,1) = PROP(2,1) = 0. Then propagator 0 attaches to vertices
C 1 and 2: VERT(0,1) = 2, VERT(0,2) = 1.
C----------
C
      INTEGER NIN,NOUT
      COMMON /IOUNIT/ NIN,NOUT
      INTEGER NLOOPS,NPROPS,NVERTS,CUTMAX
      COMMON /SIZES/ NLOOPS,NPROPS,NVERTS,CUTMAX
C
      INTEGER C(2*SIZE,3),COUNT(2*SIZE)
      INTEGER NUSED(2*SIZE),VA,VB
      INTEGER V,VV,I,P,NPERMS
      LOGICAL ONEPI,OK
      LOGICAL FAIL,NEWSTART,UP
      DATA NEWSTART/.TRUE./
      SAVE
C
C Initializations.
C
      IF (NEWSTART) THEN
        DO  VV = 1,NVERTS
          COUNT(VV) = 0
          DO I = 1,3
            C(VV,I) = 0
          ENDDO
        ENDDO
        C(1,1) = -1
        COUNT(1) = 1
        C(2,1) = -2
        COUNT(2) = 1
        V = 1
        UP = .TRUE.
      ENDIF
C
C Move from level to level in tree structure of choices. When UP
C is true, we have moved to a higher V; when UP is false, we have
C moved to a smaller V.
C
      DO WHILE (.TRUE.)
C
      IF (UP) THEN
        CALL NEWCHOICE(C,COUNT,V,FAIL)
      ELSE
        CALL NEXTCHOICE(C,COUNT,V,FAIL)
      ENDIF
      IF (FAIL) THEN
C
C If we couldn't find connectections for vertex V, then we should
C step back and look for the next connections for vertex V-1.  But if
C V is currently 1, then we can't step back, so we have found all
C the graphs.
C
        IF (V.GT.1) THEN
          V = V - 1
          UP = .FALSE.
        ELSE
          NEWSTART = .TRUE.
          GRAPHFOUND = .FALSE.
          DO P = 0,NPROPS
           DO I = 1,2
            VRTX(P,I) = 0
           ENDDO
          ENDDO
          RETURN
        ENDIF
C
C If we did find connections for vertex V, then we should step onward
C and look for new connections for vertex V+1.  But if V is currently
C equal to NVERTS, then we must have found a graph.  We check for
C validity. If it is valid, we exit with the results, setting V and UP
C so that the next time the subroutine is called we will start looking
C for the next connections for vertex V-1.  If our graph is not valid
C (eg. one particle reducible) then we step back to look for new
C connections for vertex V-1 right away.
C
      ELSE
        IF (V.LT.NVERTS) THEN
          V = V + 1
          UP = .TRUE.
        ELSE
          V = V - 1
          UP = .FALSE.
          IF (ONEPI(C)) THEN
          CALL CHECK(C,NPERMS,OK)
            IF (OK) THEN
              NEWSTART = .FALSE.
              GRAPHFOUND = .TRUE.
C
C Exit. We translate the results for C(V,I) into VRTX(P,I), I = 1,2,
C and PROP(P,I), I = 1,2,3. Here NUSED(V) denotes how many propagators
C we have so far assigned connecting to vertex V.
C
              DO VV = 1,NVERTS
                NUSED(VV) = 0
              ENDDO
              VRTX(0,1) = 2
              VRTX(0,2) = 1
              PROP(1,1) = 0
              NUSED(1)  = 1
              PROP(2,1) = 0
              NUSED(2)  = 1
              P = 1
              DO VV = 1,NVERTS
              DO I = 1,3
                IF (C(VV,I).GT.VV) THEN
                  VA = VV
                  VB = C(VV,I)
                  VRTX(P,1) = VA
                  NUSED(VA) = NUSED(VA) + 1
                  PROP(VA,NUSED(VA)) = P
                  VRTX(P,2) = VB
                  NUSED(VB) = NUSED(VB) + 1
                  PROP(VB,NUSED(VB)) = P
                  P = P+1
                ENDIF
              ENDDO
              ENDDO
              IF (P.NE.NPROPS+1) THEN
                WRITE(NOUT,*)'SNAFU in NEWGRAPH',P-1,NPROPS
                STOP
              ENDIF
              DO VV = 1,NVERTS
                IF (NUSED(VV).NE.3) THEN
                  WRITE(NOUT,*)'Problem in NEWWGRAPH',VV,NUSED(VV)
                  STOP
                ENDIF
              ENDDO
C
C We also need to report which propagators P connect to a vertex
C that is part of a one loop self-energy subgraph.
C
              DO P = 1,NPROPS
                SELFPROP(P) = .FALSE.
              ENDDO
              DO VV = 1,NVERTS
                IF ((C(VV,1).EQ.C(VV,2)).OR.(C(VV,2).EQ.C(VV,3))) THEN
                  DO I = 1,3
                    SELFPROP(PROP(VV,I)) = .TRUE.
                  ENDDO
                ENDIF
              ENDDO
C
C OK. We are ready to return.
C
              RETURN
C
            ENDIF
          ENDIF
        ENDIF
      ENDIF
C
C End main loop "DO WHILE (.TRUE.)"
C
      ENDDO
C
      END
C
C
C23456789012345678901234567890123456789012345678901234567890123456789012
C23456789012345678901234567890123456789012345678901234567890123456789012
C
      SUBROUTINE NEWCHOICE(C,COUNT,V,FAIL)
C
      INTEGER SIZE
      PARAMETER (SIZE = 3)
C
      INTEGER C(2*SIZE,3),COUNT(2*SIZE)
      INTEGER V
      LOGICAL FAIL
C
      INTEGER NIN,NOUT
      COMMON /IOUNIT/ NIN,NOUT
      INTEGER NLOOPS,NPROPS,NVERTS,CUTMAX
      COMMON /SIZES/ NLOOPS,NPROPS,NVERTS,CUTMAX
C
      INTEGER VV,K
      LOGICAL FOUND
      SAVE
C
C If COUNT(V) = 3, then we don't need any more connections to this
C vertex. If COUNT(V) = 0, then we appear to be starting to make a
C vacuum graph after having completed a graph with too few loops, so
C we should just quit.
C
      IF (COUNT(V).EQ.3) THEN
          FAIL = .FALSE.
          RETURN
      ELSE IF (COUNT(V).EQ.0) THEN
          FAIL = .TRUE.
          RETURN
      ENDIF
C
C Generate starting choice for new vertices to connect to V. We connect
C to the vertices with the smallest possible indices.
C
      FAIL = .FALSE.
      VV = V + 1
      DO K = (COUNT(V) + 1),3
        FOUND = .FALSE.
        DO WHILE (.NOT.FOUND)
          IF (VV.GT.NVERTS) THEN
              FAIL = .TRUE.
              RETURN
          ENDIF
          IF ( COUNT(VV).LT.3) THEN
            COUNT(V) = COUNT(V) + 1
            C(V,K) = VV
            COUNT(VV) = COUNT(VV) + 1
            C(VV,COUNT(VV)) = V
            FOUND = .TRUE.
          ELSE
            VV = VV + 1
          ENDIF
        ENDDO
      ENDDO
C
      RETURN
      END
C
C23456789012345678901234567890123456789012345678901234567890123456789012
C
      SUBROUTINE NEXTCHOICE(C,COUNT,V,FAIL)
C
      INTEGER SIZE
      PARAMETER (SIZE = 3)
C
      INTEGER C(2*SIZE,3),COUNT(2*SIZE)
      INTEGER V
      LOGICAL FAIL
C
      INTEGER NIN,NOUT
      COMMON /IOUNIT/ NIN,NOUT
      INTEGER NLOOPS,NPROPS,NVERTS,CUTMAX
      COMMON /SIZES/ NLOOPS,NPROPS,NVERTS,CUTMAX
C
      INTEGER VV,VVV,V2,V3,I
      LOGICAL FOUND
      SAVE
C
C First, erase any connections among higher index vertices.
C
      DO VV = V+1,NVERTS
      DO VVV = V+1,NVERTS
      DO I = 1,3
        IF (C(VV,I).EQ.VVV) THEN
          C(VV,I) = 0
          COUNT(VV) = COUNT(VV) - 1
        ENDIF
      ENDDO
      ENDDO
      ENDDO
C
C Next, get the next connection set for vertex V.
C First, we try to find a new third connection for V.
C
      V3 = C(V,3)
C If third connection was to a lower index vertex, we can't change it.
      IF (V3.LE.V) THEN
        FAIL = .TRUE.
        RETURN
      ENDIF
C Erase third connection:
      C(V,3) = 0
      C(V3,COUNT(V3)) = 0
      COUNT(V) = COUNT(V) - 1
      COUNT(V3) = COUNT(V3) - 1
C Look for a new one:
      DO WHILE (V3.LT.NVERTS)
        V3 = V3 + 1
        IF ((COUNT(V3-1).GT.0).AND.(COUNT(V3).LT.3)) THEN
          COUNT(V) = COUNT(V) + 1
          COUNT(V3) = COUNT(V3) + 1
          C(V,3) = V3
          C(V3,COUNT(V3)) = V
          FAIL = .FALSE.
          RETURN
        ENDIF
      ENDDO
C
C We have failed to find a new third connection for V, so
C try for a second connection.  
C
      V2 = C(V,2)
C If second connection was to a lower index vertex, we can't change it.
      IF (V2.LE.V) THEN
        FAIL = .TRUE.
        RETURN
      ENDIF
C Erase second connection:
      C(V,2) = 0
      C(V2,COUNT(V2)) = 0
      COUNT(V) = COUNT(V) - 1
      COUNT(V2) = COUNT(V2) - 1
C Look for a new one:
      DO WHILE (V2.LT.NVERTS)
        V2 = V2 + 1
        IF ((COUNT(V2-1).GT.0).AND.(COUNT(V2).LT.3)) THEN
          COUNT(V) = COUNT(V) + 1
          COUNT(V2) = COUNT(V2) + 1
          C(V,2) = V2
          C(V2,COUNT(V2)) = V
C We found a new second connection.  Now get a third connection.
C---
          V3 = V2
          FOUND = .FALSE.
          DO WHILE (.NOT.FOUND)
            IF ( COUNT(V3).LT.3) THEN
              COUNT(V) = COUNT(V) + 1
              COUNT(V3) = COUNT(V3) + 1
              C(V,3) = V3
              C(V3,COUNT(V3)) = V
              FOUND = .TRUE.
            ELSE
              V3 = V3 + 1
              IF (V3.GT.NVERTS) THEN
                FAIL = .TRUE.
                RETURN
              ENDIF
            ENDIF
          ENDDO
C---
C We have found a good third connection also, so we are done!
          FAIL = .FALSE.
          RETURN
        ENDIF
      ENDDO
C We couldn't find a second connection
C
      FAIL = .TRUE.
      RETURN
      END
C
C23456789012345678901234567890123456789012345678901234567890123456789012
C
      LOGICAL FUNCTION ONEPI(CIN)
C
      INTEGER SIZE
      PARAMETER (SIZE = 3)
      INTEGER CIN(2*SIZE,3)
C
C Checks that the graph is connected and 1 particle irreducible.
C Modified 26 July 1994.
C
      INTEGER NIN,NOUT
      COMMON /IOUNIT/ NIN,NOUT
      INTEGER NLOOPS,NPROPS,NVERTS,CUTMAX
      COMMON /SIZES/ NLOOPS,NPROPS,NVERTS,CUTMAX
C
      LOGICAL LEFT(2*SIZE),CHANGE
      INTEGER C(2*SIZE,3)
      INTEGER V,I,V1,V2,I1,I2
      SAVE
C
C Initialize
C
      ONEPI = .TRUE.
C
      DO V = 1,NVERTS
      DO I = 1,3
        C(V,I) = CIN(V,I)
      ENDDO
      ENDDO
C
C Set up loops to successively erase each propagator.
C
      DO V1 = 1,NVERTS
      DO I1 = 1,3
      V2 = C(V1,I1)
      IF (V2.GT.V1) THEN
      DO I2 = 1,3
      IF (C(V2,I2).EQ.V1) THEN
        C(V1,I1) = 0
        C(V2,I2) = 0
C--We have now erased the propagator from V1 to V2.  Let's see if
C  the remaining graph is connected.
        DO V = 1,NVERTS
         LEFT(V) = .FALSE.
        ENDDO
C Construct Left set.
        LEFT(1) = .TRUE.
        CHANGE = .TRUE.
        DO WHILE (CHANGE)
          CHANGE = .FALSE.
          DO V = 1,NVERTS
          DO I = 1,3
            IF ( (1.LE.C(V,I)).AND.(C(V,I).LE.NVERTS) ) THEN
             IF ( LEFT(V) .AND. (.NOT.LEFT(C(V,I))) ) THEN
              CHANGE = .TRUE.
              LEFT(C(V,I)) = .TRUE.
             ENDIF
            ENDIF
          ENDDO
          ENDDO
        ENDDO
C Check for connectedness
        DO V = 1,NVERTS
         IF ( .NOT.LEFT(V) ) THEN
            ONEPI = .FALSE.
            RETURN
         ENDIF
        ENDDO
C--OK, that remaining graph was OK.  Restore the graph.
        C(V1,I1) = V2
        C(V2,I2) = V1
      ENDIF
      ENDDO
      ENDIF
      ENDDO
      ENDDO
C
      RETURN
      END
C
C23456789012345678901234567890123456789012345678901234567890123456789012
C
      SUBROUTINE CHECK(CIN,NPERMS,OK)
C
      INTEGER SIZE
      PARAMETER (SIZE = 3)
      INTEGER CIN(2*SIZE,3),NPERMS
      LOGICAL OK
C
      INTEGER NIN,NOUT
      COMMON /IOUNIT/ NIN,NOUT
      INTEGER NLOOPS,NPROPS,NVERTS,CUTMAX
      COMMON /SIZES/ NLOOPS,NPROPS,NVERTS,CUTMAX
      INTEGER C(2*SIZE,3),V(2*SIZE)
      INTEGER L,I,VV
C
      DO VV = 1,NVERTS
      DO I = 1,3
       C(VV,I) = CIN(VV,I)
      ENDDO
      ENDDO
      L = NVERTS
      NPERMS = 0
      OK = .TRUE.
C
C "CALL PERMUTATIONS(L,C)"
C
C-----
C "SUBROUTINE PERMUTATIONS(L,C)"
C Mock subroutine that generates each element of the permutation
C group S_(L-2), applies it to C, and calls CHECKOUT(C,CIN,N,OK).
C If OK = False is returned, the graph C was no good and we exit
C from CHECK immediately.
C
  1   CONTINUE
      IF (L.EQ.4) THEN
         CALL CHECKOUT(C,CIN,NPERMS,OK)
         IF (.NOT.OK) RETURN
         CALL EXCHANGE(3,4,C)
         CALL CHECKOUT(C,CIN,NPERMS,OK)
         IF (.NOT.OK) RETURN
         CALL EXCHANGE(3,4,C)
C "RETURN"
         GO TO 2
      ENDIF
C "DO V(L) = L,3,-1"
      V(L) = L
  3   CONTINUE
      CALL EXCHANGE(V(L),L,C)
      L = L - 1
C "CALL PERMUTATIONS(L,C)"
      GO TO 1
C Return from mock subroutine comes here:
  2   CONTINUE
      L = L + 1
      CALL EXCHANGE(V(L),L,C)
      V(L) = V(L) - 1
C "ENDDO"
      IF (V(L).GE.3) THEN
        GO TO 3
      ENDIF
C "RETURN" Executed from level L as long as L < NVERTS, else we are done.
      IF (L.LT.NVERTS) THEN
        GO TO 2
      ENDIF
C-----
C Come to here if graph CIN was OK, with NPERMS = number of permutations
C that gave a distinct numbering of the vertices.
C
      RETURN
      END
C
C23456789012345678901234567890123456789012345678901234567890123456789012
C
      SUBROUTINE CHECKOUT(C,CIN,NPERMS,OK)
C
      INTEGER SIZE
      PARAMETER (SIZE = 3)
      INTEGER C(2*SIZE,3),CIN(2*SIZE,3),NPERMS
      LOGICAL OK
C
C Test if graph C (with vertices permuted) is "less than," or
C "greater than," or equal to the original graph CIN using
C the standard ordering of graphs.  If C > CIN we leave unchanged the
C count NPERMS of how many vertex interchanges give the same graph and
C return OK = True.  If C = CIN, we add one to NPERMS, and we
C still return OK = True.  If C < CIN, then we should not have
C generated CIN, so we return OK = False.
C
      INTEGER NIN,NOUT
      COMMON /IOUNIT/ NIN,NOUT
      INTEGER NLOOPS,NPROPS,NVERTS,CUTMAX
      COMMON /SIZES/ NLOOPS,NPROPS,NVERTS,CUTMAX
      INTEGER V,I
C
      DO V = 1,NVERTS
      DO I = 1,3
       IF (C(V,I).LT.CIN(V,I)) THEN
         OK = .FALSE.
         RETURN
       ELSE IF (C(V,I).GT.CIN(V,I)) THEN
         OK = .TRUE.
         RETURN
       ENDIF
      ENDDO
      ENDDO
C
C Come to here if the new graph C is the same as CIN.
C
      OK = .TRUE.
      NPERMS = NPERMS + 1
      RETURN
      END
C
C23456789012345678901234567890123456789012345678901234567890123456789012
C
      SUBROUTINE EXCHANGE(V1,V2,C)
C
      INTEGER SIZE
      PARAMETER (SIZE = 3)
      INTEGER C(2*SIZE,3)
      INTEGER V1,V2
C
      INTEGER NIN,NOUT
      COMMON /IOUNIT/ NIN,NOUT
      INTEGER NLOOPS,NPROPS,NVERTS,CUTMAX
      COMMON /SIZES/ NLOOPS,NPROPS,NVERTS,CUTMAX
      INTEGER TEMP1,TEMP2,I,V
      LOGICAL CHANGE
C
      DO I = 1,3
        TEMP1 = C(V1,I)
        TEMP2 = C(V2,I)
        C(V1,I) = TEMP2
        C(V2,I) = TEMP1
      ENDDO
C
      DO V = 1,NVERTS
        CHANGE = .FALSE.
        DO I = 1,3
          IF (C(V,I).EQ.V1) THEN
            C(V,I) = V2
            CHANGE = .TRUE.
          ELSE IF (C(V,I).EQ.V2) THEN
            C(V,I) = V1
            CHANGE = .TRUE.
          ENDIF
        ENDDO
        IF (CHANGE) THEN
C
C Put vertices connected to vertex V in order
C--
          IF (C(V,2).LT.C(V,1)) THEN
            TEMP1 = C(V,1)
            TEMP2 = C(V,2)
            C(V,1) = TEMP2
            C(V,2) = TEMP1
          ENDIF
          IF (C(V,3).LT.C(V,1)) THEN
            TEMP1 = C(V,1)
            TEMP2 = C(V,3)
            C(V,1) = TEMP2
            C(V,3) = TEMP1
          ELSE IF (C(V,3).LT.C(V,2)) THEN
            TEMP1 = C(V,2)
            TEMP2 = C(V,3)
            C(V,2) = TEMP2
            C(V,3) = TEMP1
          ENDIF
C--
        ENDIF
      ENDDO
C
      RETURN
      END
C
C23456789012345678901234567890123456789012345678901234567890123456789012
C            End of subroutines associated with NEWGRAPH               C
C23456789012345678901234567890123456789012345678901234567890123456789012
C
      SUBROUTINE FINDA(VRTX,Q,NQ,A,QOK)
C
      INTEGER SIZE
      PARAMETER (SIZE = 3)
C In:
      INTEGER VRTX(0:3*SIZE-1,2),Q(0:SIZE),NQ
C Out:
      INTEGER A(0:3*SIZE-1,0:SIZE)
      LOGICAL QOK
C
C Finds matrix A relating propagator momenta to loop momenta.
C
C VRTX(P,N)     specifies the graph considered
C Q(L)           specifies the propagators to be considered independent
C NQ             specifies how many entries of Q should be considered
C   NQ = NLOOPS  all the entries in Q should be considered.
C                If Q(0),Q(1),...,Q(NLOOPS) are independent then
C                FINDA generates the matrix A and sets QOK = .TRUE.
C                Otherwise the generation of A fails and QOK = .FALSE.
C   NQ < NLOOPS  only first NQ entries in Q should be considered.
C                If Q(0),Q(1),...,Q(NQ) are independent then
C                FINDA sets QOK = .TRUE.
C                Otherwise  QOK = .FALSE.
C                In either case, a complete A is not generated.
C
C L              index of loop momenta, L = 0,1,...,NLOOPS.
C                L = 0 normally denontes the virtual photon momentum.
C P              index of propagator, P = 0,1,...,NPROPS.
C                P = 0 denotes the virtual photon momentum.
C V              index of vertices, V = 1,...,NVERTS
C A(P,L)         matrix relating propagator momenta to loop momenta.
C                K(P) = Sum_L A(P,L) K(Q(L)).
C VRTX(P,1) = V means that the vertex connected to the tail of
C                propagator P is V.
C VRTX(P,2) = V means that the vertex connected to the head of
C                propagator P is V.
C Q(L) = P       means that we consider the Lth loop momentum to
C                be that carried by propagator P.
C CONNECTED(V,J) = P means that the Jth propagator connected to
C                vertex V is P.
C FIXED(P) = True means that we have determined the momentum carried
C                by propagator P.
C FINISHED(V) = True means that we have determined the momenta carried
C                by all the propagators connected to vertex V.
C PROPSIGN(VRTX,P,V) is a function that returns +1 if the head of
C                propagator P is at V, -1 if the tail is at V.
C COUNT          is the number of propagators connected to the vertex
C                under consideration such that FIXED(P) = True. If
C                COUNT = 2, then we can fix another propagator momentum.
C 3 July 1994
C 19 December 1995
C
      INTEGER NIN,NOUT
      COMMON /IOUNIT/ NIN,NOUT
      INTEGER NLOOPS,NPROPS,NVERTS,CUTMAX
      COMMON /SIZES/ NLOOPS,NPROPS,NVERTS,CUTMAX
C
      INTEGER L,P,V,J,L1,L2
      INTEGER CONNECTED(2*SIZE,3)
      LOGICAL FIXED(0:3*SIZE-1),FINISHED(2*SIZE)
      LOGICAL CHANGE
      INTEGER PROPSIGN,SIGN
      INTEGER SUM(0:SIZE)
      INTEGER COUNT
      INTEGER PTOFIX
C
      IF((NQ.LT.1).OR.(NQ.GT.NLOOPS)) THEN
        WRITE(NOUT,*)'NQ out of range in FINDA'
      ENDIF
C
C First check to see that the same propagator hasn't been
C assigned to two loop variables.
C
      DO L1 = 0,NQ-1
       DO L2 = L1+1,NQ
        IF (Q(L1).EQ.Q(L2)) THEN
          QOK = .FALSE.
          RETURN
        ENDIF
       ENDDO
      ENDDO
C
C Initialization.
C
      QOK = .FALSE.
C
      DO V = 1,NVERTS
        J = 1
        DO P = 0,NPROPS
         IF( (VRTX(P,1).EQ.V).OR.(VRTX(P,2).EQ.V) ) THEN
          CONNECTED(V,J) = P
          J = J+1
         ENDIF
        ENDDO
      ENDDO
C
      DO P = 0,NPROPS
        DO L = 0,NLOOPS
          A(P,L) = 0
        ENDDO
      ENDDO
      DO L = 0,NQ
        A(Q(L),L) = 1
      ENDDO
C
      DO P = 0,NPROPS
        FIXED(P) = .FALSE.
      ENDDO
      DO L = 0,NQ
        FIXED(Q(L)) = .TRUE.
      ENDDO
C
      DO V = 1,NVERTS
        FINISHED(V) = .FALSE.
      ENDDO
C
      CHANGE = .TRUE.
C
C Start.
C
      DO WHILE (CHANGE)
      CHANGE = .FALSE.
C
      DO V = 1,NVERTS
      IF (.NOT.FINISHED(V)) THEN
        COUNT = 0
        DO J = 1,3
         P = CONNECTED(V,J)
         IF ( FIXED(P) ) THEN
           COUNT = COUNT + 1
         ENDIF
        ENDDO
C
C There are 3 already fixed propagators conencted to this vertex, so
C we must check to see if the momenta coming into the vertex sum to
C zero.
C
        IF (COUNT.EQ.3) THEN
          DO L = 0,NQ
            SUM(L) = 0
          ENDDO
          DO J = 1,3
            P = CONNECTED(V,J)
            SIGN = PROPSIGN(VRTX,P,V)
            DO L = 0,NQ
              SUM(L) = SUM(L) + SIGN * A(P,L)
            ENDDO
          ENDDO
          DO L = 0,NQ
C
C Dependent propagators given to FINDA.
C
            IF (SUM(L).NE.0) THEN
              QOK = .FALSE.
              RETURN
C
            ENDIF
          ENDDO
          FINISHED(V) = .TRUE.
          CHANGE = .TRUE.
C
C There are two already fixed propagators connected to this vertex,
C so we should determine the momentum carried by the remaining,
C unfixed, propagator.
C
        ELSE IF (COUNT.EQ.2) THEN
          DO L = 0,NQ
            SUM(L) = 0
          ENDDO
          DO J = 1,3
            P = CONNECTED(V,J)
            IF ( FIXED(P) ) THEN
              SIGN = PROPSIGN(VRTX,P,V)
              DO L = 0,NQ
                SUM(L) = SUM(L) + SIGN * A(P,L)
              ENDDO
            ELSE
              PTOFIX = P
            ENDIF
          ENDDO
          SIGN = PROPSIGN(VRTX,PTOFIX,V)
          DO L = 0,NQ
            A(PTOFIX,L) =  - SIGN * SUM(L)
          ENDDO
          FIXED(PTOFIX) = .TRUE.
          FINISHED(V) = .TRUE.
          CHANGE = .TRUE.
        ENDIF
C
C Close loop DO V = 1,NVERTS ; IF (.NOT.FINISHED(V)) THEN.
C
      ENDIF
      ENDDO
C
C Close loop DO WHILE (CHANGE)
C
      ENDDO
C
C At this point, we have not found a contradiction with momentum
C conservation, so the Q's must have been OK:
C
      QOK = .TRUE.
C
C If we had been given a complete set of Q's, then we should have
C fixed each propagator at each vertex. Check just to make sure.
C
      IF (NQ.EQ.NLOOPS) THEN
      DO V = 1,NVERTS
        IF (.NOT.FINISHED(V) ) THEN
          WRITE(NOUT,*)'SNAFU in FINDA'
          write(nout,*)'v = ',v,'  nq =',nq
          write(nout,*)'q =',q(0),q(1),q(2),q(3)
          STOP
        ENDIF
      ENDDO
      ENDIF
C
      RETURN
C
      END
C
C23456789012345678901234567890123456789012345678901234567890123456789012
C
      INTEGER FUNCTION PROPSIGN(VRTX,P,V)
C
      INTEGER SIZE
      PARAMETER (SIZE = 3)
C In:
      INTEGER VRTX(0:3*SIZE-1,2)
      INTEGER P,V
C
      INTEGER NIN,NOUT
      COMMON /IOUNIT/ NIN,NOUT
      INTEGER NLOOPS,NPROPS,NVERTS,CUTMAX
      COMMON /SIZES/ NLOOPS,NPROPS,NVERTS,CUTMAX
C
      IF ( VRTX(P,1).EQ.V ) THEN
          PROPSIGN = -1
          RETURN
      ELSE IF ( VRTX(P,2).EQ.V ) THEN
          PROPSIGN = 1
          RETURN
      ELSE
          WRITE(NOUT,*)'PROPSIGN called for P not connected to V.'
          STOP
      ENDIF
      END
C
C23456789012345678901234567890123456789012345678901234567890123456789012
C23456789012345678901234567890123456789012345678901234567890123456789012
C
      SUBROUTINE FINDTYPES(VRTX,PROP,NMAPS,QS,QSIGNS,MAPTYPES)
C
      INTEGER SIZE,MAXMAPS
      PARAMETER (SIZE = 3)
      PARAMETER (MAXMAPS = 64)
C In:
      INTEGER VRTX(0:3*SIZE-1,2),PROP(2*SIZE,3)
C Out:
      INTEGER NMAPS,QS(MAXMAPS,0:SIZE),QSIGNS(MAXMAPS,0:SIZE)
      CHARACTER*6 MAPTYPES(MAXMAPS)
C
C Given a graph specified by VRTX and PROP, this
C subroutine finds the characteristics of each map,
C labelled by an index MAPNUMBER. For a given map, it
C finds labels Q of 'special' propagators and
C the corresponding signs QSIGN and the MAPTYPE.
C The subroutine does finds the total number
C of maps, NMAPS, and fills the corresponding
C arrays QS, QSIGNS, and MAPTYPES, each of which
C carries a MAPNUMBER index.
C
C The possibilities for maptypes are as follows:
C 
C 1) T2TO2T used for k1 + k2 -> p1 + p2 with a virtual parton
C with momentum q exchanged, q = k1 - p1. Then P(Q(1)) = q,
C P(Q(2)) = p1, P(Q(3)) = p2. We will generate points with
C CHOOSET2TO2T(p1,p2,q,ok).
C 
C 2) T2TO2S used for k1 + k2 -> p1 + p2 with a *no* virtual
C parton with momentum q = k1 - p1 exchanged. Then P(Q(1)) = k1,
C P(Q(2)) = p1, P(Q(3)) = p2. We will generate points with
C CHOOSE2TO2S(p1,p2,k1,ok).
C 
C 3) T2TO3  used for k1 + k2 -> p1 + p2 + p3 with k2 = - k1.
C Then P(Q(1)) = k1, P(Q(2)) = p1, P(Q(3)) = p2. We will generate
C points with CHOOSET2TO3(p1,p2,k1,ok).
C 
C 4) T2TO1  used for k1 + k2 -> p1 on shell. We will choose points
C with CHOOSEST2TO1(p1,p2,k1,ok).
C
C 20 December 2000
C 20 March  2001
C
      INTEGER NIN,NOUT
      COMMON /IOUNIT/ NIN,NOUT
      INTEGER NLOOPS,NPROPS,NVERTS,CUTMAX
      COMMON /SIZES/ NLOOPS,NPROPS,NVERTS,CUTMAX
C
      INTEGER MAPNUMBER
      LOGICAL MORENEEDED
C Newcut variables
      LOGICAL NEWCUTINIT
      INTEGER CUTINDEX(SIZE+1),CUTSIGN(SIZE+1),NCUT
      INTEGER ISIGN(3*SIZE-1)
      LOGICAL LEFTLOOP,RIGHTLOOP
      INTEGER LOOPINDEX(SIZE+1),LOOPSIGN(SIZE+1),NINLOOP
      LOGICAL CUTFOUND
C
      INTEGER L,P,KJ,K1,K2,KDIRECT,PTEST,PLEAVING,PP1,PP2
      INTEGER L1,SIGNL1,LTESTA,LTESTB
      INTEGER I,J,JFOUND1,JFOUND2
      INTEGER V1,V2,V3,VOTHER,VV1,VV2
      INTEGER SIGN0,SIGN1,SIGN2
      INTEGER TIMESFOUND1,TIMESFOUND2,TIMESFOUND
      LOGICAL NOTINLOOP
C
C----------------------------------
C
      MAPNUMBER = 0
      MORENEEDED = .TRUE.
      NEWCUTINIT = .TRUE.
      DO WHILE (MORENEEDED)
      CALL NEWCUT(VRTX,NEWCUTINIT,NCUT,ISIGN,
     >            CUTINDEX,CUTSIGN,LEFTLOOP,RIGHTLOOP,
     >            NINLOOP,LOOPINDEX,LOOPSIGN,CUTFOUND)
      IF (CUTFOUND) THEN
C
C We want to do something only if there is a virtual loop:
C
      IF (NCUT.EQ.(CUTMAX-1)) THEN
C---
C Case of 4 propagators in the loop
C---
      IF (NINLOOP.EQ.4) THEN
C
C For a 4 propagator loop there are two ellipse maps (T2T02T) and 
C one circle map (T2TO3). We do the two ellipse maps first.
C
      DO L = 2,3
          MAPNUMBER = MAPNUMBER + 1
          P = LOOPINDEX(L)
          V1 = VRTX(P,1)
          V2 = VRTX(P,2)
C
C We find the cut propagators K1 and K2 connected to V1 and V2 along
C with the sign = +1 if the cut propagator Kj is leaving vertex Vj
C and sign = -1 if the cut propagator Kj is entering vertex Vj. Just
C as a check, we define FOUNDJ to see if we find K1 and K2 exactly
C once.
C
          TIMESFOUND1 = 0
          TIMESFOUND2 = 0
          DO J = 1,3
            KJ = CUTINDEX(J)
            IF (VRTX(KJ,1).EQ.V1) THEN
              K1 = KJ
              SIGN1 = +1
              TIMESFOUND1 = TIMESFOUND1+1
            ELSE IF (VRTX(KJ,2).EQ.V1) THEN
              K1 = KJ
              SIGN1 = -1
              TIMESFOUND1 = TIMESFOUND1+1
            ELSE IF (VRTX(KJ,1).EQ.V2) THEN
              K2 = KJ
              SIGN2 = +1
              TIMESFOUND2 = TIMESFOUND2+1
            ELSE IF (VRTX(KJ,2).EQ.V2) THEN
              K2 = KJ
              SIGN2 = -1
              TIMESFOUND2 = TIMESFOUND2+1
            ENDIF
          ENDDO
          IF ((TIMESFOUND1.NE.1).OR.(TIMESFOUND2.NE.1)) THEN
            WRITE(NOUT,*) 'Failure in FINDTYPES'
            STOP
          ENDIF
C
C Now we record this information.
C
          QS(MAPNUMBER,0) = 0
          QSIGNS(MAPNUMBER,0) = +1
          QS(MAPNUMBER,1) = P
          QSIGNS(MAPNUMBER,1) = +1
          QS(MAPNUMBER,2) = K1
          QSIGNS(MAPNUMBER,2) = SIGN1
          QS(MAPNUMBER,3) = K2
          QSIGNS(MAPNUMBER,3) = SIGN2
          MAPTYPES(MAPNUMBER) = 'T2TO2T'
C
C End  DO L = 2,3 for the choice of two ellipse maps.
C
      ENDDO
C
C Now we do the circle map.
C Our definition for the circle map T2TO3E is that Q(1) is
C LOOPINDEX(1) the first propagator in the loop starting from the
C current vertex.  Then Q(2) is the label of the propagator that
C enters the final state and connects to the vertex at the head
C of propagator Q(1).  Then Q(3) is the label of the propagator
C that enters the final state and connects to the propagator with
C label LOOPINDEX(4), the last propagator in the loop.  The sign
C QSIGN(1) = +1 since this propagator always points *from* the
C current vertex. For QSIGN(2) and QSIGN(3) a plus sign indicates
C that the propagator points toward the final state, a minus sign
C indicates the opposite.
C
      IF(LOOPSIGN(1).NE.1) THEN
        WRITE(NOUT,*)'LOOPSIGN(1) not 1 in FINDTYPES'
        STOP
      ELSE IF(LOOPSIGN(4).NE.-1) THEN
        WRITE(NOUT,*)'LOOPSIGN(4) not -1 in FINDTYPES'
        STOP
      ENDIF
C
      V1 = VRTX(LOOPINDEX(1),2)
      V2 = VRTX(LOOPINDEX(4),2)
      TIMESFOUND1 = 0
      TIMESFOUND2 = 0
      DO J = 1,3
        KJ = CUTINDEX(J)
        IF (VRTX(KJ,1).EQ.V1) THEN
          K1 = KJ
          SIGN1 = +1
          TIMESFOUND1 = TIMESFOUND1+1
        ELSE IF (VRTX(KJ,2).EQ.V1) THEN
          K1 = KJ
          SIGN1 = -1
          TIMESFOUND1 = TIMESFOUND1+1
        ELSE IF (VRTX(KJ,1).EQ.V2) THEN
          K2 = KJ
          SIGN2 = +1
          TIMESFOUND2 = TIMESFOUND2+1
        ELSE IF (VRTX(KJ,2).EQ.V2) THEN
          K2 = KJ
          SIGN2 = -1
          TIMESFOUND2 = TIMESFOUND2+1
        ENDIF
      ENDDO
      IF ((TIMESFOUND1.NE.1).OR.(TIMESFOUND2.NE.1)) THEN
        WRITE(NOUT,*) 'Oops, failure in FINDTYPES',
     >                 TIMESFOUND1,TIMESFOUND2
        STOP
      ENDIF
C
      MAPNUMBER = MAPNUMBER + 1
      QS(MAPNUMBER,0) = 0
      QSIGNS(MAPNUMBER,0) = +1
      QS(MAPNUMBER,1) = LOOPINDEX(1)
      QSIGNS(MAPNUMBER,1) = +1
      QS(MAPNUMBER,2) = K1
      QSIGNS(MAPNUMBER,2) = SIGN1
      QS(MAPNUMBER,3) = K2
      QSIGNS(MAPNUMBER,3) = SIGN2
      MAPTYPES(MAPNUMBER) = 'T2TO3E'
C
C---
C Case of 3 propagators in the loop
C---
      ELSE IF (NINLOOP.EQ.3) THEN
C
C We are not sure which of two possibilities we have, but we proceed
C as if we had the case of a virtual loop that connects to two 
C propagators that go into the final state.
C
      MAPNUMBER = MAPNUMBER + 1
      P = LOOPINDEX(2)
      V1 = VRTX(P,1)
      V2 = VRTX(P,2)
C
C We find the cut propagators K1 and K2 connected to V1 and V2 along
C with the sign = +1 if the cut propagator Kj is leaving vertex Vj
C and sign = -1 if the cut propagator Kj is entering vertex Vj. We
C check using FOUNDJ to see if we find K1 and K2 exactly once.
C
      TIMESFOUND1 = 0
      TIMESFOUND2 = 0
      DO J = 1,3
        KJ = CUTINDEX(J)
        IF (VRTX(KJ,1).EQ.V1) THEN
          K1 = KJ
          SIGN1 = +1
          TIMESFOUND1 = TIMESFOUND1+1
        ELSE IF (VRTX(KJ,2).EQ.V1) THEN
          K1 = KJ
          SIGN1 = -1
          TIMESFOUND1 = TIMESFOUND1+1
        ELSE IF (VRTX(KJ,1).EQ.V2) THEN
          K2 = KJ
          SIGN2 = +1
          TIMESFOUND2 = TIMESFOUND2+1
        ELSE IF (VRTX(KJ,2).EQ.V2) THEN
          K2 = KJ
          SIGN2 = -1
          TIMESFOUND2 = TIMESFOUND2+1
        ENDIF
      ENDDO
C
C Now we figure out what to do based on what we found.
C
      IF ((TIMESFOUND1.GT.1).OR.(TIMESFOUND2.GT.1)) THEN
         WRITE(NOUT,*) 'Failure in FINDTYPES'
         STOP
      ELSE IF ((TIMESFOUND1.LT.1).AND.(TIMESFOUND2.LT.1)) THEN
         WRITE(NOUT,*) 'Failure in FINDTYPES'
         STOP
C
      ELSE IF ((TIMESFOUND1.EQ.1).AND.(TIMESFOUND2.EQ.1)) THEN
C
C This is the case we were looking for. Now we record the information.
C
      QS(MAPNUMBER,0) = 0
      QSIGNS(MAPNUMBER,0) = +1
      QS(MAPNUMBER,1) = P
      QSIGNS(MAPNUMBER,1) = +1
      QS(MAPNUMBER,2) = K1
      QSIGNS(MAPNUMBER,2) = SIGN1
      QS(MAPNUMBER,3) = K2
      QSIGNS(MAPNUMBER,3) = SIGN2
      MAPTYPES(MAPNUMBER) = 'T2TO2T'
C
      ELSE
C
C Either Found1 = 1 and Found2 = 0 or Found2 = 1 and Found1 = 0.
C In these cases our loop does *not* connect to two propagators 
C that go to the final state. The label of the propagator
C that enters the final state will be called Kdirect and the 
C vertex that does not connect to this propagator will be called
C Vother. We take sign0 = +1 if our loop propagator points from
C Kdirect to the s-channel propagator that splits into two
C propagators that go to the final state. Otherwise sign0 = -1.
C
      IF ((TIMESFOUND1.EQ.1).AND.(TIMESFOUND2.EQ.0)) THEN
        KDIRECT = K1
        VOTHER = V2
        SIGN0 = +1
      ELSE IF ((TIMESFOUND1.EQ.0).AND.(TIMESFOUND2.EQ.1)) THEN
        KDIRECT = K2
        VOTHER = V1
        SIGN0 = -1
      ENDIF
C
C Now we deal with this case.
C
      IF (CUTINDEX(1).EQ.KDIRECT) THEN
        K1 = CUTINDEX(2)
        K2 = CUTINDEX(3)
      ELSE IF (CUTINDEX(2).EQ.KDIRECT) THEN
        K1 = CUTINDEX(3)
        K2 = CUTINDEX(1)
      ELSE IF (CUTINDEX(3).EQ.KDIRECT) THEN
        K1 = CUTINDEX(1)
        K2 = CUTINDEX(2)
      ELSE
        WRITE(NOUT,*)'We are in real trouble here.'
        STOP
      ENDIF
C
C We have K1 and K2, but we need the corresponding signs.
C Find the index Pleaving of the propagator leaving the loop toward
C the final state.
C
      TIMESFOUND = 0
      DO J = 1,3
        PTEST = PROP(VOTHER,J)
        NOTINLOOP = .TRUE.
        DO I = 1,3
           IF (PTEST.EQ.LOOPINDEX(I)) THEN
              NOTINLOOP = .FALSE.
           ENDIF
        ENDDO
        IF (NOTINLOOP) THEN
           PLEAVING = PTEST
           TIMESFOUND = TIMESFOUND + 1
        ENDIF
      ENDDO
      IF (TIMESFOUND.NE.1) THEN
         WRITE(NOUT,*)'Pleaving not found or found twice.'
         STOP
      ENDIf
C
C Let V3 be the vertex not in the loop at the end of propagator
C Pleaving. Two propagators in the final state must connect to this
C vertex.
C
      V3 = VRTX(PLEAVING,1)
      IF (V3.EQ.VOTHER) THEN
        V3 = VRTX(PLEAVING,2)
      ENDIF
C
C We use V3 to get the proper signs.
C
      IF (VRTX(K1,1).EQ.V3) THEN
          SIGN1 = +1
      ELSE IF (VRTX(K1,2).EQ.V3) THEN
          SIGN1 = -1
      ELSE
        WRITE(NOUT,*)'Yikes, this is bad.'
        STOP
      ENDIF
      IF (VRTX(K2,1).EQ.V3) THEN
          SIGN2 = +1
      ELSE IF (VRTX(K2,2).EQ.V3) THEN
          SIGN2 = -1
      ELSE
        WRITE(NOUT,*)'Yikes, this is also bad.'
        STOP
      ENDIF
C
C Now we record the information.
C Recall that P = LOOPINDEX(2) and that SIGN0 = +1 if propagator P
C points toward propagator Pleaving.
C
      QS(MAPNUMBER,0) = 0
      QSIGNS(MAPNUMBER,0) = +1
      QS(MAPNUMBER,1) = P
      QSIGNS(MAPNUMBER,1) = SIGN0
      QS(MAPNUMBER,2) = K1
      QSIGNS(MAPNUMBER,2) = SIGN1
      QS(MAPNUMBER,3) = K2
      QSIGNS(MAPNUMBER,3) = SIGN2
      MAPTYPES(MAPNUMBER) = 'T2TO2S'
C
C But we are not done, because in this case we need a circle map too.
C Our definition for the circle map T2TO3D is that Q(1) is
C LOOPINDEX(1) or LOOPINDEX(3), one of the two propagators that
C connects to a propagator that connects to the current vertex.
C We take the one that connects to vertex Vother that connects to
C a propagator Pleaving that connects vertex V3 that, finally,
C connects to to two propagators that enter the final state. Then
C Q(3) and Q(3) are these two propagators that enter the final
C state from vertex V3.  For QSIGN(2) and QSIGN(3) a plus sign
C indicates that the propagator points toward the final state, a
C minus sign indicates the opposite.  The sign QSIGN(1) is + 1 if
C this propagator points toward the final state, -1 in the
C opposite circumstance.
C
      IF (LOOPSIGN(1).NE.1) THEN
        WRITE(NOUT,*)'LOOPSIGN not 1 in FINDTYPES'
        STOP
      ENDIF
C
C The loop momentum with label L1 is the one that
C attaches to VOTHER (the vertex that connects to a propagator
C that splits before going to the final state.) We take
C SIGNL1 = +1 if this propagator points towards VOTHER.
C
      LTESTA = LOOPINDEX(1)
      LTESTB = LOOPINDEX(3)
      IF (VRTX(LTESTA,2).EQ.VOTHER) THEN
          L1 = LTESTA
          SIGNL1 = +1
      ELSE IF (VRTX(LTESTA,1).EQ.VOTHER) THEN
          L1 = LTESTA
          SIGNL1 = -1
      ELSE IF (VRTX(LTESTB,2).EQ.VOTHER) THEN
          L1 = LTESTB
          SIGNL1 = +1
      ELSE IF (VRTX(LTESTB,1).EQ.VOTHER) THEN
          L1 = LTESTB
          SIGNL1 = -1
      ELSE
        WRITE(NOUT,*)'Cannot seem to find L1'
        STOP
      ENDIF
C
      MAPNUMBER = MAPNUMBER + 1
      QS(MAPNUMBER,0) = 0
      QSIGNS(MAPNUMBER,0) = +1
      QS(MAPNUMBER,1) = L1
      QSIGNS(MAPNUMBER,1) = SIGNL1
      QS(MAPNUMBER,2) = K1
      QSIGNS(MAPNUMBER,2) = SIGN1
      QS(MAPNUMBER,3) = K2
      QSIGNS(MAPNUMBER,3) = SIGN2
      MAPTYPES(MAPNUMBER) = 'T2TO3D'
C
C Close the IF structure for the case of three particles in the loop,
C IF (FOUND1.GT.1).OR.(FOUND2.GT.1) THEN ...
C
      ENDIF
C---
C Case of 2 propagators in the loop
C---
      ELSE IF (NINLOOP.EQ.2) THEN
C
C We are not sure which of two possibilities we have, but we proceed
C as if we had the case of a virtual loop that connects to two 
C propagators that go into the final state.
C
      MAPNUMBER = MAPNUMBER + 1
      P = LOOPINDEX(1)
      V1 = VRTX(P,1)
      V2 = VRTX(P,2)
C
C We find the cut propagators K1 or K2 connected to V1 or V2 along
C with the sign = +1 if the cut propagator Kj is leaving vertex Vj
C and sign = -1 if the cut propagator Kj is entering vertex Vj. We
C check using FOUNDJ to see if we find K1 or K2 exactly once.
C
      TIMESFOUND1 = 0
      TIMESFOUND2 = 0
      DO J = 1,3
        KJ = CUTINDEX(J)
        IF (VRTX(KJ,1).EQ.V1) THEN
          K1 = KJ
          SIGN1 = +1
          JFOUND1 = J
          TIMESFOUND1 = TIMESFOUND1+1
        ELSE IF (VRTX(KJ,2).EQ.V1) THEN
          K1 = KJ
          SIGN1 = -1
          JFOUND1 = J
          TIMESFOUND1 = TIMESFOUND1+1
        ELSE IF (VRTX(KJ,1).EQ.V2) THEN
          K2 = KJ
          SIGN2 = +1
          JFOUND2 = J
          TIMESFOUND2 = TIMESFOUND2+1
        ELSE IF (VRTX(KJ,2).EQ.V2) THEN
          K2 = KJ
          SIGN2 = -1
          JFOUND2 = J
          TIMESFOUND2 = TIMESFOUND2+1
        ENDIF
      ENDDO
C
C Now we figure out what to do based on what we found.
C
      IF ((TIMESFOUND1.GT.1).OR.(TIMESFOUND2.GT.1)) THEN
         WRITE(NOUT,*) 'Failure in FINDTYPES'
         STOP
C
      ELSE IF ((TIMESFOUND1.EQ.1).AND.(TIMESFOUND2.EQ.0)) THEN
C
C This is one of the cases we were looking for. Now we record the 
C information. The propagator Q(3) is one of the propagators
C in the final state other than that connected to our loop. The
C corresponding sign is +1 if this propagator crosses the final
C state cut in the same direction as the propagator connected to
C our loop.
C
      QS(MAPNUMBER,0) = 0
      QSIGNS(MAPNUMBER,0) = +1
      QS(MAPNUMBER,1) = P
      QSIGNS(MAPNUMBER,1) = -1
      QS(MAPNUMBER,2) = K1
      QSIGNS(MAPNUMBER,2) = SIGN1
      IF (CUTINDEX(1).NE.K1) THEN
        QS(MAPNUMBER,3) = CUTINDEX(1)
        QSIGNS(MAPNUMBER,3) = CUTSIGN(1)*CUTSIGN(JFOUND1)*SIGN1
      ELSE
        QS(MAPNUMBER,3) = CUTINDEX(2)
        QSIGNS(MAPNUMBER,3) = CUTSIGN(2)*CUTSIGN(JFOUND1)*SIGN1
      ENDIF
      MAPTYPES(MAPNUMBER) = 'T2TO1 '
C
      ELSE IF ((TIMESFOUND1.EQ.0).AND.(TIMESFOUND2.EQ.1)) THEN
C
C This is one of the cases we were looking for. Now we record the 
C information. The propagator Q(3) is one of the propagators
C in the final state other than that connected to our loop. The
C corresponding sign is +1 if this propagator crosses the final
C state cut in the same direction as the propagator connected to
C our loop.
C
      QS(MAPNUMBER,0) = 0
      QSIGNS(MAPNUMBER,0) = +1
      QS(MAPNUMBER,1) = P
      QSIGNS(MAPNUMBER,1) = +1
      QS(MAPNUMBER,2) = K2
      QSIGNS(MAPNUMBER,2) = SIGN2
      IF (CUTINDEX(1).NE.K2) THEN
        QS(MAPNUMBER,3) = CUTINDEX(1)
        QSIGNS(MAPNUMBER,3) = CUTSIGN(1)*CUTSIGN(JFOUND2)*SIGN2
      ELSE
        QS(MAPNUMBER,3) = CUTINDEX(2)
        QSIGNS(MAPNUMBER,3) = CUTSIGN(2)*CUTSIGN(JFOUND2)*SIGN2
      ENDIF
      MAPTYPES(MAPNUMBER) = 'T2TO1 '
C
      ELSE IF ((TIMESFOUND1.EQ.0).AND.(TIMESFOUND2.EQ.0)) THEN
C
C Here TimesFound1 = 0 and TimesFound2 = 0, so our loop does *not* 
C connect to a propagator that goes to the final state.
C Find the indices PP1 and PP2 of the propagators connected to
C our loop.
C
      TIMESFOUND = 0
      DO J = 1,3
        PTEST = PROP(V1,J)
        NOTINLOOP = .TRUE.
        DO I = 1,2
           IF (PTEST.EQ.LOOPINDEX(I)) THEN
              NOTINLOOP = .FALSE.
           ENDIF
        ENDDO
        IF (NOTINLOOP) THEN
           PP1 = PTEST
           TIMESFOUND = TIMESFOUND + 1
        ENDIF
      ENDDO
      IF (TIMESFOUND.NE.1) THEN
         WRITE(NOUT,*)'PP1 not found or found twice.'
         STOP
      ENDIf
      TIMESFOUND = 0
C
      DO J = 1,3
        PTEST = PROP(V2,J)
        NOTINLOOP = .TRUE.
        DO I = 1,2
           IF (PTEST.EQ.LOOPINDEX(I)) THEN
              NOTINLOOP = .FALSE.
           ENDIF
        ENDDO
        IF (NOTINLOOP) THEN
           PP2 = PTEST
           TIMESFOUND = TIMESFOUND + 1
        ENDIF
      ENDDO
      IF (TIMESFOUND.NE.1) THEN
         WRITE(NOUT,*)'PP2 not found or found twice.'
         STOP
      ENDIf
C
C Let VV1 and VV2 be the vertices not in the loop at the end of
C propagators PP1 and PP2 respectively. Two propagators in the final 
C state must connect to one of these vertices.
C
      VV1 = VRTX(PP1,1)
      IF (VV1.EQ.V1) THEN
        VV1 = VRTX(PP1,2)
      ENDIF
      VV2 = VRTX(PP2,1)
      IF (VV2.EQ.V2) THEN
        VV2 = VRTX(PP2,2)
      ENDIF
C
C We have VV1 and VV2. A slight hitch is that one of them might be
C the vertex 1 or 2 that connect to the photon. In this case, 
C in the next step we do *not* want to find the final state
C propagator that connects to this vertex. A cure is to set the
C vertex number to something impossible.
C
      IF ((VV1.EQ.1).OR.(VV1.EQ.2)) THEN
        VV1 = -17
      ENDIF
      IF ((VV2.EQ.1).OR.(VV2.EQ.2)) THEN
        VV2 = -17
      ENDIF
C
C Now we find two final state propagators connected to VV1 or
C else two final state propagators connected to VV2.
C
      TIMESFOUND = 0
      DO J = 1,3
        KJ = CUTINDEX(J)
        IF (VRTX(KJ,1).EQ.VV1) THEN
          IF(TIMESFOUND.EQ.0) THEN
            K1 = KJ
            SIGN1 = +1
          ELSE
            K2 = KJ
            SIGN2 = +1
          ENDIF
          SIGN0 = -1
          TIMESFOUND = TIMESFOUND+1
        ELSE IF (VRTX(KJ,1).EQ.VV2) THEN
          IF(TIMESFOUND.EQ.0) THEN
            K1 = KJ
            SIGN1 = +1
          ELSE
            K2 = KJ
            SIGN2 = +1
          ENDIF
          SIGN0 = +1
          TIMESFOUND = TIMESFOUND+1
        ELSE IF (VRTX(KJ,2).EQ.VV1) THEN
          IF(TIMESFOUND.EQ.0) THEN
            K1 = KJ
            SIGN1 = -1
          ELSE
            K2 = KJ
            SIGN2 = -1
          ENDIF
          SIGN0 = -1
          TIMESFOUND = TIMESFOUND+1
        ELSE IF (VRTX(KJ,2).EQ.VV2) THEN
          IF(TIMESFOUND.EQ.0) THEN
            K1 = KJ
            SIGN1 = -1
          ELSE
            K2 = KJ
            SIGN2 = -1
          ENDIF
          SIGN0 = +1
          TIMESFOUND = TIMESFOUND+1
        ENDIF
      ENDDO
      IF (TIMESFOUND.NE.2) THEN
         WRITE(NOUT,*)'Where are those tricky propagators?',TIMESFOUND
         STOP
      ENDIf
C
C Now we record the information.
C Recall that P = LOOPINDEX(1) and that SIGN0 = +1 if propagator P
C points toward propagators in the final state.
C
      QS(MAPNUMBER,0) = 0
      QSIGNS(MAPNUMBER,0) = +1
      QS(MAPNUMBER,1) = P
      QSIGNS(MAPNUMBER,1) = SIGN0
      QS(MAPNUMBER,2) = K1
      QSIGNS(MAPNUMBER,2) = SIGN1
      QS(MAPNUMBER,3) = K2
      QSIGNS(MAPNUMBER,3) = SIGN2
      MAPTYPES(MAPNUMBER) = 'T2TO2S'
C
C Close IF ((TIMESFOUND1.GT.1).OR.(TIMESFOUND2.GT.1)) THEN ...
C
      ENDIF
C
C Case of less than 2 propagators in the loop
C
      ELSE
        WRITE(NOUT,*) 'Looped the loop in FINDDQS'
        STOP
C
C End IF (NINLOOP.EQ. n ) series
C
      ENDIF
C
C End IF (there is a virtual loop) THEN ... 
C
      ENDIF
C
C End IF (CUTFOUND) THEN ... If the cut was not found, then we are done.
C
      ELSE
          MORENEEDED = .FALSE.
      ENDIF
C
C End main loop: DO WHILE (MORENEEDED)
C
      ENDDO
C
      NMAPS = MAPNUMBER
      RETURN
      END
C
C23456789012345678901234567890123456789012345678901234567890123456789012
C23456789012345678901234567890123456789012345678901234567890123456789012
C
      SUBROUTINE NEWPOINT(A,QSIGN,MAPTYPE,K,ABSK,BADPOINT)
C
      INTEGER SIZE
      PARAMETER (SIZE = 3)
C In:
      INTEGER A(0:3*SIZE-1,0:SIZE),QSIGN(0:SIZE)
      CHARACTER*6 MAPTYPE
C Out:
      REAL*8 K(0:3*SIZE-1,0:3),ABSK(0:3*SIZE-1)
      LOGICAL BADPOINT
C
C Chooses a new Monte Carlo point in the space of loop 3-momenta.
C 4 March 1993
C 12 July 1993
C 17 July 1994
C  2 May 1996
C  5 February 1997
C  4 February 1999
C 10 March 1999
C  9 April 1999
C 20 August 1999
C 21 December 2000 <cut based method for choosing points>
C 20 March 2001
C
      INTEGER NIN,NOUT
      COMMON /IOUNIT/ NIN,NOUT
      INTEGER NLOOPS,NPROPS,NVERTS,CUTMAX
      COMMON /SIZES/ NLOOPS,NPROPS,NVERTS,CUTMAX
C
      REAL*8 P1(3),P2(3),P3(3),ELL1(3)
      INTEGER P,MU
      REAL*8 TEMP,KSQ
      LOGICAL OK
C
C------------
C
      BADPOINT = .FALSE.
C
C Our notation: 
C   special loop momentum, to become QSIGN(1)*ELL(1,mu), is ELL1(mu);
C   first final state parton, to become QSIGN(2)*ELL(2,mu), is P1(mu);
C   second final state parton, to become QSIGN(3)*ELL(3,mu), is P2(mu);
C   third final state parton, not reported, is P3(mu).
C We use {ELL1,P1,P2} directly to generate the K(P,mu).
C
C First, we need to generate a three parton final state.
C Abort if we get a not OK signal.
C
      CALL CHOOSE3(P1,P2,P3,OK)
      IF(.NOT.OK) THEN
        DO P = 1,NPROPS
        DO MU = 0,3
          K(P,MU) = 0.0D0
        ENDDO
        ENDDO
        BADPOINT = .TRUE.
        RETURN
      ENDIF
C
C Then we generate the loop momentum, ell1.
C Abort if we get a not OK signal.
C
      IF (MAPTYPE.EQ.'T2TO3D') THEN
        CALL CHOOSE2TO3D(P1,P2,ELL1,OK)
        IF(.NOT.OK) THEN
          DO P = 1,NPROPS
          DO MU = 0,3
            K(P,MU) = 0.0D0
          ENDDO
          ENDDO
          BADPOINT = .TRUE.
          RETURN
        ENDIF
      ELSE IF (MAPTYPE.EQ.'T2TO3E') THEN
        CALL CHOOSE2TO3E(P1,P2,ELL1,OK)
        IF(.NOT.OK) THEN
          DO P = 1,NPROPS
          DO MU = 0,3
            K(P,MU) = 0.0D0
          ENDDO
          ENDDO
          BADPOINT = .TRUE.
          RETURN
        ENDIF
      ELSE IF (MAPTYPE.EQ.'T2TO2T') THEN
        CALL CHOOSE2TO2T(P1,P2,ELL1,OK)
        IF(.NOT.OK) THEN
          DO P = 1,NPROPS
          DO MU = 0,3
            K(P,MU) = 0.0D0
          ENDDO
          ENDDO
          BADPOINT = .TRUE.
          RETURN
        ENDIF
      ELSE IF (MAPTYPE.EQ.'T2TO2S') THEN
        CALL CHOOSE2TO2S(P1,P2,ELL1,OK)
        IF(.NOT.OK) THEN
          DO P = 1,NPROPS
          DO MU = 0,3
            K(P,MU) = 0.0D0
          ENDDO
          ENDDO
          BADPOINT = .TRUE.
          RETURN
        ENDIF
      ELSE IF (MAPTYPE.EQ.'T2TO1 ') THEN
        CALL CHOOSE2TO1(P1,P2,ELL1,OK)
        IF(.NOT.OK) THEN
          DO P = 1,NPROPS
          DO MU = 0,3
            K(P,MU) = 0.0D0
          ENDDO
          ENDDO
          BADPOINT = .TRUE.
          RETURN
        ENDIF
      ELSE 
        WRITE(NOUT,*)'Bad MAPTYPE in NEWPOINT'
        STOP
      ENDIF
C
C Now we have ELL1(mu), P1(mu), and P2(mu) and we need to translate to 
C the propagator momenta K(P,MU).
C
      DO P = 1,NPROPS
        KSQ = 0.0D0
        DO MU = 1,3
          TEMP =  A(P,1)*QSIGN(1)*ELL1(MU) 
     >          + A(P,2)*QSIGN(2)*P1(MU) 
     >          + A(P,3)*QSIGN(3)*P2(MU)
          K(P,MU) = TEMP
          KSQ = KSQ + TEMP**2
      ENDDO
        ABSK(P) = SQRT(KSQ)
        K(P,0) = 0.0D0
      ENDDO
      DO MU = 0,3
        K(0,MU) = 0.0D0
      ENDDO
      ABSK(0) = 0.0D0
C
      RETURN
      END
C
C23456789012345678901234567890123456789012345678901234567890123456789012
C23456789012345678901234567890123456789012345678901234567890123456789012
C
      SUBROUTINE CHECKPOINT(K,ABSK,PROP,BADNESS)
C
      INTEGER SIZE
      PARAMETER (SIZE = 3)
C In:
      REAL*8 K(0:3*SIZE-1,0:3),ABSK(0:3*SIZE-1)
      INTEGER PROP(2*SIZE,3)
C Out:
      REAL*8 BADNESS
C
C Calculates the BADNESS of a point chosen by NEWPOINT. If there
C are very collinear particles meeting at a vertex or of there is a
C very soft particle, then the badness is big. Specifically, for
C each vertex V we order the momenta entering the vertex Kmin, Kmid
C Kmax in order of |K|. Then
C
C   Kmin (Kmin + Kmid - Kmax )/Kmax^2
C
C is the 1/badness^2 for that vertex. The badness for the point is the 
C largest of the badness values of all the vertices.
C 
C Revised 13 may 1998
C
      INTEGER NIN,NOUT
      COMMON /IOUNIT/ NIN,NOUT
      INTEGER NLOOPS,NPROPS,NVERTS,CUTMAX
      COMMON /SIZES/ NLOOPS,NPROPS,NVERTS,CUTMAX
C
      REAL*8 SMALLNESSV,SMALLNESS
      INTEGER V
      REAL*8 KMIN,KMID,KMAX,K1,K2,K3
C
      SMALLNESS = 1.0D0
      DO V = 3,NVERTS
       K1 = ABSK(PROP(V,1))
       K2 = ABSK(PROP(V,2))
       K3 = ABSK(PROP(V,3))
       IF (K1.LT.K2) THEN
         KMIN = K1
         KMAX = K2
       ELSE  
         KMIN = K2
         KMAX = K1
       ENDIF
       IF (K3.LT.KMIN) THEN
         KMID = KMIN
         KMIN = K3
       ELSE IF (K3.GT.KMAX) THEN
         KMID = KMAX
         KMAX = K3
       ELSE
         KMID = K3
       ENDIF
       SMALLNESSV = KMIN * (KMIN + KMID - KMAX) /KMAX**2
       IF( SMALLNESSV .LT. SMALLNESS ) THEN
         SMALLNESS = SMALLNESSV
       ENDIF
       ENDDO
       IF (SMALLNESS.LT.1.0D-30) THEN
         BADNESS = 1.0D15
       ELSE
         BADNESS = SQRT(1.0D0/SMALLNESS)
       ENDIF
       RETURN
       END
C
C23456789012345678901234567890123456789012345678901234567890123456789012
C
      SUBROUTINE AXES(EA,EB,EC)
C
C In:
      REAL*8 EA(3)
C Out:
      REAL*8 EB(3),EC(3)
C
C Given a unit vector E_a(mu), generates unit vectors E_b(mu) and
C E_c(mu) such that E_a*E_b = E_b*E_c = E_c*E_a = 0.
C
C The vector E_b will lie in the plane formed by the z-axis and
C E_a unless E_a itself is nearly aligned along the z-axis, in which
C case E_b will lie in the plane formed by the x-axis and E_a.
C
C 18 April 1996
C
      INTEGER NIN,NOUT
      COMMON /IOUNIT/ NIN,NOUT
C
      REAL*8 COSTHETASQ,SINTHETAINV
C
C For check
C
      INTEGER MU
      REAL*8 DOTAA,DOTBB,DOTCC,DOTAB,DOTAC,DOTBC
C
      COSTHETASQ = EA(3)**2
      IF(COSTHETASQ.LT.0.9D0) THEN
        SINTHETAINV = 1.0D0/SQRT(1.0D0 - COSTHETASQ)
        EC(1) = - EA(2)*SINTHETAINV
        EC(2) =   EA(1)*SINTHETAINV
        EC(3) =   0.0D0
      ELSE
        COSTHETASQ = EA(1)**2
        SINTHETAINV = 1.0D0/SQRT(1.0D0 - COSTHETASQ)
        EC(1) =   0.0D0
        EC(2) = - EA(3)*SINTHETAINV
        EC(3) =   EA(2)*SINTHETAINV
      ENDIF
      EB(1) = EC(2)*EA(3) - EC(3)*EA(2)
      EB(2) = EC(3)*EA(1) - EC(1)*EA(3)
      EB(3) = EC(1)*EA(2) - EC(2)*EA(1) 
C     
C Check:
C
      DOTAA = 0.0D0
      DOTBB = 0.0D0
      DOTCC = 0.0D0
      DOTAB = 0.0D0
      DOTAC = 0.0D0
      DOTBC = 0.0D0
      DO MU = 1,3
        DOTAA = DOTAA + EA(MU)*EA(MU)
        DOTBB = DOTBB + EB(MU)*EB(MU)
        DOTCC = DOTCC + EC(MU)*EC(MU)
        DOTAB = DOTAB + EA(MU)*EB(MU)
        DOTAC = DOTAC + EA(MU)*EC(MU)
        DOTBC = DOTBC + EB(MU)*EC(MU)
      ENDDO
      IF (ABS(DOTAA - 1.0D0).GT.1.0D20) THEN
        WRITE(NOUT,*)'DOTAA messed up in AXES'
        STOP
      ELSE IF (ABS(DOTBB - 1.0D0).GT.1.0D20) THEN
        WRITE(NOUT,*)'DOTBB messed up in AXES'
        STOP
      ELSE IF (ABS(DOTCC - 1.0D0).GT.1.0D20) THEN
        WRITE(NOUT,*)'DOTCC messed up in AXES'
        STOP
      ELSE IF (ABS(DOTAB).GT.1.0D20) THEN
        WRITE(NOUT,*)'DOTAB messed up in AXES'
        STOP
      ELSE IF (ABS(DOTAC).GT.1.0D20) THEN
        WRITE(NOUT,*)'DOTAC messed up in AXES'
        STOP
      ELSE IF (ABS(DOTBC).GT.1.0D20) THEN
        WRITE(NOUT,*)'DOTBC messed up in AXES'
        STOP
      ENDIF
C
      RETURN
      END
C
C23456789012345678901234567890123456789012345678901234567890123456789012
C                 Subroutine to calculate integrand                    C
C23456789012345678901234567890123456789012345678901234567890123456789012
C
      SUBROUTINE CALCULATE(VRTX,SELFPROP,GRAPHNUMBER,KIN,ABSKIN,
     >               QS,QSIGNS,MAPTYPES,NMAPS,VALUE,MAXPART,VALUECHK)
C
      INTEGER SIZE,MAXMAPS
      PARAMETER (SIZE = 3)
      PARAMETER (MAXMAPS = 64)
C In:
      INTEGER VRTX(0:3*SIZE-1,2)
      LOGICAL SELFPROP(3*SIZE-1)
      INTEGER GRAPHNUMBER
      REAL*8 KIN(0:3*SIZE-1,0:3),ABSKIN(0:3*SIZE-1)
      INTEGER QS(MAXMAPS,0:SIZE),QSIGNS(MAXMAPS,0:SIZE)
      CHARACTER*6 MAPTYPES(MAXMAPS)
      INTEGER NMAPS
C Out:
      COMPLEX*16 VALUE,VALUECHK
      REAL*8 MAXPART
C
C Calculates the value of the graph specified by VRTX at the point K,
C returning result in VALUE, which includes the division by the density
C of points and the jacobian for deforming the contour. Also reports
C MAXPART, the biggest absolute value of the contributions to Re(VALUE).
C This helps us to keep track of cancellations and thus to abort the 
C calculation if too much cancellation among terms will be required.
C
C*********************
C
C Max number of graphs, for array size:
      INTEGER MAXGRAPHS
      PARAMETER (MAXGRAPHS = 10)
C
      INTEGER NIN,NOUT
      COMMON /IOUNIT/ NIN,NOUT
      INTEGER NLOOPS,NPROPS,NVERTS,CUTMAX
      COMMON /SIZES/ NLOOPS,NPROPS,NVERTS,CUTMAX
      REAL*8 MUOVERRTS
      COMMON /RENORMALIZE/ MUOVERRTS
      LOGICAL REPORT,DETAILS
      COMMON /CALCULOOK/ REPORT,DETAILS
      REAL*8 NC,NF
      COMMON /COLORFACTORS/ NC,NF
C How many graphs and how many cuts and maps for each:
      INTEGER NUMBEROFGRAPHS
      INTEGER NUMBEROFCUTS(MAXGRAPHS)
      INTEGER NUMBEROFMAPS(MAXGRAPHS)
      COMMON /GRAPHCOUNTS/ NUMBEROFGRAPHS,NUMBEROFCUTS,NUMBEROFMAPS
C Labels:
      INTEGER QE(0:SIZE)
C Momenta:
      REAL*8 K(0:3*SIZE-1,0:3),ABSK(0:3*SIZE-1)
      REAL*8 KINLOOP(SIZE+1,0:3)
      REAL*8 KCUT(SIZE+1,0:3)
      COMPLEX*16 NEWKINLOOP(0:3)
      COMPLEX*16 KC(0:3*SIZE-1,0:3)
      COMPLEX*16 ELLSQ,ELL
      REAL*8 E(0:SIZE),RTS
C Renormalization:
      REAL*8 MUMSBAR
C Matrices:
      INTEGER AE(0:3*SIZE-1,0:SIZE)
C FINDA variable:
      LOGICAL QOK
C DENSITY variables:
      REAL*8 JACNEWPOINT,DENSITY
C NEWCUT variables:
      INTEGER CUTINDEX(SIZE+1),CUTSIGN(SIZE+1),NCUT
      INTEGER ISIGN(3*SIZE-1)
      LOGICAL LEFTLOOP,RIGHTLOOP
      INTEGER LOOPINDEX(SIZE+1),LOOPSIGN(SIZE+1),NINLOOP
C Loopcut variables:
      LOGICAL CALCMORE
      INTEGER JCUT
C DEFORM variables:
      COMPLEX*16 JACDEFORM
C Functions:
      REAL*8 CALS0,SMEAR
      REAL*8 XXREAL,XXIMAG
      COMPLEX*16 COMPLEXSQRT
C Index variables:
      INTEGER P,MU,I,J,CUTNUMBER
C Propagator properties and momenta
      LOGICAL INLOOP(3*SIZE-1),CUT(3*SIZE-1)
      LOGICAL LOOPCUT(3*SIZE-1)
      INTEGER CUTSIGNP(3*SIZE-1)
      INTEGER LOOPLABEL(3*SIZE-1)
      REAL*8 KSQ
      COMPLEX*16 KCSQ
C Results variables:
      REAL*8 CALSVAL
      REAL*8 WEIGHT,MAXWEIGHT
      COMPLEX*16 NUMERATOR,RNUMERATOR,FEYNMAN
      COMPLEX*16 LOOPDENOM
      REAL*8 PLAINDENOM
      REAL*8 PREFACTOR
      COMPLEX*16 INTEGRAND
      COMPLEX*16 CHECK
C Useful constants:
      REAL*8 METRIC(0:3)
      DATA METRIC /+1.0D0,-1.0D0,-1.0D0,-1.0D0/
      REAL*8 PI
      DATA PI /3.1415926535898D0/
C      
C-----------------------------------------------------------------------
C Latest revision: 11 May 1996
C                  24 October 1996 (call to CHECKDEFORM)
C                  15 November 1996 (remove finite 'i epsilon')
C                  18 November 1996 (add CHECKVALUE)
C                  22 November 1996 Bug fixed.
C                  27 November 1996 (complex checkvalue)
C                  29 November 1996 (branchcut check; better checkvalue)
C                  27 February 1997 renormalization; reporting
C                  25 July     1997 renormalization; self-energy graphs
C                  17 September 1997 more renormalization & self-energy
C                  21 September 1997 finish DEFORM
C                  24 September 1997 fix bugs
C                  19 October 1997 fix cutsign bug
C                  22 October 1997 fix renormalization sign bug
C                   6 November 1997 improvements for deformation
C                  28 November 1997 more work on deformation
C                   2 December 1997 more precision in "report" numbers
C                   4 January  1998 revisions for self-energy graphs
C                  11 January 1998 renormalizaion for self-energy graphs
C                  27 February 1998 use Hrothgar for output
C                   5 March 1998 integrate Hrothgar
C                  14 March 1998 restore checks of deformation direction
C                  24 July 1998 use countfactor(graphnumber)
C                   4 August 1998 better CHECKDEFORM
C                   5 August 1998 change to groupsize(graphnumber)
C                  22 August 1998 add color factors
C                  22 December 1998 precalculate cut structure in RENO
C                  26 April 1999 omit REFLECT except as option
C                  22 December 2000 omit REFLECT entirely
C                  22 December 2000 change method of choosing points
C----------------------------------
C
C We do not want to change the value of KIN and ABSKIN, even though
C K and ABSK get changed by the reflection feature of the subroutine.
C
      DO P = 1,NPROPS
        ABSK(P) = ABSKIN(P)
        DO MU = 0,3
          K(P,MU) = KIN(P,MU)
        ENDDO
      ENDDO
C
C Initialize contribution to integral from this point. Also initialize
C BIGGEST, which will be the biggest absolute value of the contributions
C to VALUE. This helps us to keep track of cancellations and thus to
C abort the calculation if too much cancellation among terms will be
C required.
C
      MAXPART = 0.0D0
      VALUE = (0.0D0,0.0D0)
      VALUECHK = (0.0D0,0.0D0)
C
C Calculate jacobian.
C
      JACNEWPOINT = 
     >     1.0D0/DENSITY(GRAPHNUMBER,K,QS,QSIGNS,MAPTYPES,NMAPS)
C
C Loop over cuts.
C
      DO CUTNUMBER = 1,NUMBEROFCUTS(GRAPHNUMBER)
      CALL GETCUTINFO(GRAPHNUMBER,CUTNUMBER,NCUT,ISIGN,
     >            CUTINDEX,CUTSIGN,LEFTLOOP,RIGHTLOOP,
     >            NINLOOP,LOOPINDEX,LOOPSIGN)
C....
      IF (REPORT) THEN
        WRITE(NOUT,301)NCUT,CUTINDEX(1),CUTINDEX(2),
     >                    CUTINDEX(3),CUTINDEX(4)
301     FORMAT('Ncut =',I2,' CUTINDEX =',4I2)
      ENDIF
C''''
C
C Calculate Sqrt(s) and the renormalization scale MUMSBAR.
C
      RTS = 0.0D0
      DO J=1,NCUT
        RTS = RTS + ABSK(CUTINDEX(J))
      ENDDO
      MUMSBAR = MUOVERRTS * RTS
C
C Calculate final state momenta.
C Then we can also calculate CALSVAL and the PREFACTOR.
C
      DO I = 1,NCUT
        KCUT(I,0) = ABSK(CUTINDEX(I))
        DO MU = 1,3
          KCUT(I,MU) = CUTSIGN(I) * K(CUTINDEX(I),MU)
        ENDDO
      ENDDO
      CALSVAL = CALS0(NCUT,KCUT)
      PREFACTOR = 1.0D0 / (NC * RTS**2 * (2.0D0 * PI)**NLOOPS )
C
C Calculate momenta around loop (if any). In case NINLOOP = 0, this
C DO loop is skipped.
C
      DO J = 1,NINLOOP
        DO MU = 1,3
         KINLOOP(J,MU) = LOOPSIGN(J) * K(LOOPINDEX(J),MU)
        ENDDO
      ENDDO
C
C Please note that at this point the energy in the loop, KINLOOP(J,0),
C is not calculated. We have to wait until we have a loop cut to
C do this.
C
C Now KINLOOP(J,MU) gets an imaginary part for MU = 1,2,3.
C DEFORM calculates NEWKINLOOP and the associated jacobian, JACDEFORM.
C In case NINLOOP = 0, this subroutine just returns NEWKINLOOP(MU) = 0
C and JACDEFORM = 1.
C
      CALL DEFORM(VRTX,LOOPINDEX,RTS,LEFTLOOP,RIGHTLOOP,
     >      NINLOOP,KINLOOP,NEWKINLOOP,JACDEFORM)
C
C If there is a loop, we need to go around the loop and generate
C a "loopcut."  There are three cases.
C 1) NINLOOP = 0, with NCUT = CUTMAX.  
C Then we are ready to proceed, and we should calculate only once 
C before going back to NEWCUT. Therefore we set CALCMORE to .FALSE. 
C so that we do not enter this code again.
C 2) NINLOOP = 2, with NCUT = CUTMAX - 1.
C Then the loop is a self-energy subgraph and, with our dispersive
C treatment of these graphs, there is one term, with JCUT = 1.
C 3) NINLOOP > 2, with NCUT = CUTMAX - 1  
C Then we should loop over JCUT = 1,2,...,NINLOOP and
C set CUTINDEX(CUTMAX) = LOOPINDEX(JCUT).  When we are done with this
C we set CALCMORE to .FALSE. .
C
C If we need to renormalize this loop, we will do it
C as part of the JCUT = 1 calculation.
C
C We initialize the weight, then add to it for the renormalization
C counterterm and for each loopcut.
C
      WEIGHT = 0.0D0
      MAXWEIGHT = 0.0D0
C
      JCUT = 0
      CALCMORE = .TRUE.
      DO WHILE (CALCMORE)
      IF (NINLOOP.EQ.0) THEN
         CALCMORE = .FALSE.
      ELSE
         JCUT = JCUT + 1
         CUTINDEX(CUTMAX) = LOOPINDEX(JCUT)
         CUTSIGN(CUTMAX) = LOOPSIGN(JCUT)
         IF ((NINLOOP.EQ.2).OR.(JCUT.EQ.NINLOOP)) THEN
             CALCMORE = .FALSE.
         ENDIF
      ENDIF
C
C Calculate matrix AE(P,I) relating propagator energies to energies of 
C cut lines. NOTE that the index I here is displaced by 1.
C
      DO I = 0,NLOOPS
        QE(I) = CUTINDEX(I+1)
      ENDDO
      CALL FINDA(VRTX,QE,NLOOPS,AE,QOK)
      IF (.NOT.QOK) THEN
         WRITE(NOUT,*)'AE not found'
         STOP
      ENDIF
C
C Find which propagators are which. A propagator can be exactly
C on shell even if it isn't cut if it is linked by a self-energy
C correction to a propagator that is cut. The matrix AE(P,I) will
C tell us.
C
C Define logical and sign variables:
C CUT(P) = .TRUE. if propagator P crosses the final state cut.
C CUTSIGNP(P) = CUTSIGN(I) for P = CUTINDEX(I).
C LOOPCUT(P) = .TRUE. propagator P crosses the loopcut.
C INLOOP(P) = .TRUE. if it is in a virtual loop.
C LOOPLABEL(P) = label 1,2,... counting around loop for propagator P in 
C                 loop.
C SELFPROP(P) = .TRUE. if it is part of a one loop self-energy diagram
C                      or attaches to such a diagram.
C
      DO P = 1,NPROPS
         CUT(P) = .FALSE.
         LOOPCUT(P) = .FALSE.
         INLOOP(P) = .FALSE.
         CUTSIGNP(P) = 0
         LOOPLABEL(P) = 0
      ENDDO
      DO I = 1,CUTMAX
         CUT(CUTINDEX(I)) = .TRUE.
         CUTSIGNP(CUTINDEX(I)) = CUTSIGN(I)
      ENDDO
      IF (NINLOOP.GT.0) THEN
         CUT(CUTINDEX(CUTMAX)) = .FALSE.
         LOOPCUT(CUTINDEX(CUTMAX)) = .TRUE.
      ENDIF
      DO J = 1,NINLOOP
         INLOOP(LOOPINDEX(J)) = .TRUE.
         LOOPLABEL(LOOPINDEX(J)) = J
      ENDDO
C
C Calculate part of the energies of cut lines corresponding to the
C real part of the loop three-momenta. NOTE that I is displaced by 1
C in order to work with the matrix AE(P,I).
C
      DO I = 0,NLOOPS
        E(I) = CUTSIGN(I+1) * ABSK(CUTINDEX(I+1))
      ENDDO
C
C Calculate part of the propagator energies corresponding to the
C real part of the loop three-momenta.
C
      DO P = 0,NPROPS
         K(P,0) = 0.0D0
         DO I = 0,NLOOPS
           K(P,0) = K(P,0) + AE(P,I) * E(I)
         ENDDO
      ENDDO
      IF ( ABS(RTS - K(0,0)).GT.1.0D-8 ) THEN
         WRITE(NOUT,*)'Oops, the calculation of RTS did not work'
         STOP
      ENDIF
C
C Calculate the added complex loop energy. Check that we do not
C cross the cut of Sqrt(ELLSQ) by using COMPLEXSQRT(ELLSQ).
C
      IF (NINLOOP.GT.0) THEN
        KINLOOP(JCUT,0) =  LOOPSIGN(JCUT) * K(LOOPINDEX(JCUT),0)
        ELLSQ = (0.0D0,0.0D0)
        DO MU = 1,3
          ELLSQ = ELLSQ + ( KINLOOP(JCUT,MU) + NEWKINLOOP(MU) )**2
        ENDDO
        ELL = COMPLEXSQRT(ELLSQ)
        NEWKINLOOP(0) = ELL - KINLOOP(JCUT,0)
      ELSE
        NEWKINLOOP(0) = (0.0D0,0.0D0)
      ENDIF
C....
      IF (REPORT) THEN
        IF( DETAILS .AND. (NINLOOP.GT.0) ) THEN
          WRITE(NOUT,340)NEWKINLOOP(0),XXIMAG(NEWKINLOOP(1)),
     >                 XXIMAG(NEWKINLOOP(2)),XXIMAG(NEWKINLOOP(3))
340       FORMAT('NEWKINLOOP =',2(1P G12.3),' AND',3(1P G12.3))
        ENDIF
      ENDIF
C''''
C Now we calculate the complex propagator momenta.
C
      DO P = 0,NPROPS
        DO MU = 0,3
         KC(P,MU) = K(P,MU)
        ENDDO
      ENDDO
C
      DO J = 1,NINLOOP
       DO MU = 0,3
         KC(LOOPINDEX(J),MU) = KC(LOOPINDEX(J),MU)
     >       + LOOPSIGN(J) * NEWKINLOOP(MU)
        ENDDO
      ENDDO
C
C Calculate denominator.
C
C We calculate two denominators: the part from the loop propagators
C (LOOPDENOM) and the part from the other propagators (PLAINDENOM).
C The renormalization counterterm uses only PLAINDENOM.
C
C Propagators that are part of a one loop self-energy subgraph, or
C attached to such a subgraph, do not contribute to the denominator
C factor at all. The functions QPROP and GPROP, called by NUMERATOR,
C take care of the factors associated with these propagators.
C
      PLAINDENOM = 1.0D0
      LOOPDENOM = (1.0D0,0.0D0)
      DO P = 1,NPROPS
C
      IF (.NOT.SELFPROP(P)) THEN
        IF (INLOOP(P)) THEN
C
C P is in the loop:
C
          IF(LOOPCUT(P)) THEN
            LOOPDENOM   = LOOPDENOM * 2.0D0 * CUTSIGNP(P) * KC(P,0)
          ELSE
            KCSQ = 0.0D0
            DO MU = 0,3
              KCSQ = KCSQ + METRIC(MU) * KC(P,MU)**2
            ENDDO
            CALL CHECKDEFORM(KCSQ,LEFTLOOP,RIGHTLOOP,LOOPLABEL(P),JCUT,
     >                       GRAPHNUMBER,CUT,LOOPCUT)
            LOOPDENOM = LOOPDENOM * KCSQ
          ENDIF
C
        ELSE
C
C P is not in the loop:
C
          IF (CUT(P)) THEN
            PLAINDENOM = PLAINDENOM * 2.0D0 * CUTSIGNP(P) * K(P,0)
          ELSE
            KSQ = 0.0D0
            DO MU = 0,3
              KSQ = KSQ + METRIC(MU) * K(P,MU)**2
            ENDDO
            PLAINDENOM = PLAINDENOM * KSQ
          ENDIF
C
C End   IF (INLOOP(P)) ... ELSE ...
C End IF (.NOT.SELFPROP(P)) ...
C
        ENDIF
      ENDIF
C........ 
        IF (REPORT.AND.DETAILS) THEN
        IF (.NOT.SELFPROP(P)) THEN
          IF (INLOOP(P)) THEN
            IF(LOOPCUT(P)) THEN
              WRITE(NOUT,350)P,CUTSIGNP(P) * KC(P,0)
350           FORMAT('Loopcut propagator',I3,' Energy =',2(1P G12.3))
            ELSE
              WRITE(NOUT,351)P,KCSQ
351           FORMAT('Loop propagator',I3,' KCSQ =',2(1P G12.3))
            ENDIF
          ELSE
            IF (CUT(P)) THEN
              WRITE(NOUT,352)P,CUTSIGNP(P) * K(P,0)
352           FORMAT('Cut propagator',I3,' Energy =',(1P G12.3))
            ELSE
              WRITE(NOUT,353)P,KSQ
353           FORMAT('Tree propagator',I3,' KSQ =',(1P G12.3))
            ENDIF
          ENDIF
        ENDIF
        ENDIF
C'''''''
C
C End DO P = 1,NPROPS
C
      ENDDO
C
C Calculate graph.
C
C Add to contribution for this point.
C
C If we have a virtual loop with 2 or 3 lines, then we need the
C renormalization counter term. We are in a loop over JCUT. We
C will include the counter term when JCUT = 1.
C
C There is a minus sign because this is the counter term and we
C want to subtract it.
C      
      IF (((NINLOOP.EQ.2).OR.(NINLOOP.EQ.3)).AND.(JCUT.EQ.1)) THEN
C
        FEYNMAN = RNUMERATOR(GRAPHNUMBER,KC,MUMSBAR,CUT)/PLAINDENOM
C
        INTEGRAND = - PREFACTOR * JACNEWPOINT * JACDEFORM 
     >              * FEYNMAN * SMEAR(RTS)
        MAXWEIGHT = MAX(MAXWEIGHT,ABS(XXREAL(INTEGRAND)))
        WEIGHT = WEIGHT + XXREAL(INTEGRAND)
C
        INTEGRAND = INTEGRAND * CALSVAL
        MAXPART = MAX(MAXPART,ABS(XXREAL(INTEGRAND)))
        VALUE = VALUE + INTEGRAND
C....
        IF (REPORT) THEN
         IF (DETAILS) THEN
          WRITE(NOUT,360)
360       FORMAT('PREFACTOR * JACNEWPOINT * (JACDEFORM-R JACDEFORM-I)',
     >           ' (FEYNMAN-R FEYNMAN-I) * CALSVAL * SMEAR(RTS)')
          WRITE(NOUT,361)PREFACTOR,JACNEWPOINT,JACDEFORM,
     >                 FEYNMAN,CALSVAL,SMEAR(RTS)
361       FORMAT(8(1P G12.3))
         ENDIF
         WRITE(NOUT,362)INTEGRAND
362      FORMAT('Contribution (CT):',2(1P G18.10))
         IF (DETAILS) THEN
           WRITE(NOUT,*)' '
         ENDIF
        ENDIF
C''''
      ENDIF
C
C Done with the counter term (if any), now we do the main term.
C
      FEYNMAN = NUMERATOR(GRAPHNUMBER,KC,CUT)/LOOPDENOM /PLAINDENOM
      INTEGRAND = PREFACTOR * JACNEWPOINT * JACDEFORM 
     >            * FEYNMAN * SMEAR(RTS)
      MAXWEIGHT = MAX(MAXWEIGHT,ABS(XXREAL(INTEGRAND)))
      WEIGHT = WEIGHT + XXREAL(INTEGRAND)
C
      INTEGRAND = INTEGRAND * CALSVAL
      MAXPART = MAX(MAXPART,ABS(XXREAL(INTEGRAND)))
      VALUE = VALUE + INTEGRAND
C....
      IF (REPORT) THEN
        IF (DETAILS) THEN
          WRITE(NOUT,370)
370       FORMAT('PREFACTOR * JACNEWPOINT * (JACDEFORM-R JACDEFORM-I)',
     >           ' (FEYNMAN-R FEYNMAN-I) * CALSVAL * SMEAR(RTS)')
          WRITE(NOUT,371)PREFACTOR,JACNEWPOINT,JACDEFORM,
     >                 FEYNMAN,CALSVAL,SMEAR(RTS)
371       FORMAT(8(1P G12.3))
        ENDIF
        IF (NINLOOP.GT.0) THEN
          WRITE(NOUT,373)LOOPINDEX(JCUT),INTEGRAND
373       FORMAT(I3,'  Contribution:',2(1P G18.10))
        ELSE
          WRITE(NOUT,374)INTEGRAND
374       FORMAT('     Contribution:',2(1P G18.10))        
        ENDIF
        IF (DETAILS) THEN
          WRITE(NOUT,*)' '
        ENDIF
      ENDIF
C''''
C
C Compute a known integral to see if we have it right.
C Subroutine CHECKCALC calculates CHECK.
C
      CALL 
     >  CHECKCALC(GRAPHNUMBER,CUTINDEX,KC,JACNEWPOINT,JACDEFORM,CHECK)
      VALUECHK = VALUECHK + CHECK
C
C End of loop DO WHILE (CALCMORE) that runs over loopcuts.
C
      ENDDO
C
C We are ready to call Hrothgar to process the result for this cut.
C
      CALL HROTHGAR(NCUT,KCUT,WEIGHT,1,'NEWRESULT ')
C
C Close loop DO CUTNUMBER = 1,NUMBEROFCUTS(GRAPHNUMBER)
C
      ENDDO
C
      RETURN
      END
C
C23456789012345678901234567890123456789012345678901234567890123456789012
C             End of subroutine to calculate integrand                 C
C23456789012345678901234567890123456789012345678901234567890123456789012
C
      SUBROUTINE 
     >  CHECKCALC(GRAPHNUMBER,CUTINDEX,KC,JACNEWPOINT,JACDEFORM,CHECK)
C
      INTEGER SIZE
      PARAMETER (SIZE = 3)
C In:
      INTEGER GRAPHNUMBER,CUTINDEX(SIZE+1)
      COMPLEX*16 KC(0:3*SIZE-1,0:3)
      REAL*8 JACNEWPOINT
      COMPLEX*16 JACDEFORM
C Out:
      COMPLEX*16 CHECK
C
C Compute a known integral to see if we have it right.
C This subroutine calculates the integrand.
C The check is based on 
C   Int d^3 p [p^2 + M^2]^(-3) = Pi^2/ (4 M^3).
C   Int d^3 p [p^2 (p^2 + M^2)]^(-1) = 2 Pi^2 /M
C Note that we look at just one term in the sum over cuts
C and loopcuts: 
C For graph 10, we take Cutindex = (7,5,4,1);
C For graph 8, we take Cutindex = (8,6,4,1), etc.
C
C Latest modification: 24 December 1998
C
C Max number of graphs, cuts, maps for array sizes:
      INTEGER MAXGRAPHS,MAXMAPS
      PARAMETER (MAXGRAPHS = 10)
      PARAMETER (MAXMAPS = 64)
C Input and output units.
      INTEGER NIN,NOUT
      COMMON /IOUNIT/ NIN,NOUT
C How many graphs and how many cuts and maps for each:
      INTEGER NUMBEROFGRAPHS
      INTEGER NUMBEROFCUTS(MAXGRAPHS)
      INTEGER NUMBEROFMAPS(MAXGRAPHS)
      COMMON /GRAPHCOUNTS/ NUMBEROFGRAPHS,NUMBEROFCUTS,NUMBEROFMAPS
C Reno size and counting variables:
      INTEGER GROUPSIZE(MAXGRAPHS,MAXMAPS)
      INTEGER GROUPSIZEGRAPH(MAXGRAPHS)
      INTEGER GROUPSIZETOTAL
      COMMON /MONTECARLO/GROUPSIZE,GROUPSIZEGRAPH,GROUPSIZETOTAL
C
      REAL*8 MM
      DATA MM /3.0D-1/
      REAL*8 PI
      DATA PI /3.1415926535898D0/
      COMPLEX*16 TEMP1,TEMP2,TEMP3
      INTEGER MU
C
C If it is not the right graph and the right cut, this default
C value will be returned.
C
      CHECK = (0.0D0,0.0D0)
C
      TEMP1 = 0.0D0
      TEMP2 = 0.0D0
      TEMP3 = 0.0D0
C
      IF (GRAPHNUMBER.EQ.10) THEN
C
        IF (   (CUTINDEX(1).EQ.7).AND.(CUTINDEX(2).EQ.5)
     >    .AND.(CUTINDEX(3).EQ.4).AND.(CUTINDEX(4).EQ.1) ) THEN
          DO MU = 1,3
            TEMP1 = TEMP1 + KC(7,MU)*KC(7,MU)
            TEMP2 = TEMP2 + KC(6,MU)*KC(6,MU)
            TEMP3 = TEMP3 + KC(1,MU)*KC(1,MU) 
          ENDDO
        ELSE
          RETURN
        ENDIF
C
      ELSE IF (GRAPHNUMBER.EQ.9) THEN
C
        IF (   (CUTINDEX(1).EQ.8).AND.(CUTINDEX(2).EQ.7)
     >    .AND.(CUTINDEX(3).EQ.3).AND.(CUTINDEX(4).EQ.5) ) THEN
          DO MU = 1,3
            TEMP1 = TEMP1 + KC(8,MU)*KC(8,MU)
            TEMP2 = TEMP2 + KC(6,MU)*KC(6,MU)
            TEMP3 = TEMP3 + KC(1,MU)*KC(1,MU) 
          ENDDO
        ELSE
          RETURN
        ENDIF
C      
      ELSE IF (GRAPHNUMBER.EQ.8) THEN
C
        IF (   (CUTINDEX(1).EQ.8).AND.(CUTINDEX(2).EQ.6)
     >    .AND.(CUTINDEX(3).EQ.4).AND.(CUTINDEX(4).EQ.1) ) THEN
          DO MU = 1,3
            TEMP1 = TEMP1 + KC(5,MU)*KC(5,MU)
            TEMP2 = TEMP2 + KC(8,MU)*KC(8,MU)
            TEMP3 = TEMP3 + KC(1,MU)*KC(1,MU) 
          ENDDO
        ELSE
          RETURN
        ENDIF
C
      ELSE IF (GRAPHNUMBER.EQ.7) THEN
C
        IF (   (CUTINDEX(1).EQ.5).AND.(CUTINDEX(2).EQ.4)
     >    .AND.(CUTINDEX(3).EQ.1).AND.(CUTINDEX(4).EQ.6) ) THEN
          DO MU = 1,3
            TEMP1 = TEMP1 + KC(5,MU)*KC(5,MU)
            TEMP2 = TEMP2 + KC(7,MU)*KC(7,MU)
            TEMP3 = TEMP3 + KC(1,MU)*KC(1,MU) 
          ENDDO
        ELSE
          RETURN
        ENDIF
C
      ELSE IF (GRAPHNUMBER.EQ.6) THEN
C
        IF (   (CUTINDEX(1).EQ.8).AND.(CUTINDEX(2).EQ.5)
     >    .AND.(CUTINDEX(3).EQ.3).AND.(CUTINDEX(4).EQ.6) ) THEN
          DO MU = 1,3
            TEMP1 = TEMP1 + KC(5,MU)*KC(5,MU)
            TEMP2 = TEMP2 + KC(7,MU)*KC(7,MU)
            TEMP3 = TEMP3 + KC(1,MU)*KC(1,MU) 
          ENDDO
        ELSE
          RETURN
        ENDIF
C 
      ELSE IF (GRAPHNUMBER.EQ.5) THEN
C
        IF (   (CUTINDEX(1).EQ.8).AND.(CUTINDEX(2).EQ.7)
     >    .AND.(CUTINDEX(3).EQ.3).AND.(CUTINDEX(4).EQ.1) ) THEN
          DO MU = 1,3
            TEMP1 = TEMP1 + KC(5,MU)*KC(5,MU)
            TEMP2 = TEMP2 + KC(8,MU)*KC(8,MU)
            TEMP3 = TEMP3 + KC(1,MU)*KC(1,MU) 
          ENDDO
        ELSE
          RETURN
        ENDIF
C
      ELSE IF (GRAPHNUMBER.EQ.4) THEN
C
        IF (   (CUTINDEX(1).EQ.6).AND.(CUTINDEX(2).EQ.4)
     >    .AND.(CUTINDEX(3).EQ.1).AND.(CUTINDEX(4).EQ.7) ) THEN
          DO MU = 1,3
            TEMP1 = TEMP1 + KC(7,MU)*KC(7,MU)
            TEMP2 = TEMP2 + KC(5,MU)*KC(5,MU)
            TEMP3 = TEMP3 + KC(1,MU)*KC(1,MU) 
          ENDDO
        ELSE
          RETURN
        ENDIF
C 
      ELSE IF (GRAPHNUMBER.EQ.3) THEN
C
        IF (   (CUTINDEX(1).EQ.5).AND.(CUTINDEX(2).EQ.4)
     >    .AND.(CUTINDEX(3).EQ.1).AND.(CUTINDEX(4).EQ.6) ) THEN
          DO MU = 1,3
            TEMP1 = TEMP1 + KC(7,MU)*KC(7,MU)
            TEMP2 = TEMP2 + KC(5,MU)*KC(5,MU)
            TEMP3 = TEMP3 + KC(1,MU)*KC(1,MU) 
          ENDDO
        ELSE
          RETURN
        ENDIF
C 
      ELSE IF (GRAPHNUMBER.EQ.2) THEN
C
        IF (   (CUTINDEX(1).EQ.7).AND.(CUTINDEX(2).EQ.6)
     >    .AND.(CUTINDEX(3).EQ.1).AND.(CUTINDEX(4).EQ.4) ) THEN
          DO MU = 1,3
            TEMP1 = TEMP1 + KC(7,MU)*KC(7,MU)
            TEMP2 = TEMP2 + KC(5,MU)*KC(5,MU)
            TEMP3 = TEMP3 + KC(1,MU)*KC(1,MU) 
          ENDDO
        ELSE
          RETURN
        ENDIF
C
      ELSE IF (GRAPHNUMBER.EQ.1) THEN
C
        IF (   (CUTINDEX(1).EQ.5).AND.(CUTINDEX(2).EQ.4)
     >    .AND.(CUTINDEX(3).EQ.1).AND.(CUTINDEX(4).EQ.7) ) THEN
          DO MU = 1,3
            TEMP1 = TEMP1 + KC(7,MU)*KC(7,MU)
            TEMP2 = TEMP2 + KC(5,MU)*KC(5,MU)
            TEMP3 = TEMP3 + KC(1,MU)*KC(1,MU) 
          ENDDO
        ELSE
          RETURN
        ENDIF
C 
      ELSE
        WRITE(NOUT,*)'Problem with graph number in CHECKCALC'
        STOP
      ENDIF
C
C Here is an infrared sensitive check integral:
      CHECK =         TEMP1 * (TEMP1 + MM**2)
      CHECK = CHECK * TEMP2 * (TEMP2 + MM**2)
      CHECK = CHECK * (TEMP3 + MM**2)**3
      CHECK = (MM**5/PI**6) /CHECK
C Here is a nice smooth check integral:
C     CHECK =         (TEMP1 + MM**2)**3
C     CHECK = CHECK * (TEMP2 + MM**2)**3
C     CHECK = CHECK * (TEMP3 + (2.0D0*MM)**2)**3
C     CHECK = (512.0D0 * MM**9 / PI**6) /CHECK
C
      CHECK = JACDEFORM * JACNEWPOINT * CHECK
C
C Weight according to the number of points devoted to the current
C graph.
C
      CHECK = CHECK * GROUPSIZEGRAPH(GRAPHNUMBER)/GROUPSIZETOTAL
C
      RETURN
      END
C
C23456789012345678901234567890123456789012345678901234567890123456789012
C23456789012345678901234567890123456789012345678901234567890123456789012
C
      REAL*8 FUNCTION
     >    DENSITY(GRAPHNUMBER,K,QS,QSIGNS,MAPTYPES,NMAPS)
C
      INTEGER SIZE,MAXMAPS
      PARAMETER (SIZE = 3)
      PARAMETER (MAXMAPS = 64)
C In:
      INTEGER GRAPHNUMBER
      REAL*8 K(0:3*SIZE-1,0:3)
      INTEGER QS(MAXMAPS,0:SIZE),QSIGNS(MAXMAPS,0:SIZE)
      CHARACTER*6 MAPTYPES(MAXMAPS)
      INTEGER NMAPS
C
C Density of Monte Carlo points as a function of |K(p)|'s.
C 
C 29 June 1993
C 12 July 1993
C 17 July 1994
C  4 May 1996
C 21 November 1996
C  5 December 1996
C  5 February 1997
C 15 December 1998
C 23 December 1998
C  9 February 1999
C 10 March    1999
C 20 August   1999
C 21 December 2000
C 20 March    2001
C
      INTEGER NIN,NOUT
      COMMON /IOUNIT/ NIN,NOUT
      INTEGER NLOOPS,NPROPS,NVERTS,CUTMAX
      COMMON /SIZES/ NLOOPS,NPROPS,NVERTS,CUTMAX
      INTEGER MAXGRAPHS
      PARAMETER (MAXGRAPHS = 10)
      INTEGER GROUPSIZE(MAXGRAPHS,MAXMAPS)
      INTEGER GROUPSIZEGRAPH(MAXGRAPHS)
      INTEGER GROUPSIZETOTAL
      COMMON /MONTECARLO/GROUPSIZE,GROUPSIZEGRAPH,GROUPSIZETOTAL

C
      INTEGER MAPNUMBER,L,MU
      REAL*8 P1(3),P2(3),ELL1(3),ABSP1,ABSP2,ABSP3
      REAL*8 TEMP1,TEMP2,TEMP3,P1SQ,P2SQ,P3SQ
      CHARACTER*6 MAPTYPE
      INTEGER QSIGN(0:SIZE),Q(0:SIZE)
      REAL*8 RHO3,RHO2TO3D,RHO2TO3E,RHO2TO2T,RHO2TO2S,RHO2TO1
      REAL*8 RHOTHREE,RHOLOOP
C
C We construct the density as a sum.
C
      DENSITY = 0.0D0
      DO MAPNUMBER = 1,NMAPS
C
      MAPTYPE = MAPTYPES(MAPNUMBER)
      DO L = 0,NLOOPS
       Q(L) = QS(MAPNUMBER,L)
       QSIGN(L) = QSIGNS(MAPNUMBER,L)
      ENDDO
C
C First, we need the kinematic variables for this map.
C
      P1SQ = 0.0D0
      P2SQ = 0.0D0
      P3SQ = 0.0D0
      DO MU = 1,3
        ELL1(MU) = QSIGN(1)*K(Q(1),MU)
        TEMP1 = QSIGN(2)*K(Q(2),MU)
        TEMP2 = QSIGN(3)*K(Q(3),MU)
        TEMP3 = - TEMP1 - TEMP2
        P1(MU) = TEMP1
        P1SQ = P1SQ + TEMP1**2
        P2(MU) = TEMP2
        P2SQ = P2SQ + TEMP2**2
        P3SQ = P3SQ + TEMP3**2
      ENDDO
      ABSP1 = SQRT(P1SQ)
      ABSP2 = SQRT(P2SQ)
      ABSP3 = SQRT(P3SQ)
C
C Now, there are two factors, one for the 'final state momenta' and
C one for the 'loop momentum.'
C
      RHOTHREE = RHO3(ABSP1,ABSP2,ABSP3)
C
      IF (MAPTYPE.EQ.'T2TO3D') THEN
        RHOLOOP = RHO2TO3D(P1,P2,ELL1)
      ELSE IF (MAPTYPE.EQ.'T2TO3E') THEN
        RHOLOOP = RHO2TO3E(P1,P2,ELL1)
      ELSE IF (MAPTYPE.EQ.'T2TO2T') THEN
        RHOLOOP = RHO2TO2T(P1,P2,ELL1)
      ELSE IF (MAPTYPE.EQ.'T2TO2S') THEN
        RHOLOOP = RHO2TO2S(P1,P2,ELL1)
      ELSE IF (MAPTYPE.EQ.'T2TO1 ') THEN
        RHOLOOP = RHO2TO1(P1,P2,ELL1)
      ELSE 
        WRITE(NOUT,*)'Bad MAPTYPE in DENSITY'
        STOP
      ENDIF
C
      DENSITY = DENSITY
     >          + RHOTHREE*RHOLOOP*GROUPSIZE(GRAPHNUMBER,MAPNUMBER)
C
C Close  DO MAPNUMBER = 1,NMAPS 
C
      ENDDO
      RETURN
      END
C
C23456789012345678901234567890123456789012345678901234567890123456789012
C           Subroutines associated with NEWPOINT and DENSITY           2
C   CHOOSEx and RHOx where x = 3, 2to2T, 2to2S, 2to3D, 2to3E, 2to1     2
C23456789012345678901234567890123456789012345678901234567890123456789012
C
      SUBROUTINE CHOOSE3(P1,P2,P3,OK)
C
C Out:
      REAL*8 P1(3),P2(3),P3(3)
      LOGICAL OK
C
C Generates momenta P1(mu),P2(mu),P3(mu) for a three body final 
C state with a distribution in momentum fractions x1,x2,x3 
C proportional to
C
C        [max(1-x1,1-x2,1-x3)]^B/[(1-x1)*(1-x2)*(1-x3)]^B.
C
C 28 December 2000
C 16 January 2001
C
      REAL*8 BADNESSLIMIT,CANCELLIMIT,THRUSTCUT
      COMMON /LIMITS/ BADNESSLIMIT,CANCELLIMIT,THRUSTCUT
      REAL*8 ONETHIRD,TWOTHIRDS,PI
      PARAMETER(ONETHIRD  = 0.3333333333333333333D0)
      PARAMETER(TWOTHIRDS = 0.6666666666666666667D0)
      PARAMETER (PI = 3.141592653589793239D0)
C
C The parameter E3PAR should match between CHOOSE3 and RHO3.
C
      REAL*8 E3PAR
      PARAMETER(E3PAR = 1.5D0)
C
C The parameters A, B, and C need to match between CHOOSE3 and RHO3.
C CHOOSE3 uses A, while RHO3 uses B and C. The relation is
C B = 1 - 1/A and then C is the normalization factor and is
C a rather complicated function of B.
C
C Some soft and collinear points:
C
       REAL*8 A,B,C
       PARAMETER(A = 2.0D0)
       PARAMETER(B = 0.5D0)
       PARAMETER(C = 0.0036376552621307193655D0)
C
C Lots of soft and collinear points:
C
C      REAL*8 A,B,C
C      PARAMETER(A = 4.0D0)
C      PARAMETER(B = 0.75D0)
C      PARAMETER(C = 0.00058417226323428314253D0)
C
      REAL*8 X,RANDOM
      LOGICAL DONE
      INTEGER MU
      REAL*8 EMAX
      REAL*8 X1,X2,X3,Y1,Y2,Y3
      REAL*8 EA(3),EB(3),EC(3),ED(3)
      REAL*8 PHI,COSTHETA,SINTHETA
      REAL*8 K1(3),K2(3),K3(3)
C
C----------
C
      OK = .TRUE.
C
C We will generate vectors K1(mu), K2(mu), K3(mu) with |K1| > |K3| and
C |K2| > |K3|. At the end, we will associate each Ki(mu) with a Pj(mu)
C with the index j of the Pj(mu) that matches K3(mu) chosen at random.
C 
C We choose y1, y2, y3 in 0< y_i < 1 with y1 + y2 + y3 = 1. The y_i are
C related to the momentum fractions x_i by y_i = 1 - x_i. For the y_i,
C we want y3 to be the largest, with no specification about whether y1
C or y2 is larger.  We want to choose y1 and y2 with a 1/sqrt(y1*y2)
C distribution. Then y3 = 1 - y1 - y2.  We must insure that y3 > y1 and
C y3 > y2 for the point to be valid. Note that the allowed region is
C inside the region 0 < y1 < 1/2, 0 < y2 < 1/2. If we choose a random
C variable x in 0 < x < 1 and define y = x**2/2 then the density dx/dy
C is proportional to 1/sqrt(y) and 0 < y < 1/2.
C
C We loop until we are "done" choosing a valid point.
C
      DONE = .FALSE.
      DO WHILE (.NOT.DONE)
        X = RANDOM(1)
        Y1 = 0.5D0 * X**A
        X = RANDOM(1)
        Y2 = 0.5D0 * X**A
        Y3 = 1.0D0 - Y1 - Y2
        IF ((Y1 .LT. Y3).AND.(Y2.LT.Y3)) THEN
          DONE = .TRUE.
        ENDIF
      ENDDO
      X1 = 1.0D0 - Y1
      X2 = 1.0D0 - Y2
      X3 = Y1 + Y2
C
C If the chosen point is too soft or collinear, we will not be able
C to compute the kinematics for the rest of this subroutine
C or the other CHOOSEx subroutines, so we just abort.
C
      IF ( Y1*Y2.LT.(100.0D0*BADNESSLIMIT)**(-2) ) THEN
        DO MU = 1,3
          P1(MU) = 0.0D0
          P2(MU) = 0.0D0
          P3(MU) = 0.0D0
        ENDDO
        OK = .FALSE.
        RETURN
      ENDIF
C
C Choose Emax = sum_i |p_i| /2.
C
      X = RANDOM(1)
      EMAX = E3PAR * ( 1.0D0/X - 1.0D0 )**ONETHIRD
C
C Choose a direction EA(mu) at random on the unit sphere.
C
      X = RANDOM(1)
      COSTHETA = 2.0D0*X - 1.0D0
      SINTHETA = SQRT(1.0D0 - COSTHETA**2)
      X = RANDOM(1)
      PHI = 2.0D0 * PI * X
      EA(1) = SINTHETA * COS(PHI)
      EA(2) = SINTHETA * SIN(PHI)
      EA(3) = COSTHETA
C
C Generate vectors EB and EC that form a right handed basis set with EA.
C
      CALL AXES(EA,EB,EC)
C
C Generate a unit vector ED at a with a random azimuthal angle around
C the EA axis in this basis.
C
      X = RANDOM(1)
      PHI = 2.0D0 * PI * X
      DO MU = 1,3
        ED(MU) = COS(PHI)*EB(MU) + SIN(PHI)*EC(MU)
      ENDDO
C
C Now construct the momenta. P1(mu) is directed in the random direction
C EA(mu) with magnitude determined from Emax and X1. Then P3(mu)
C is in the plane of EA(mu) and ED(mu) with angle THETA to P1(mu) 
C determined from the Xi and magnitude determined by X2.
C
      COSTHETA = 1.0D0 - 2.0D0*Y2/X1/X3
      SINTHETA = 2.0D0*SQRT(Y1*Y2*Y3)/X1/X3
      DO MU = 1,3
        K1(MU) = X1*EMAX*EA(MU)
        K3(MU) = X3*EMAX*(COSTHETA*EA(MU) + SINTHETA*ED(MU))
        K2(MU) = - K1(MU) - K3(MU)
      ENDDO
C
C Match K3(mu) to one of the Pi(mu) at random.
C
      X = RANDOM(1)
      IF (X.GT.TWOTHIRDS) THEN
        DO MU = 1,3
          P1(MU) = K1(MU)
          P2(MU) = K2(MU)
          P3(MU) = K3(MU)
        ENDDO
      ELSE IF (X.GT.ONETHIRD) THEN
        DO MU = 1,3
          P1(MU) = K2(MU)
          P2(MU) = K3(MU)
          P3(MU) = K1(MU)
        ENDDO
      ELSE
        DO MU = 1,3
          P1(MU) = K3(MU)
          P2(MU) = K1(MU)
          P3(MU) = K2(MU)
        ENDDO
      ENDIF
C
      RETURN
      END
C
C23456789012345678901234567890123456789012345678901234567890123456789012
C23456789012345678901234567890123456789012345678901234567890123456789012
C
      REAL*8 FUNCTION RHO3(ABSP1,ABSP2,ABSP3)
C
C In:
      REAL*8 ABSP1,ABSP2,ABSP3
C
      INTEGER NIN,NOUT
      COMMON /IOUNIT/ NIN,NOUT
C
C Density of points for points chosen with CHOOSE3(p1,p2,p3,ok).
C 16 January 2001
C
      REAL*8 EMAX,X1,X2,X3
      REAL*8 E03,EMAX3,FACTOR,DENOM
C
C The parameter E3PAR should match between CHOOSE3 and RHO3.
C 
      REAL*8 E3PAR
      PARAMETER(E3PAR = 1.5D0)
C
C The parameters A, B, and C need to match between CHOOSE3 and RHO3.
C CHOOSE3 uses A, while RHO3 uses B and C. The relation is
C B = 1 - 1/A and then C is the normalization factor and is
C a rather complicated function of B.
C
C Some soft and collinear points:
C
       REAL*8 A,B,C
       PARAMETER(A = 2.0D0)
       PARAMETER(B = 0.5D0)
       PARAMETER(C = 0.0036376552621307193655D0)
C
C Lots of soft and collinear points:
C
C      REAL*8 A,B,C
C      PARAMETER(A = 4.0D0)
C      PARAMETER(B = 0.75D0)
C      PARAMETER(C = 0.00058417226323428314253D0)
C
      EMAX = 0.5D0*(ABSP1 + ABSP2 + ABSP3)
      X1 = ABSP1/EMAX
      X2 = ABSP2/EMAX
      X3 = ABSP3/EMAX
C
      IF (X1.LT.X2) THEN
        IF (X1.LT.X3) THEN
C         X1 is smallest: X1<X2<X3 or X1<X3<X2
          FACTOR = ((1.0D0-X2)*(1.0D0-X3))**B
        ELSE
C         X3 is smallest: X3<X1<X2
          FACTOR = ((1.0D0-X1)*(1.0D0-X2))**B
        ENDIF
      ELSE
        IF (X2.LT.X3) THEN
C         X2 is smallest: X2<X1<X3 or X2<X3<X1
          FACTOR = ((1.0D0-X3)*(1.0D0-X1))**B
        ELSE
C         X3 is smallest: X3<X2<X1
          FACTOR = ((1.0D0-X1)*(1.0D0-X2))**B
        ENDIF
      ENDIF
      IF(FACTOR.LT.1.0D-15) THEN
        WRITE(NOUT,*)'FACTOR too small in RHO3',X1,X2,X3
        STOP
      ENDIF
C
      E03 = E3PAR**3
      EMAX3 = EMAX**3
      DENOM = E03 * EMAX3 * (EMAX3/E03 + 1.0D0)**2 
      DENOM = DENOM * X1 * X2 * X3 * FACTOR
      RHO3 = C/DENOM
C
      RETURN
      END
C
C23456789012345678901234567890123456789012345678901234567890123456789012
C23456789012345678901234567890123456789012345678901234567890123456789012
C
      SUBROUTINE CHOOSE2TO2S(PA,PB,ELL,OK)
C In:
      REAL*8 PA(3),PB(3)
C Out:
      REAL*8 ELL(3)
      LOGICAL OK
C
C Generates a point ell(mu) using an elliptical coordinate system
C based on the vectors p_A(mu) and p_B(mu). The points are concentrated
C near Ell(mu) = 0 and near the ellipse |ell| + |q - el| = |p_A| + |p_B|
C where q = p_A + p_B.
C 15 December 2000
C 14 Maarch 2001
C
      REAL*8 PI
      PARAMETER (PI = 3.141592653589793239D0)
C
C Input and output units.
C  
      INTEGER NIN,NOUT
      COMMON /IOUNIT/ NIN,NOUT
C
C Function RANDOM(1) give a random number in the range 0<x<1.
C
      REAL*8 RANDOM,X
C
C Function EXPM1(z) gives exp(z) - 1 while SQRTM1(z) gives sqrt(1+z) -1.
C
      REAL*8 EXPM1,SQRTM1
C
      REAL*8 ABSPA,ABSPB
      REAL*8 SUMAB(3),CROSS(3)
      REAL*8 TWOKAPPA,KAPPA,ABSCROSS
      REAL*8 NZ(3),NY(3),NX(3)
      REAL*8 PHI,AMINUS,APLUS
      REAL*8 SPLUS,TAU,CM1,C,N
      REAL*8 LAMBDA,LOGA,LOGB,SPLUSL,SMINUSL,NORM,CPLUS,TEMP
      REAL*8 LZ,LT,LX,LY
C
C-----------------------------------------------------------------------
C
      OK = .TRUE.
C
C First we calculate the unit vectors n_x, n_y, n_z used to define
C the orientation of the elliptical coordinate system. For later
C use, the variable |p_a + p_b| gets a special name, 2 kappa.
C
      ABSPA = SQRT(PA(1)**2 + PA(2)**2 + PA(3)**2)
      ABSPB = SQRT(PB(1)**2 + PB(2)**2 + PB(3)**2)
C
      SUMAB(1) = PA(1) + PB(1)
      SUMAB(2) = PA(2) + PB(2)
      SUMAB(3) = PA(3) + PB(3)
      TWOKAPPA = SQRT(SUMAB(1)**2 + SUMAB(2)**2 + SUMAB(3)**2)
      KAPPA = 0.5D0*TWOKAPPA
      CROSS(1) = PB(2)*PA(3) - PB(3)*PA(2)
      CROSS(2) = PB(3)*PA(1) - PB(1)*PA(3)
      CROSS(3) = PB(1)*PA(2) - PB(2)*PA(1) 
      ABSCROSS = SQRT(CROSS(1)**2 + CROSS(2)**2 + CROSS(3)**2)
      IF (TWOKAPPA**2 .LT. 1D-16 * ABSPA * ABSPB) THEN
        WRITE(NOUT,*) 'TWOKAPPA too small in CHOOSE2TO2S',KAPPA
        STOP
      ENDIF
      IF (ABSCROSS .LT. 1D-8 * ABSPA * ABSPB) THEN
        WRITE(NOUT,*) 'ABSCROSS too small in CHOOSE2TO2S',ABSCROSS
        WRITE(NOUT,*) 'PA is    ',PA
        WRITE(NOUT,*) 'PB is    ',PB
        WRITE(NOUT,*) 'CROSS is ',CROSS
        OK = .FALSE.
      ENDIF
      NZ(1) = SUMAB(1)/TWOKAPPA
      NZ(2) = SUMAB(2)/TWOKAPPA
      NZ(3) = SUMAB(3)/TWOKAPPA
      NY(1) = CROSS(1)/ABSCROSS
      NY(2) = CROSS(2)/ABSCROSS
      NY(3) = CROSS(3)/ABSCROSS
      NX(1) = NY(2)*NZ(3) - NY(3)*NZ(2)
      NX(2) = NY(3)*NZ(1) - NY(1)*NZ(3)
      NX(3) = NY(1)*NZ(2) - NY(2)*NZ(1)
C
C Choose phi.
C
      X = RANDOM(1)
      PHI = PI * (2.0D0*X - 1.0D0)
C
C Choose A-.
C Here N is N-/C.
C
      SPLUS  = (ABSPA + ABSPB)/TWOKAPPA
      TAU = (SPLUS - 1.0D0)/SPLUS
      CM1 = SQRTM1(TAU)
      C = CM1 + 1.0D0
      N = 1.0D0/LOG((C + 1.0D0)/CM1)
      X = RANDOM(1)
      TEMP = EXPM1((2.0D0*X - 1.0D0)/N)
      AMINUS = C * TEMP/(TEMP + 2.0D0)
C
C Choose A+.
C
      LAMBDA = (SPLUS - 1.0D0)**2/SPLUS
      LOGA = LOG(SPLUS*(SPLUS - 1.0D0 + LAMBDA)/LAMBDA)
      LOGB = LOG(SPLUS/LAMBDA)
      SPLUSL = SPLUS + LAMBDA
      SMINUSL = SPLUS - LAMBDA
C
      NORM  = 1.0D0/( LOGA/SPLUSL + LOGB/SMINUSL )
      CPLUS = 1.0D0/(1.0D0 + SPLUSL*LOGB/(SMINUSL*LOGA) )
C
      X = RANDOM(1)
      IF (X .GT. CPLUS) THEN
        TEMP = 1.0D0 - LAMBDA/SPLUS * EXP( SMINUSL*(X - CPLUS)/NORM )
        IF (TEMP.LT.1.0D-15) THEN
          WRITE(NOUT,*)'There could be a problem in CHOOSE2TO2S'
          OK = .FALSE.
        ENDIF
        APLUS = SMINUSL/TEMP
      ELSE
        TEMP = 1.0D0 + LAMBDA/SPLUS * EXP( SPLUSL*(CPLUS - X)/NORM )
        APLUS = SPLUSL/TEMP
      ENDIF
C
C We now have A+, A-, and phi so we find ell(mu).
C
      LZ = KAPPA*(1.0D0 + APLUS*AMINUS)
      LT = KAPPA*SQRT((APLUS**2 - 1.0D0)*(1.0D0 - AMINUS**2))
      LX = LT*COS(PHI)
      LY = LT*SIN(PHI)
      ELL(1) = LX*NX(1) + LY*NY(1) + LZ*NZ(1)
      ELL(2) = LX*NX(2) + LY*NY(2) + LZ*NZ(2)
      ELL(3) = LX*NX(3) + LY*NY(3) + LZ*NZ(3)
C
      RETURN
      END
C
C23456789012345678901234567890123456789012345678901234567890123456789012
C23456789012345678901234567890123456789012345678901234567890123456789012
C
      REAL*8 FUNCTION RHO2TO2S(PA,PB,ELL)
      REAL*8 PA(3),PB(3),ELL(3)
C
C Density function for points ell chosen by CHOOSE2TO2S(p_a,p_b,ell,ok).
C 15 December 2000
C 20 March 2001
C
C Input and output units.
C      
      INTEGER NIN,NOUT
      COMMON /IOUNIT/ NIN,NOUT
C
      REAL*8 PI
      PARAMETER (PI = 3.141592653589793239D0)
C
C Function SQRTM1(z) gives sqrt(1+z) -1.
C
      REAL*8 SQRTM1
C
      INTEGER MU
      REAL*8 SUMABSQ,ELLSQ,ELLPRIMESQ,PASQ,PBSQ
      REAL*8 TWOKAPPA,KAPPA,ABSELL,ABSELLPRIME,ABSPA,ABSPB
      REAL*8 SPLUS,APLUS,AMINUS
      REAL*8 TAU,CM1,C
      REAL*8 LAMBDA,LOGA,LOGB,SPLUSL,SMINUSL
      REAL*8 TEMP,DENOM
C
C-----------------------------------------------------------------------
C We start with the absolute values of combinations of vectors.
C
      SUMABSQ = 0.0D0
      ELLSQ = 0.0D0
      ELLPRIMESQ = 0.0D0
      PASQ = 0.0D0
      PBSQ = 0.0D0
      DO MU = 1,3
        TEMP = PA(MU) + PB(MU)
        SUMABSQ = SUMABSQ + TEMP**2
        ELLSQ =  ELLSQ + ELL(MU)**2
        ELLPRIMESQ =  ELLPRIMESQ  + (ELL(MU) - TEMP)**2
        PASQ = PASQ + PA(MU)**2
        PBSQ = PBSQ + PB(MU)**2
      ENDDO
      TWOKAPPA = SQRT(SUMABSQ)
      KAPPA = 0.5D0*TWOKAPPA
      ABSELL = SQRT(ELLSQ)
      ABSELLPRIME = SQRT(ELLPRIMESQ)
      ABSPA = SQRT(PASQ)
      ABSPB = SQRT(PBSQ)
C
C Now S+ and S- and A+ and A-.
C
      SPLUS  = (ABSPA + ABSPB)/TWOKAPPA
      APLUS =  (ABSELL + ABSELLPRIME)/TWOKAPPA
      AMINUS = (ABSELL - ABSELLPRIME)/TWOKAPPA
C
C Finally some auxilliary parameters.
C
      TAU = (SPLUS - 1.0D0)/SPLUS
      CM1 = SQRTM1(TAU)
      C = CM1 + 1.0D0
      LAMBDA = (SPLUS - 1.0D0)**2/SPLUS
      LOGA = LOG(SPLUS*(SPLUS - 1.0D0 + LAMBDA)/LAMBDA)
      LOGB = LOG(SPLUS/LAMBDA)
      SPLUSL = SPLUS + LAMBDA
      SMINUSL = SPLUS - LAMBDA
      DENOM = APLUS**2 - AMINUS**2
      IF (DENOM.LT.1.0D0-12) THEN
        WRITE(NOUT,*)'DENOM too small in RHO2TO2S',APLUS,AMINUS
        STOP
      ENDIF
C
C RHO is the inverse of the product of several factors.
C Phi factor:
      TEMP = 2.0D0*PI
C A- normalization:
      TEMP = TEMP*LOG((C + 1.0D0)/CM1)/C
C A- factor:
      TEMP = TEMP*(C**2 - AMINUS**2)
C A+ normalization:
      TEMP = TEMP*(LOGA/SPLUSL + LOGB/SMINUSL)
C A+ factor:
      TEMP = TEMP*APLUS*( ABS(APLUS - SPLUS) + LAMBDA )
C Jacobian for k to A+,A-,phi:
      TEMP = TEMP*KAPPA**3*DENOM
C Invert this:
      RHO2TO2S = 1.0D0/TEMP
C
      RETURN
      END
C
C23456789012345678901234567890123456789012345678901234567890123456789012
C23456789012345678901234567890123456789012345678901234567890123456789012
C
      SUBROUTINE CHOOSE2TO2T(PA,PB,ELL,OK)
C In:
      REAL*8 PA(3),PB(3)
C Out:
      REAL*8 ELL(3)
      LOGICAL OK
C
C Generates a point ell(mu) using an elliptical coordinate system
C based on the vectors p_A(mu) and p_B(mu). The points are concentrated
C near ell(mu) = 0 and near the ellipse |p_A + ell| + |p_B - ell| =
C |p_A| + |p_B|.
C 28 December 1999
C 13 March 2001
C
C The parameters A2T02T, B2TO2T, C2TO2T 
C must match between CHOOSE2TO2T and RHO2TO2T. 
C
      REAL*8 A2TO2T,B2TO2T,C2TO2T
      PARAMETER (A2TO2T = 0.3D0)
      PARAMETER (B2TO2T = 0.3D0)
      PARAMETER (C2TO2T = 6.0D0)
C
C The parameter CONST is a derived number, equal to 
C (1/g)*(log(1/g)**2 - 2*log(1/g) + 2) where g = C2TO2T. For
C g = 6, this constant is 1.4656534890040852295.
C
      REAL*8 CONST
      PARAMETER (CONST = 1.4656534890040852295D0)
C
      REAL*8 PI
      PARAMETER (PI = 3.141592653589793239D0)
C
C Input and output units.
C
      INTEGER NIN,NOUT
      COMMON /IOUNIT/ NIN,NOUT
C
C Function SQRTM1(x) gives sqrt(1+x) - 1.
C Function EXPM1(x)  gives exp(x) - 1.
C Function INVLOGSQINT(w) = y <==> w = y*(Log(y)**2 - 2*Log(y) + 2).
C Function RANDOM(1) give a random number in the range 0<x<1.
C
      REAL*8 SQRTM1,EXPM1,INVLOGSQINT,RANDOM
C
      REAL*8 SUMAB(3),CROSS(3)
      REAL*8 TWOKAPPA,ABSCROSS
      REAL*8 NZ(3),NY(3),NX(3)
      REAL*8 KAPPA,ABSPA,ABSPB,SPLUS,SMINUS
      REAL*8 X,TWOXM1,W,PHI,PHIPI
      REAL*8 LPHI2,LP2,LM2,NMINUS,CMINUS,ROOT,DAMINUS,DX,AMINUS
      REAL*8 OMEGA,SPLUSM1,BSW,A,B,NPLUS,CPLUS,FACTOR,APLUS
      REAL*8 DAPLUS,TEMP
      REAL*8 UPLUS,UMINUS,VPLUS,VMINUS,LT0,Z0,ZETA,DLT,DZ
      REAL*8 SINPHI,COSPHI,COSPHIM1,LX,LY
C
      REAL*8 TAU,CM1,C,N
      REAL*8 LAMBDA,LOGA,LOGB,SPLUSL,SMINUSL,NORM
      REAL*8 LZ,LT
C
C
C-----------------------------------------------------------------------
C
      OK = .TRUE.
C
C First we calculate the unit vectors n_x, n_y, n_z used to define
C the orientation of the elliptical coordinate system. For later
C use, the variable |p_a + p_b| gets a special name, 2 kappa.
C
      ABSPA = SQRT(PA(1)**2 + PA(2)**2 + PA(3)**2)
      ABSPB = SQRT(PB(1)**2 + PB(2)**2 + PB(3)**2)
C
      SUMAB(1) = PA(1) + PB(1)
      SUMAB(2) = PA(2) + PB(2)
      SUMAB(3) = PA(3) + PB(3)
      TWOKAPPA = SQRT(SUMAB(1)**2 + SUMAB(2)**2 + SUMAB(3)**2)
      KAPPA = 0.5D0*TWOKAPPA
      CROSS(1) = PB(2)*PA(3) - PB(3)*PA(2)
      CROSS(2) = PB(3)*PA(1) - PB(1)*PA(3)
      CROSS(3) = PB(1)*PA(2) - PB(2)*PA(1) 
      ABSCROSS = SQRT(CROSS(1)**2 + CROSS(2)**2 + CROSS(3)**2)
      IF (TWOKAPPA**2 .LT. 1D-16 * ABSPA * ABSPB) THEN
        WRITE(NOUT,*) 'TWOKAPPA too small in CHOOSE2TO2T',KAPPA
        OK = .FALSE.
      ENDIF
      IF (ABSCROSS .LT. 1D-8 * ABSPA * ABSPB) THEN
        WRITE(NOUT,*) 'ABSCROSS too small in CHOOSE2TO2T',ABSCROSS
        WRITE(NOUT,*) 'PA is    ',PA
        WRITE(NOUT,*) 'PB is    ',PB
        WRITE(NOUT,*) 'CROSS is ',CROSS
        OK = .FALSE.
      ENDIF
      NZ(1) = SUMAB(1)/TWOKAPPA
      NZ(2) = SUMAB(2)/TWOKAPPA
      NZ(3) = SUMAB(3)/TWOKAPPA
      NY(1) = CROSS(1)/ABSCROSS
      NY(2) = CROSS(2)/ABSCROSS
      NY(3) = CROSS(3)/ABSCROSS
      NX(1) = NY(2)*NZ(3) - NY(3)*NZ(2)
      NX(2) = NY(3)*NZ(1) - NY(1)*NZ(3)
      NX(3) = NY(1)*NZ(2) - NY(2)*NZ(1)
C
C At this point, we have two options. There is a main way to choose
C our point and there is a subsidiary way, which is the same as in
C CHOOSE2TO2S. If a random variable X is greater than A2TO2T, 
C we choose the main way, otherwise we choose the subsidiary way.
C
      IF (RANDOM(1).GT.A2TO2T) THEN
C
C Now some variables that depend on p_a and p_b, namely S+ and S-.
C
      SPLUS  = (ABSPA + ABSPB)/TWOKAPPA
      SMINUS = (ABSPA - ABSPB)/TWOKAPPA
C
C Choose phi.
C The standard function SIGN(xx,X) gives xx with the sign of X.
C For later use, we define PHIPI = |Phi|/Pi.
C
      X = RANDOM(1)
      TWOXM1 = 2.0D0*X - 1.0D0
      W = CONST*ABS(TWOXM1)
      PHIPI = C2TO2T*INVLOGSQINT(W)
      PHI = PI * SIGN(PHIPI,TWOXM1)
C
C Choose A-.
C There are three versions for the same thing, chosen according to
C which one should give the most accurate result, especially if
C |phi| is small.
C
      X = RANDOM(1)
      LPHI2 = LOG(C2TO2T/PHIPI)**2
      LP2   = LOG(C2TO2T/(1.0D0 + SMINUS + PHIPI))**2
      LM2   = LOG(C2TO2T/(1.0D0 - SMINUS + PHIPI))**2
      NMINUS = 1.0D0/(LPHI2 - 0.5D0 *(LP2 + LM2))
      CMINUS = 0.5D0 * NMINUS * (LPHI2 - LP2)
      IF (X .LT. 0.5D0*CMINUS) THEN
        ROOT = SQRT( 2.0D0*X/NMINUS + LP2 )
        DAMINUS =  - C2TO2T * EXP(- ROOT) + PHIPI
      ELSE IF ((1.0D0 - X) .LT. 0.5D0*(1.0D0 - CMINUS)) THEN
        ROOT = SQRT( 2.0D0*(1.0D0 - X)/NMINUS + LM2 )
        DAMINUS =  C2TO2T * EXP(- ROOT) - PHIPI
      ELSE
        DX = X - CMINUS
        ROOT = SQRT( LPHI2 - 2.0D0*ABS(DX)/NMINUS )
        DAMINUS = SIGN( C2TO2T*EXP(-ROOT) - PHIPI, DX)
      ENDIF
      AMINUS = SMINUS + DAMINUS
C
C Choose A+.
C There are three versions for the same thing, chosen according to
C which one should give the most accurate result, especially if
C omega is small. 
C EXPM1(X) gives exp(x) - 1.
C
      X = RANDOM(1)
      OMEGA = ABS(DAMINUS) + PHIPI
      SPLUSM1 = SPLUS - 1.0D0
      BSW = B2TO2T*SPLUS*OMEGA
      A = SPLUSM1 + BSW*OMEGA
      B = SPLUSM1 + C2TO2T*BSW
      NPLUS = 1.0D0 /LOG( C2TO2T**2 * A /(OMEGA**2 * B) )
      CPLUS = NPLUS * LOG( C2TO2T * A /(OMEGA * B) )
      IF (X .LT. 0.5D0*CPLUS) THEN
        FACTOR = - EXPM1( - X/NPLUS)
        APLUS = 1.0D0 
     >      + A*B*FACTOR/(BSW*(C2TO2T - OMEGA) + A*FACTOR)
        DAPLUS = APLUS - SPLUS
      ELSE IF ((1.0D0 - X) .LT. 0.5D0*(1.0D0 - CPLUS)) THEN
        FACTOR = - EXPM1((X - 1.0D0)/NPLUS)
        DAPLUS = BSW
     >     *(C2TO2T - OMEGA  - C2TO2T*FACTOR)/FACTOR
        APLUS = SPLUS + DAPLUS
      ELSE
        DX = X - CPLUS
        FACTOR = EXPM1( ABS(DX)/NPLUS )
        TEMP = BSW*OMEGA*FACTOR
        TEMP = TEMP/(1.0D0 - OMEGA/C2TO2T*(FACTOR + 1.0D0))
        DAPLUS = SIGN(TEMP,DX)
        APLUS = SPLUS + DAPLUS
      ENDIF
C
C We now have A+, A-, and phi so we find ell(mu).
C Use SQRTM1(zeta) = SQRT(1+zeta) - 1.
C
        UPLUS  = SPLUS**2 -1.0D0
        UMINUS = 1.0D0 - SMINUS**2
        VPLUS  = DAPLUS *(2.0D0*SPLUS  + DAPLUS) /UPLUS
        VMINUS = DAMINUS*(2.0D0*SMINUS + DAMINUS)/UMINUS
        LT0 = KAPPA * SQRT( UPLUS*UMINUS )
        Z0 = KAPPA * SPLUS * SMINUS
        ZETA = VPLUS - VMINUS - VPLUS*VMINUS
        DLT = LT0*SQRTM1(ZETA)
        DZ = KAPPA*(DAPLUS*SMINUS + DAMINUS*SPLUS + DAPLUS*DAMINUS)
        SINPHI = SIN(PHI)
        COSPHI = COS(PHI)
        IF (COSPHI .GT. 0.9 ) THEN
          COSPHIM1 = SQRTM1(-SINPHI**2)
        ELSE
          COSPHIM1 = COSPHI - 1.0D0
        ENDIF
        LX = DLT*COSPHI + LT0*COSPHIM1
        LY = (LT0 + DLT)*SINPHI
        ELL(1) = LX*NX(1) + LY*NY(1) + DZ*NZ(1)
        ELL(2) = LX*NX(2) + LY*NY(2) + DZ*NZ(2)
        ELL(3) = LX*NX(3) + LY*NY(3) + DZ*NZ(3)
C
      RETURN
C
C Recall, we had two options. There was a main way to choose
C our point and there was a subsidiary way, which is the same as in
C CHOOSE2TO2S. If a random variable X was greater than A2TO2T, 
C we choose the main way, otherwise we get to here and choose the 
C subsidiary way, from CHOOSE2TO2S.
C
      ELSE
C-----
C Choose phi.
C
      X = RANDOM(1)
      PHI = PI * (2.0D0*X - 1.0D0)
C
C Choose A-.
C Here N is N-/C.
C
      SPLUS  = (ABSPA + ABSPB)/TWOKAPPA
      TAU = (SPLUS - 1.0D0)/SPLUS
      CM1 = SQRTM1(TAU)
      C = CM1 + 1.0D0
      N = 1.0D0/LOG((C + 1.0D0)/CM1)
      X = RANDOM(1)
      TEMP = EXPM1((2.0D0*X - 1.0D0)/N)
      AMINUS = C * TEMP/(TEMP + 2.0D0)
C
C Choose A+.
C
      LAMBDA = (SPLUS - 1.0D0)**2/SPLUS
      LOGA = LOG(SPLUS*(SPLUS - 1.0D0 + LAMBDA)/LAMBDA)
      LOGB = LOG(SPLUS/LAMBDA)
      SPLUSL = SPLUS + LAMBDA
      SMINUSL = SPLUS - LAMBDA
C
      NORM  = 1.0D0/( LOGA/SPLUSL + LOGB/SMINUSL )
      CPLUS = 1.0D0/(1.0D0 + SPLUSL*LOGB/(SMINUSL*LOGA) )
C
      X = RANDOM(1)
      IF (X .GT. CPLUS) THEN
        TEMP = 1.0D0 - LAMBDA/SPLUS * EXP( SMINUSL*(X - CPLUS)/NORM )
        IF (TEMP.LT.1.0D-15) THEN
          WRITE(NOUT,*)'There could be a problem in CHOOSE2TO2T'
          OK = .FALSE.
        ENDIF
        APLUS = SMINUSL/TEMP
      ELSE
        TEMP = 1.0D0 + LAMBDA/SPLUS * EXP( SPLUSL*(CPLUS - X)/NORM )
        APLUS = SPLUSL/TEMP
      ENDIF
C
C We now have A+, A-, and phi so we find ell(mu).
C
      LZ = KAPPA*(1.0D0 + APLUS*AMINUS)
      LT = KAPPA*SQRT((APLUS**2 - 1.0D0)*(1.0D0 - AMINUS**2))
      LX = LT*COS(PHI)
      LY = LT*SIN(PHI)
C-----
C Here we copy the construction of ell from CHOOSE2TO2S except
C that ell_T = ell_S - PA, so we have to subtract PA.
C
      ELL(1) = LX*NX(1) + LY*NY(1) + LZ*NZ(1) - PA(1)
      ELL(2) = LX*NX(2) + LY*NY(2) + LZ*NZ(2) - PA(2)
      ELL(3) = LX*NX(3) + LY*NY(3) + LZ*NZ(3) - PA(3)
C
      RETURN
C
C End of IF (RANDOM(1).GT.A2TO2T) THEN ... ELSE ...
C
      ENDIF

      END
C
C23456789012345678901234567890123456789012345678901234567890123456789012
C23456789012345678901234567890123456789012345678901234567890123456789012
C
      REAL*8 FUNCTION RHO2TO2T(PA,PB,ELL)
      REAL*8 PA(3),PB(3),ELL(3)
C
C Density function for points ell chosen by CHOOSE2TO2T(p_a,p_b,ell,ok).
C 28 December 1999
C 13 March 2001
C
C Input and output units.      
      INTEGER NIN,NOUT
      COMMON /IOUNIT/ NIN,NOUT
C
C The parameters A2TO2T, B2TO2T, C2TO2T
C must match between CHOOSE2TO2T and RHO2TO2T. 
C
      REAL*8 A2TO2T,B2TO2T,C2TO2T
      PARAMETER (A2TO2T = 0.3D0)
      PARAMETER (B2TO2T = 0.3D0)
      PARAMETER (C2TO2T = 6.0D0)
C
C The parameter CNST is a derived number, equal to 
C 1/[2 Pi*(log(1/g)**2 - 2*log(1/g) + 2)] where g = C2TO2T. For
C g = 6, this constant is 0.0180982913407954142662161898.
C
      REAL*8 CNST
      PARAMETER (CNST = 0.0180982913407954142662D0)
C
      REAL*8 PI
      PARAMETER (PI = 3.141592653589793239D0)
C
C Function SQRTM1(x) gives sqrt(1+x) - 1.
C
      REAL*8 SQRTM1
C
      REAL*8 SUMAB(3),CROSS(3)
      REAL*8 TWOKAPPA,ABSCROSS
      REAL*8 NZ(3),NY(3),NX(3)
      REAL*8 KAPPA,PASQ,PBSQ,ABSPA,ABSPB,SPLUS,SMINUS
      REAL*8 ELLSQ,DOTAL,DOTBL
      REAL*8 WA,WB,VA,VB,DAPLUS,DAMINUS,APLUS,AMINUS
      REAL*8 DOTLSTARNX,DOTLSTARNY,PHI,PHIPI
      REAL*8 RHO0
      REAL*8 RHOPHI
      REAL*8 TEMP,NMINUS,OMEGA
      REAL*8 RHOMINUS
      REAL*8 SPLUSM1,BSW,A,B,NPLUS,DENOM1,DENOM2
      REAL*8 RHOPLUS
      REAL*8 DENOM
      REAL*8 RHOMAIN,RHOEXTRA
      REAL*8 TAU,CM1,C,LAMBDA,LOGA,LOGB,SPLUSL,SMINUSL
C
C-----------------------------------------------------------------------
C First we calculate the unit vectors n_x, n_y, n_z used to define
C the orientation of the elliptical coordinate system. For later
C use, the variable |p_a + p_b| gets a special name, 2 kappa.
C
      PASQ  = PA(1)**2  + PA(2)**2  + PA(3)**2
      PBSQ  = PB(1)**2  + PB(2)**2  + PB(3)**2
      ABSPA = SQRT(PASQ)
      ABSPB = SQRT(PBSQ)
C
      SUMAB(1) = PA(1) + PB(1)
      SUMAB(2) = PA(2) + PB(2)
      SUMAB(3) = PA(3) + PB(3)
      TWOKAPPA = SQRT(SUMAB(1)**2 + SUMAB(2)**2 + SUMAB(3)**2)
      KAPPA = 0.5D0*TWOKAPPA
      CROSS(1) = PB(2)*PA(3) - PB(3)*PA(2)
      CROSS(2) = PB(3)*PA(1) - PB(1)*PA(3)
      CROSS(3) = PB(1)*PA(2) - PB(2)*PA(1) 
      ABSCROSS = SQRT(CROSS(1)**2 + CROSS(2)**2 + CROSS(3)**2)
      IF (TWOKAPPA**2 .LT. 1D-16 * ABSPA * ABSPB) THEN
        WRITE(NOUT,*) 'TWOKAPPA too small in RHOELLIPSE',KAPPA
        STOP
      ENDIF
      IF (ABSCROSS .LT. 1D-9 * ABSPA * ABSPB) THEN
        WRITE(NOUT,*) 'ABSCROSS too small in RHO2TO2T',ABSCROSS
        WRITE(NOUT,*) 'PA is    ',PA
        WRITE(NOUT,*) 'PB is    ',PB
        WRITE(NOUT,*) 'CROSS is ',CROSS
        STOP
      ENDIF
      NZ(1) = SUMAB(1)/TWOKAPPA
      NZ(2) = SUMAB(2)/TWOKAPPA
      NZ(3) = SUMAB(3)/TWOKAPPA
      NY(1) = CROSS(1)/ABSCROSS
      NY(2) = CROSS(2)/ABSCROSS
      NY(3) = CROSS(3)/ABSCROSS
      NX(1) = NY(2)*NZ(3) - NY(3)*NZ(2)
      NX(2) = NY(3)*NZ(1) - NY(1)*NZ(3)
      NX(3) = NY(1)*NZ(2) - NY(2)*NZ(1)
C
C Now some further variables that do not depend on ell, namely
C S+ and S-.
C
      SPLUS  = (ABSPA + ABSPB)/TWOKAPPA
      SMINUS = (ABSPA - ABSPB)/TWOKAPPA
C
C Next, the dot products of ell.
C
      ELLSQ = ELL(1)**2 + ELL(2)**2 + ELL(3)**2
      DOTAL = PA(1)*ELL(1) + PA(2)*ELL(2) + PA(3)*ELL(3)
      DOTBL = PB(1)*ELL(1) + PB(2)*ELL(2) + PB(3)*ELL(3)
C
C With these we can calculate A+ and A-. We first calculate
C DA+ = A+ - S+ and DA- = A- - S- since these  variables appear in
C the density functions and they are small when ell
C is small. Thus we want to know these separately from A+ and A-.
C We use the function SQRTM1(x) = sqrt(1+x) - 1 to define 
C V_a = |p_a + ell| - |p_a| and V_b = |p_b - ell| - |p_b|.
C
      WA = (  2.0D0*DOTAL + ELLSQ)/PASQ
      WB = ( -2.0D0*DOTBL + ELLSQ)/PBSQ
      VA = ABSPA*SQRTM1(WA)
      VB = ABSPB*SQRTM1(WB)
      DAPLUS  = (VA + VB)/TWOKAPPA
      DAMINUS = (VA - VB)/TWOKAPPA
      APLUS  = DAPLUS  + SPLUS
      AMINUS = DAMINUS + SMINUS
C
C We can also calculate phi. For this, we need the the dot products of 
C L* with the unit vectors n_x and n_y. Here L* = Pa + Ell. Note that
C NY is orthogonal to Pa so for an accurate calculation in the case
C that Ell is small, we drop Pa from L* when dotting into Ny. For later
C use, we define PHIPI = |Phi|/Pi.

      DOTLSTARNX =  (PA(1) + ELL(1))*NX(1) + (PA(2) + ELL(2))*NX(2)
     >            + (PA(3) + ELL(3))*NX(3)
      DOTLSTARNY = ELL(1)*NY(1) + ELL(2)*NY(2) + ELL(3)*NY(3)
      PHI = ATAN2(DOTLSTARNY,DOTLSTARNX)
      PHIPI = ABS(PHI/PI)
C
C Now we are ready to calculate the density. First the factor rho0
C that gives the jacobian for the change of variables from 
C {ell(1), ell(2), ell(3)} to {A+,A-,phi}.
C
      DENOM = APLUS**2 - AMINUS**2
      IF (DENOM.LT.1.0D0-12) THEN
        WRITE(NOUT,*)'DENOM too small in RHO2TO2T',APLUS,AMINUS
        STOP
      ENDIF
      RHO0 = 1.0D0/(KAPPA**3 * DENOM)
C
C Next the factor for our choice of phi.
C
      RHOPHI = CNST * LOG(PHIPI/C2TO2T)**2
C
C Next the factor for our choice of A-.
C
      TEMP = LOG(C2TO2T/PHIPI)**2 
     >      -0.5D0*( LOG(C2TO2T/(1.0D0 + SMINUS + PHIPI))**2
     >             + LOG(C2TO2T/(1.0D0 - SMINUS + PHIPI))**2 )
      NMINUS = 1.0D0/TEMP
      OMEGA = ABS(DAMINUS) + PHIPI
      RHOMINUS = NMINUS * LOG(C2TO2T/OMEGA)/OMEGA
C
C Finally the factor for our choice of A+.
C
      SPLUSM1 = SPLUS - 1.0D0
      BSW = B2TO2T*SPLUS*OMEGA
      A = SPLUSM1 + BSW*OMEGA
      B = SPLUSM1 + C2TO2T*BSW
      NPLUS = 1.0D0/LOG( C2TO2T**2 * A /(OMEGA**2 * B) )
      TEMP = ABS(DAPLUS)
      DENOM1 = TEMP + BSW*OMEGA
      DENOM2 = TEMP + C2TO2T*BSW
      RHOPLUS = NPLUS*BSW*(C2TO2T - OMEGA) /(DENOM1*DENOM2)
C
C The net density is the product of the factors just calculated.
C
      RHOMAIN = RHO0*RHOPLUS*RHOMINUS*RHOPHI
C
C Now we calculate a subsidiary rho just as in RHO2TO2S.
C---------
C
      TAU = (SPLUS - 1.0D0)/SPLUS
      CM1 = SQRTM1(TAU)
      C = CM1