PROGRAM MAKERNCTEQ4

* This is an interface to the smoothened CTEQ4 parton distribution
* routine that produces a table in the standard form.
* D. E. Soper and P. Anandam.  22 October 1997.


      DOUBLE PRECISION LAMBDA
      DOUBLE PRECISION X,LX,LXMIN,LXMAX,DELTALX
      DOUBLE PRECISION MU,LMU,LMUMIN,LMUMAX,DELTALMU
      INTEGER ISET
      EXTERNAL Ctq4Pdf, SetCtq4
      DOUBLE PRECISION Ctq4Pdf

      CHARACTER*100 HEADER
      DOUBLE PRECISION LAMBDA(13)
      INTEGER LAMBDAFLAVORS(13)
      DOUBLE PRECISION XMIN,XMAX
      INTEGER NX
      DOUBLE PRECISION MUMIN,MUMAX
      INTEGER NMU
      CHARACTER*20 FILENAME(13)
      CHARACTER*20 NAME(13)
      DOUBLE PRECISION PTN(-2:6,64,64)
      DOUBLE PRECISION A, B, AP, BP, CP, TG, ANGL, PT

* Loop variables

      INTEGER I, N,M

      DATA FILENAME / '~/dat/cteq4m.ptn', '~/dat/cteq4d.ptn',
     $     '~/dat/cteq4l.ptn', '~/dat/cteq4a1.ptn', 
     $     '~/dat/cteq4a2.ptn', '~/dat/cteq4m.ptn', 
     $     '~/dat/cteq4a4.ptn', '~/dat/cteq4a5.ptn',
     $     '~/dat/cteq4hj.ptn', '~/dat/cteq4lq.ptn',
     $     '~/dat/cteq4hq.ptn', '~/dat/cteq4f3.ptn',
     $     '~/dat/cteq4f4.ptn' /

      DATA NAME /'CTEQ4M','CTEQ4D','CTEQ4L','CTEQ4A1','CTEQ4A2',
     $     'CTEQ4M','CTEQ4A4','CTEQ4A5','CTEQ4HJ','CTEQ4LQ',
     $     'CTEQ4HQ','CTEQ4F3','CTEQ4F4'/

      DATA LAMBDA /0.202D0,0.202D0,0.181D0,0.140D0,0.169D0,0.202D0,
     $        0.239D0,0.282D0,0.206D0,0.174D0,0.202D0,0.385D0,0.282D0/

      DATA LAMBDAFLAVORS /5,5,5,5,5,5,5,5,5,5,5,3,4/

      XMIN =   0.1000000000D-04
      XMAX =   0.9500000000D0
      NX   =   32
      MUMAX =  10000D0   

      DO ISET=2,13
      IF (ISET.EQ.10) THEN 
         MUMIN = 0.7D0
         NMU = 43
      ELSE 
         MUMIN = 1.6D0
         NMU = 30
      ENDIF
      CALL SetCtq4(ISET)
      WRITE(*,*) 'Making ',FILENAME(ISET)
      OPEN(UNIT = 99, FILE = FILENAME(ISET), STATUS = 'UNKNOWN')
      HEADER = '! ' // NAME(ISET)(1:7) //
     $     ' distribution generated by DES and PA 22 October 1997'
      WRITE(99,*)HEADER
      WRITE(99,*)'!'
      WRITE(99,*)'PARTON DATA'
      WRITE(99,*)'  VERSION =        110496'
      WRITE(99,103)LAMBDA(ISET)
103   FORMAT('    LAMBDA',' = ',G20.10)
      WRITE(99,104)LAMBDAFLAVORS(ISET)
104   FORMAT('NFL_LAMBDA',' = ',I10)
      WRITE(99,105)XMIN
105   FORMAT('      XMIN',' = ',G20.10)
      WRITE(99,106)XMAX
106   FORMAT('      XMAX',' = ',G20.10)
      WRITE(99,107)NX
107   FORMAT(' N_XPOINTS',' = ',I10)
      WRITE(99,108)MUMIN
108   FORMAT('     MUMIN',' = ',G20.10)
      WRITE(99,109)MUMAX
109   FORMAT('     MUMAX',' = ',G20.10)
      WRITE(99,110)NMU
110   FORMAT('N_MUPOINTS',' = ',I10)

* Calculate lattice info.

      LXMIN  = DLOG(XMIN/(1.0D0 - XMIN))
      LXMAX  = DLOG(XMAX/(1.0D0 - XMAX))
      DELTALX = (LXMAX - LXMIN)/(NX-1)
      LMUMIN = DLOG(MUMIN)
      LMUMAX = DLOG(MUMAX)
      DELTALMU = (LMUMAX - LMUMIN)/(NMU-1)

      LX = LXMIN - DELTALX
      DO N=1,NX 
        LX = LX + DELTALX
        X = DEXP(LX)/(1.0D0+DEXP(LX))
        LMU = LMUMIN - DELTALMU   
        DO M=1,NMU
          LMU = LMU + DELTALMU
          MU = DEXP(LMU)

          PTN(-2,N,M) = Ctq4Pdf(-2, X, MU)
          PTN(-1,N,M) = Ctq4Pdf(-1, X, MU)
          PTN(0,N,M)  = Ctq4Pdf(0, X, MU)
          PTN(1,N,M)  = Ctq4Pdf(1, X, MU)
          PTN(2,N,M)  = Ctq4Pdf(2, X, MU)
          PTN(3,N,M)  = Ctq4Pdf(3, X, MU)
          PTN(4,N,M)  = Ctq4Pdf(4, X, MU)
          PTN(5,N,M)  = Ctq4Pdf(5, X, MU)
          PTN(6,N,M)  = Ctq4Pdf(6, X, MU)

* We check for negative values and set those values to zero.

	 IF (PTN(-2,N,M).LT.0D0) PTN(-2,N,M) = 0D0
	 IF (PTN(-1,N,M).LT.0D0) PTN(-1,N,M) = 0D0
	 IF (PTN(0,N,M).LT.0D0)  PTN(0,N,M)  = 0D0
	 IF (PTN(1,N,M).LT.0D0)  PTN(1,N,M)  = 0D0
	 IF (PTN(2,N,M).LT.0D0)  PTN(2,N,M)  = 0D0
	 IF (PTN(3,N,M).LT.0D0)  PTN(3,N,M)  = 0D0
	 IF (PTN(4,N,M).LT.0D0)  PTN(4,N,M)  = 0D0
	 IF (PTN(5,N,M).LT.0D0)  PTN(5,N,M)  = 0D0

         PTN(-2,N,M) = DLOG(1.0D-16 + PTN(-2,N,M))
         PTN(-1,N,M) = DLOG(1.0D-16 + PTN(-1,N,M))
         PTN(0,N,M)  = DLOG(1.0D-16 + PTN(0,N,M))
         PTN(1,N,M)  = DLOG(1.0D-16 + PTN(1,N,M))
         PTN(2,N,M)  = DLOG(1.0D-16 + PTN(2,N,M))
         PTN(3,N,M)  = DLOG(1.0D-16 + PTN(3,N,M))
         PTN(4,N,M)  = DLOG(1.0D-16 + PTN(4,N,M))
         PTN(5,N,M)  = DLOG(1.0D-16 + PTN(5,N,M))
         PTN(6,N,M)  = DLOG(1.0D-16 + 0.0D0)

       ENDDO
      ENDDO

* Smoothen out the distributions
     
      DO I = -2, 6
         DO M = 1, NMU
            DO N = 2, NX

* The distributions shouldn't fall to zero too sharply

               IF ( (PTN(I,N,M).LT.DLOG(1D-14)) . AND. 
     $              (PTN(I,N-1,M).GT.DLOG(1D-10)) )
     $              PTN(I,N,M) = (PTN(I,N-1,M)+PTN(I,N,M)) / 2D0
            ENDDO
         ENDDO
      ENDDO               

* Check for smoothness by looking at angles

      DO I = -2, 6
         DO M = 1, NMU
            DO N = 2, NX-1
               IF (PTN(I,N+1,M).GT.DLOG(1D-14)) THEN
                  A = PTN(I,N,M) - PTN(I,N-1,M)
                  B = PTN(I,N+1,M) - PTN(I,N,M)
                  ANGL = DATAN((B-A)/(1D0+B*A))/3.1415D0*180D0
                  IF (ANGL.LT.-20D0) THEN
                     TG = DTAN( DSIGN(20D0 / 180D0 * 3.14159D0,ANGL))
                     AP = - TG
                     BP = TG*PTN(I,N+1,M) + TG*PTN(I,N-1,M) + 2D0
                     CP = TG - TG*PTN(I,N+1,M)*PTN(I,N-1,M) 
     $                    - PTN(I,N+1,M) - PTN(I,N-1,M)
                     PT = (-BP + DSQRT(BP**2-4D0*AP*CP))
     $                    /(2*AP)
                     PTN(I,N,M) = PT
                  ENDIF
                  IF (ANGL.GT.15D0) THEN
                     PT = (PTN(I,N+1,M) + PTN(I,N-1,M))/2D0
                     PTN(I,N,M) = PT
                  ENDIF
               ENDIF
            ENDDO
         ENDDO
      ENDDO

* The smoothened parton distributions are written out

      DO N = 1, NX
         DO M = 1, NMU
            WRITE(99,'(9G13.5)') (PTN(I,N,M), I=-2,6)
         ENDDO
      ENDDO

      CLOSE(99)
      ENDDO

      STOP
      END