C********************************************************************
C                                                                   *
C     TITEL:       SUMREG                                           *
C                                                                   *
C     VERSION:     11.MAI.1994                                      *
C     REVISED:                                                      *
C                                                                   *
C     SUBROUTINES:  nag-routines                                    *
C                                                                   *
C     Description:  This routine calculates the two dimensional     *
C                   integral over the spectral function             *
C                   int d3p_N dE S(p_N,E)                           *
C               In addition:                                        *
C               This routine is intended to give precise hints on   *
C               how the spectral function is normalized and saved.  *
C               The spectral function is based on the three-body    *
C               wave function calculated on a finite grid with      *
C               40*30 mesh points for the (p,q=p_N)-momenta. It is  *
C               again calculated on a finite grid with (30*30) mesh *
C               points for the momentum/energy mesh (p_N,E).        *
C                                                                   * 
C                                                                   *
C     Compile with: %f77 -r8 sumreg.f                               * 
C                   %f77 -o prog sumreg.o -lnag                     *
C                                                                   *            
C********************************************************************
      PROGRAM SUMREG
C

      PARAMETER (NDE=32,NDE1=NDE-1,NDP=40,NDQ=30)
      REAL Y(NDE1),X(NDE1),ANS
C
C
      COMMON/MESH/ENER(NDE),QMESH(NDQ),QWEIGH(NDQ),NERG,NQ
      COMMON /MASS/ XM,EMDEL,HBCQ,PI
      COMMON/SPEC/R0T0(NDP,NDQ),R0T1(NDP,NDQ),
     *        R1T0(NDP,NDQ),R1T1(NDP,NDQ),R2T0(NDP,NDQ),R2T1(NDP,NDQ)

C     CONST.
C
      PI=4.0*ATAN(1.0)
C
C
      OPEN(UNIT=1,FILE='sfparis.dat')
      OPEN(UNIT=6,FILE='specout.dat')
C
C-----input of the spectral function:
C
C     The spectral function is given for the pair-isospin 
C     T=0 and T=1
C     The two-body breakup is given in the first columne, i.e.,
C     R0T0(1,J) etc.
C     R0T0(2..NERG,J) contains the three-body breakup.
C     The normalization is chosen such that the proton and 
C     neutron contribution is given by:
C
C      R0P = (F0T0 + F0T1)*4.*PI/3.
C      R0N = 2.*F0T1*4.*PI/3.
C
C     with R0P = 2/3
C     and  R0N = 1/3
C     
C     Since the spectral function is given on a finite mesh
C     the sumrule is not exactly fullfilled !
C     The mesh was chosen for a particular kinematic situation.
C
C     The momentum mesh is given for Gaussian points.  
C     
      READ(1,'(2I3)')NERG,NQ
      READ(1,'(6E13.6)')(ENER(IE),IE=1,NERG)
      READ(1,'(6E13.6)')(QMESH(IQ),IQ=1,NQ),(QWEIGH(IQ),IQ=1,NQ)
      READ(1,'(6E13.6)')((R0T0(I,J),R1T0(I,J),R2T0(I,J),
     *     J=1,NQ),I=1,NERG)
      READ(1,'(6E13.6)')((R0T1(I,J),R1T1(I,J),R2T1(I,J),
     *     J=1,NQ),I=1,NERG)
C
      F0T0=0.0
      F1T0=0.0
      F2T0=0.0
      F0T1=0.0
      F1T1=0.0
      F2T1=0.0
      DO 100 IEN = 2,NERG
100     X(IEN-1) = ENER(IEN)
C
      NERG1 = NERG - 1
C
C     SUMMATION OVER ALL MOMENTA
      DO 10000 JQ =1,NQ 
C
        DO 1000 IEN=2,NERG
          Y(IEN-1) = R0T0(IEN,JQ)
1000    CONTINUE
C
      IFAIL = 1
      CALL D01GAF(X,Y,NERG1,ANS,ER,IFAIL)
      F0T0    = (ANS + R0T0(1,JQ))*QMESH(JQ)**2*QWEIGH(JQ) + F0T0
C
        DO 1001 IEN=2,NERG
          Y(IEN-1) = R0T1(IEN,JQ)
1001    CONTINUE
C
      IFAIL = 1
      CALL D01GAF(X,Y,NERG1,ANS,ER,IFAIL)
      F0T1    = ANS*QMESH(JQ)**2*QWEIGH(JQ) + F0T1
C
        DO 1002 IEN=2,NERG
          Y(IEN-1) =R1T0(IEN,JQ)
1002    CONTINUE
C
      IFAIL = 1
      CALL D01GAF(X,Y,NERG1,ANS,ER,IFAIL)
      F1T0     =(ANS + R1T0(1,JQ))*QMESH(JQ)**2*QWEIGH(JQ)+ F1T0
C
        DO 1003 IEN=2,NERG
          Y(IEN-1) =R1T1(IEN,JQ)
1003    CONTINUE
C
      IFAIL = 1
      CALL D01GAF(X,Y,NERG1,ANS,ER,IFAIL)
      F1T1     = ANS*QMESH(JQ)**2*QWEIGH(JQ) + F1T1
C
        DO 1004 IEN=2,NERG
          Y(IEN-1) =R2T0(IEN,JQ)
1004    CONTINUE
C
      IFAIL = 1
      CALL D01GAF(X,Y,NERG1,ANS,ER,IFAIL)
      F2T0     =(ANS + R2T0(1,JQ))*QMESH(JQ)**2*QWEIGH(JQ)+ F2T0
C
        DO 1005 IEN=2,NERG
          Y(IEN-1) =R2T1(IEN,JQ)
1005    CONTINUE
C
      IFAIL = 1
      CALL D01GAF(X,Y,NERG1,ANS,ER,IFAIL)
      F2T1     = ANS*QMESH(JQ)**2*QWEIGH(JQ) + F2T1
C
C
10000 CONTINUE
C
      R0P = (F0T0 + 1.*F0T1)*4.*PI/3.
      R0N = 2.*F0T1*4.*PI/3.
      R1P = (F1T0 + 1.*F1T1)*4.*PI/3.
      R1N = 2.*F1T1*4.*PI/3.
      R2P = (F2T0 + 1.*F2T1)*4.*PI/3.
      R2N = 2.*F2T1*4.*PI/3.
C
      WRITE(6,996)
      WRITE(6,997) R0P,R0N,R0P+R0N
      WRITE(6,998) R1P,R1N
      WRITE(6,999) R2P,R2N
C
996   FORMAT('      Proton       Neutron     Normalization')
997   FORMAT(' R0:',3E13.6)
998   FORMAT(' R1:',2E13.6)
999   FORMAT(' R2:',2E13.6)
C
      END