C *********************************************************************
C  aac.f    Version 1.1                                Aug. 17, 2003
C
C  [Package for AAC polarized parton distribution functions]
C *********************************************************************
C  AAC Polarized Parton Distribution Functions.
C      (Asymmetry Analysis Collaboration)
C
C  This package contains 3 sets of parton distribution functions:
C    1) AAC Leading Order (LO),
C    2) AAC Next-to-Leading Order (NLO), Set-1,
C    3) AAC Next-to-Leading Order (NLO), Set-2.
C
C  REFERENCE:
C    "Polarized Parton Distribution Functions in the Nucleon",
C     Y. Goto, N. Hayashi, M. Hirai, H. Horikawa, S. Kumano,
C     M. Miyama, T. Morii, N. Saito, T.-A. Shibata, E. Taniguchi,
C     and T. Yamanishi (Asymmetry Analysis Collaboration),
C     Phys. Rev. D62 (2000) 34017.
C
C  QUESTIONS OR COMMENTS TO:
C    Masanori Hirai     [mhirai@rarfaxp.riken.go.jp]
C ---------------------------------------------------------------------
C  SUBROUTINE AACPDF(IREAD,ISET,Q2,X,POLPDF,STRUCT):
C
C   Subroutine AACPDF returns the values of AAC polarized parton
C   distributions, structure functions g_1, and spin asymmetries A_1
C   at specified Q^2 and x point by interpolating the grid data
C   which are given in grid files [Log(Q^2): LINEAR INTERPOLATION,
C   x: CUBIC SPLINE INTERPOLATION].
C
C   INPUT:
C
C     IREAD 1 ... read a grid file
C           2 ... do not read a grid file
C
C     ISET  1 ... AAC LO (aaclo.grid)
C           2 ... AAC NLO Set-1 (aacnlo1.grid)
C           3 ... AAC NLO Set-2 (aacnlo2.grid)
C
C     Q2, X   ... Q^2 and x values at which the distributions are
C                 calculated. Available range: 10^-9 <= X <= 1.0,
C                 1.0 <= Q2 <= 100,000 GeV^2.
C
C     NOTE: 1) The grid files should be in the working directory.
C           2) IREAD should be "1" when this subroutine is called
C              for the first time or the ISET has been changed.
C
C   OUTPUT: Arrays POLPDF(-3:3) and STRUCT(3)
C
C     POLPDF(I) --> Polarized Parton Distribution Functions.
C      I = -3 ... s-bar quark (x Delta s-bar = x Delta s)
C          -2 ... d-bar quark (x Delta d-bar = x Delta d_sea)
C          -1 ... u-bar quark (x Delta u-bar = x Delta u_sea)
C           0 ... gluon (x Delta g)
C           1 ... u-valence quark (x Delta u_v)
C           2 ... d-valence quark (x Delta d_v)
C           3 ... s quark (x Delta s = x Delta s-bar)
C
C     STRUCT(I) --> Structure Functions g_1.
C      I =  1 ... x g_1 for proton
C           2 ... x g_1 for neutron
C           3 ... x g_1 for deuteron
C
C     NOTE: The returned values for the parton distribution functions
C           and the structure functions are multiplied by x.
C *********************************************************************
      SUBROUTINE AACPDF(ISET,Q2,X,POLPDF,STRUCT)
C ---------------------------------------------------------------------
      IMPLICIT REAL*8(A-H,O-Z)
      PARAMETER (NQ=23, NX=68, ND=9, NSET=3,IFILE=11)
      DIMENSION IREAD(2), POLPDF(-3:3), STRUCT(3), QG(NQ), XG(NX),
     1          PDFG(NX,NQ,ND), PDFJ1(ND), PDFJ2(ND),
     2          BXG(NX,NQ,ND), CXG(NX,NQ,ND), DXG(NX,NQ,ND)
      SAVE IREAD, BXG, CXG, DXG

      DATA IREAD / 1, 1 /

C Q2 AND X GRID.
      DATA QG / 1.0D0, 1.3D0, 1.8D0, 2.7D0, 4.0D0, 6.4D0,
     1          1.0D1, 1.6D1, 2.5D1, 4.0D1, 6.4D1,
     2          1.0D2, 1.8D2, 3.2D2, 5.7D2,
     3          1.0D3, 1.8D3, 3.2D3, 5.7D3,
     4          1.0D4, 2.2D4, 4.6D4,
     5          1.0D5 /
      DATA XG / 1.0D-9, 1.8D-9, 3.2D-9, 5.7D-9,
     1          1.0D-8, 1.8D-8, 3.2D-8, 5.7D-8,
     2          1.0D-7, 1.8D-7, 3.2D-7, 5.7D-7,
     3          1.0D-6, 1.4D-6, 2.0D-6, 3.0D-6, 4.5D-6, 6.7D-6,
     4          1.0D-5, 1.4D-5, 2.0D-5, 3.0D-5, 4.5D-5, 6.7D-5,
     5          1.0D-4, 1.4D-4, 2.0D-4, 3.0D-4, 4.5D-4, 6.7D-4,
     6          1.0D-3, 1.4D-3, 2.0D-3, 3.0D-3, 4.5D-3, 6.7D-3,
     7          1.0D-2, 1.4D-2, 2.0D-2, 3.0D-2, 4.5D-2,
     8          0.06D0, 0.08D0, 0.1D0, 0.125D0, 0.15D0, 0.175D0,
     9          0.2D0, 0.225D0, 0.25D0, 0.275D0, 0.3D0, 0.325D0,
     A          0.35D0, 0.375D0, 0.4D0,  0.45D0, 0.5D0, 0.55D0,
     B          0.6D0, 0.65D0,  0.7D0,  0.75D0,  0.8D0, 0.85D0,
     C          0.9D0, 0.95D0, 1.0D0 /

      IF((IREAD(1).NE.1).AND.(IREAD(2).EQ.ISET)) GO TO 40

C READ GRID DATA AND CALCULATE SPLINE COEFFICIENTS. 
      IF(ISET.EQ.1) THEN
        OPEN(UNIT=IFILE,FILE='aaclo.grd',STATUS='OLD')

      ELSE IF(ISET.EQ.2) THEN
        OPEN(UNIT=IFILE,FILE='aacnlo1.grd',STATUS='OLD')

      ELSE IF(ISET.EQ.3) THEN
        OPEN(UNIT=IFILE,FILE='aacnlo2.grd',STATUS='OLD')

      ELSE
        WRITE(*,1010) ISET
 1010   FORMAT(' ','AACPDF ERROR: ISET =', I3)
        STOP
      END IF

      DO 10 J=1,NQ
        DO 10 K=1,NX-1
          READ(IFILE,1020) (PDFG(K,J,I), I=1,6)
 1020       FORMAT(1X,6(1PE12.4))
   10 CONTINUE

      DO 20 J=1,NQ
        DO 20 K=1,NX-1
          READ(IFILE,1020) (PDFG(K,J,I), I=7,9)
   20 CONTINUE
      CLOSE(IFILE)

      DO 30 I=1,ND
        DO 30 J=1,NQ
          PDFG(NX,J,I)=0.D0
          CALL SSPLINE(NX,XG,PDFG,BXG,CXG,DXG,I,J)
   30 CONTINUE

      IREAD(1)=2
      IREAD(2)=ISET

   40 CONTINUE
C CHECK X AND Q2 VALUES.
      IF((X.LT.1.D-9).OR.(X.GT.1.D0)) THEN
        WRITE(*,1030) X
 1030   FORMAT (' ','AACPDF WARNING: OUT OF RANGE --> X =', 1PE12.3)
        STOP
      ENDIF
      IF((Q2.LT.1.D0).OR.(Q2.GT.1.D5)) THEN
        WRITE(*,1040) Q2
 1040   FORMAT (' ','AACPDF WARNING: OUT OF RANGE --> Q2 =', 1PE12.3)
        STOP
      ENDIF

C INTERPOLATION.
C X: CUBIC SPLINE INTERPOLATION, LOG(Q2): LINEAR INTERPOLATION.
      J=ISERCH(NQ,QG,Q2)
      IF(J.EQ.NQ) J=NQ-1
      K=ISERCH(NX,XG,X)
      DO 50 I=1,ND
        DX=X-XG(K)
        PDFJ1(I)=PDFG(K,J,I)
     1       +DX*(BXG(K,J,I)+DX*(CXG(K,J,I)+DX*DXG(K,J,I)))
        PDFJ2(I)=PDFG(K,J+1,I)
     1       +DX*(BXG(K,J+1,I)+DX*(CXG(K,J+1,I)+DX*DXG(K,J+1,I)))
   50 CONTINUE
      T=(DLOG(Q2)-DLOG(QG(J)))/(DLOG(QG(J+1))-DLOG(QG(J)))
      DO 60 I=1,3
        POLPDF(I-1)=(1.D0-T)*PDFJ1(I)+T*PDFJ2(I)
        POLPDF(-I)=(1.D0-T)*PDFJ1(I+3)+T*PDFJ2(I+3)
        STRUCT(I)=(1.D0-T)*PDFJ1(I+6)+T*PDFJ2(I+6)
   60 CONTINUE
      POLPDF(3)=POLPDF(-3)

      RETURN
      END
C ---------------------------------------------------------------------
      SUBROUTINE SSPLINE(N,X,Y,B,C,D,I,J)
C ---------------------------------------------------------------------
C CALCULATE THE COEFFICIENTS B,C,D IN A CUBIC SPLINE INTERPOLATION.
C INTERPOLATION SUBROUTINES ARE TAKEN FROM
C G.E. FORSYTHE, M.A. MALCOLM AND C.B. MOLER,
C COMPUTER METHODS FOR MATHEMATICAL COMPUTATIONS (PRENTICE-HALL, 1977).
      IMPLICIT REAL*8(A-H,O-Z)
      PARAMETER (NQ=23, NX=68, ND=12)
      DIMENSION X(NX), Y(NX,NQ,ND),
     1          B(NX,NQ,ND), C(NX,NQ,ND), D(NX,NQ,ND)
      NM1=N-1
      IF(N.LT.2) RETURN
      IF(N.LT.3) GO TO 250
      D(1,J,I)=X(2)-X(1)
      C(2,J,I)=(Y(2,J,I)-Y(1,J,I))/D(1,J,I)
      DO 210 K=2,NM1
        D(K,J,I)=X(K+1)-X(K)
        B(K,J,I)=2.0D0*(D(K-1,J,I)+D(K,J,I))
        C(K+1,J,I)=(Y(K+1,J,I)-Y(K,J,I))/D(K,J,I)
        C(K,J,I)=C(K+1,J,I)-C(K,J,I)
  210 CONTINUE
      B(1,J,I)=-D(1,J,I)
      B(N,J,I)=-D(N-1,J,I)
      C(1,J,I)=0.0D0
      C(N,J,I)=0.0D0
      IF(N.EQ.3) GO TO 215
      C(1,J,I)=C(3,J,I)/(X(4)-X(2))-C(2,J,I)/(X(3)-X(1))
      C(N,J,I)=C(N-1,J,I)/(X(N)-X(N-2))-C(N-2,J,I)/(X(N-1)-X(N-3))
      C(1,J,I)=C(1,J,I)*D(1,J,I)**2.0D0/(X(4)-X(1))
      C(N,J,I)=-C(N,J,I)*D(N-1,J,I)**2.0D0/(X(N)-X(N-3))
  215 CONTINUE
        DO 220 K=2,N
        T=D(K-1,J,I)/B(K-1,J,I)
        B(K,J,I)=B(K,J,I)-T*D(K-1,J,I)
        C(K,J,I)=C(K,J,I)-T*C(K-1,J,I)
  220 CONTINUE
      C(N,J,I)=C(N,J,I)/B(N,J,I)
      DO 230 IB=1,NM1
        K=N-IB
        C(K,J,I)=(C(K,J,I)-D(K,J,I)*C(K+1,J,I))/B(K,J,I)
  230 CONTINUE
      B(N,J,I)=(Y(N,J,I)-Y(NM1,J,I))/D(NM1,J,I)
     1        +D(NM1,J,I)*(C(NM1,J,I)+2.0D0*C(N,J,I))
      DO 240 K=1,NM1
        B(K,J,I)=(Y(K+1,J,I)-Y(K,J,I))/D(K,J,I)
     1          -D(K,J,I)*(C(K+1,J,I)+2.0D0*C(K,J,I))
        D(K,J,I)=(C(K+1,J,I)-C(K,J,I))/D(K,J,I)
        C(K,J,I)=3.0D0*C(K,J,I)
  240 CONTINUE
      C(N,J,I)=3.0D0*C(N,J,I)
      D(N,J,I)=D(N-1,J,I)
      RETURN
  250 CONTINUE
      B(1,J,I)=(Y(2,J,I)-Y(1,J,I))/(X(2)-X(1))
      C(1,J,I)=0.0D0
      D(1,J,I)=0.0D0
      B(2,J,I)=B(1,J,I)
      C(2,J,I)=0.0D0
      D(2,J,I)=0.0D0
      RETURN
      END
C ---------------------------------------------------------------------
      INTEGER FUNCTION ISERCH(N,X,Y)
C ---------------------------------------------------------------------
C THIS FUNCTION SEARCHES "I" WHICH SATISFIES THE RELATION
C X(I) <= Y < X(I+1) BY USING A BINARY SEARCH.
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION X(70)

      MIN=1
      MAX=N+1

   10 CONTINUE
      MID=(MIN+MAX)/2
      IF(Y.LT.X(MID)) THEN
        MAX=MID
      ELSE
        MIN=MID
      END IF
      IF((MAX-MIN).GT.1) GO TO 10

      ISERCH=MIN

      RETURN
      END
C *********************************************************************
C THE END OF THE PROGRAM.
C *********************************************************************
