* V. Guzey, July 2010, JLab

* Nuclear IMPACT-PARAMETER-DEPENDENT parton distributions 
* in the leading twist theory of nuclear shadowing
* L. Frankfurt, V. Guzey and M. Strikman, Phys. Rept. (in preparation)
* Suggested abbreviation: FGS10

* For the nucleus of Ca-40

* FOR ANY GIVEN VALUE OF BJORKEN X (10^{-5} < X < 0.95),
* Q^2 (4 < Q^2 < 16,000) GEV^2 AND IMPACT PARAMETER B (0 < b < 10 fm), 
* THE CODE GIVES THE RATIOS OF NUCLEAR TO FREE
* NUCLEON NEXT-TO-LEADING ORDER (NLO) DISTRIBUTION FUNCTIONS AND ALSO 
* THE RATIO OF THE NUCLEAR TO NUCLEON
* NLO INCLUSIVE STRUCTURE FUNCTIONS F_2.
* THE WORD "NUCLEON" MEANS (Z*proton+N*neutron)/A

* FOR THE FREE PROTON NLO DISTRIBUTON FUNCTION, CTEQ5M PARAMETRIZATION IS USED.

* NUCLEAR SHADOWING FOR THE VALENCE QUARKS IS GIVEN BY THE PARAMETRIZAION
* BY ESKOLA AND COLLABORATORS
* (K.J. Eskola, V.J. Kolhinen and P.V. Ruuskanen, Nucl. Phys. B 535 (1998) 351;
* K.J. Eskola, V.J. Kolhinen and C.A. Salgado, Eur. Phys. J. C 9 (1999) 61.)

*  **** HOW TO USE THIS CODE ***

* THE MAIN SOUBROUTINE IS CALLED ''LT_IMPACT(bb,xbj,Q2,uv,dv,ubar,dbar,sbar,cbar,glue,f2)''
* WHEN THE USER SPECIFIES THE VALUES OF BJORKEN X (xbj), Q^2 (Q2), and impact param. bb (b),
* THE SUBROUTINE OUTPUTS THE FOLLOWING RATIOS AT THE GIVEN X AND Q^2:
*
* U_VAL(NUCLEUS)/[T_A(b)*(Z*U_VAL+N*D_VAL)] .... uv    (valence u-quark ratios)
* D_VAL(NUCLEUS)/[T_A(b)*(Z*D_VAL+N*U_VAL)] .... dv    (valence d-quark ratios)
* U_BAR(NUCLEUS)/[T_A(b)*(Z*U_BAR+N*D_BAR)] .... ubar  (anti u-quark ratios)
* D_BAR(NUCLEUS)/[T_A(b)*(Z*D_BAR+N*U_BAR)] .... dbar  (anti d-quark ratios)
* S_BAR(NUCLEUS)/[T_A(b)*A*S_BAR]           .... sbar  (anti s-quark ratios)
* C_BAR(NUCLEUS)/[T_A(b)A*C_BAR]            .... cbar  (anti c-quark ratios)
* GLUE(NUCLEUS)/T_A(b)*A*GLUE]              .... glue  (gluon ratios)
* F_2(NUCLEUS)/[T_A(b)*(Z*F_2p+N*F_2n)]     .... f2    (F_2 ratios)

* where T_A(b) is the nuclear optical density,  T_A(b)=\int dz \rho_A(b,z)
* T_A(b) is normalized to unity: \int d^2b T_A(b)=1
* Note that T_A(b) is calculated separametely by subroutine TA_ca40(bb,TT)
* [TT=T_A(b)] which uses date file TA_ca40.dat

* **** WHAT DOES THE CODE DO?

* THE CODE READS THE DATA FILE (QCDEvolution_XXproton_2009_IMPACT_modelX.dat) AND
* INTERPOLATES BETWEEN THE X, Q^2 and b GRID POINTS OF THE  DATA FILE 
* (subroutine polin3)  USING THE LINEAR INTERPOLATION (subroutine linear). 
* THE DATA FILE CONTAINS 90 GRID POINTS IN X, 7 GRID POINTS IN Q^2  
* (ACTUALLY, IN LOG(Q^2)/LOG(4.)), and 20 points in b.

* Note that b_min=0.0.414 fm and b_max=7.059 fm
* The results are flat in b for small b. Hence, b=0 fm is equivalent to b=b_min

      implicit double precision (a-h,o-z)
      

      b_min=0.414
      b_max=7.059
      Q2_min=4.0  

      bb=0.0
      Q2=4.0    
      xbj=0.0001

      if (bb.lt.b_min) bb=b_min  
 
      call TA_ca40(bb,TT)  
      call  LT_IMPACT(bb,xbj,Q2,uv,dv,ubar,dbar,sbar,cbar,
     >glue,f2) 
        
      write(6,7)TT,ubar,glue,cbar,f2
 7    format(9f8.5) 


      stop
      end

      subroutine TA_ca40(bb,TT)
      implicit double precision (a-h,o-z)    
      integer iread_t  
      dimension bb_array(20),TT_array(20)  
        
      data iread_t/0/
      if (iread_t.eq.0) then 

      NB=20  
      open(11,file='TA_ca40.dat',status='unknown') 
      do ib=1,NB
      read(11,*) bb_array(ib),TT_array(ib)
      enddo
      close(11)  
      endif
      iread_t=1
    
      call linear(NB,bb_array,TT_array,bb,TT)

      return
      end  


      subroutine LT_IMPACT(bb,xbj,Q2,uv,dv,ubar,dbar,sbar,cbar,
     >glue,f2)

      implicit double precision (a-h,o-z)

      integer iread  

      dimension x1a(90),x2a(7),ya1(90,7,20),ya2(90,7,20),
     >ya3(90,7,20),ya4(90,7,20),ya5(90,7,20),ya6(90,7,20)
     >,ya7(90,7,20),ya8(90,7,20),res(8)

      common/fgs03/x1a,x2a,ya1,ya2,ya3,ya4,
     >ya5,ya6,ya7,ya8

      data iread/0/
      if (iread.eq.0) then

      open(1,file=
     >'QCDEvolution_ca40proton_2009_IMPACT_model2.dat',
     >status='unknown')

      do k=1,20
      do j=1,7
      read(1,*)x2a(j)  
      do i=1,90
      read(1,*)x1a(i),ya1(i,j,k),ya2(i,j,k),
     >ya3(i,j,k),ya4(i,j,k),ya5(i,j,k)
     >,ya6(i,j,k),ya7(i,j,k),ya8(i,j,k)

      enddo
      enddo
      enddo

      close(1)

      endif
      iread=1

      call polin3(bb,xbj,q2,res)

      dv=res(1)
      uv=res(2)
      ubar=res(3)
      dbar=res(4)
      sbar=res(5)
      cbar=res(6)
      glue=res(7)
      f2=res(8)

      return
      end

      subroutine polin3(bb,x1,x2,res)
      
      implicit double precision (a-h,o-z)

      integer j,k,iparton,NMAX,MMAX,NB

      dimension x1a(90),x2a(7),
     > B(100),C(100),D(100),res(8)

      dimension yntmp1(7),ymtmp1(90),yntmp2(7),
     >ymtmp2(90),yntmp3(7),ymtmp3(90),yntmp4(7),
     >ymtmp4(90),yntmp5(7),ymtmp5(90),yntmp6(7),
     >ymtmp6(90),yntmp7(7),ymtmp7(90),yntmp8(7),
     >ymtmp8(90)

      dimension ya1(90,7,20),ya2(90,7,20),
     >ya3(90,7,20),ya4(90,7,20),ya5(90,7,20),ya6(90,7,20)
     >,ya7(90,7,20),ya8(90,7,20)

      dimension Bold(20),wBold(20),b1a(20),res1a(20),res2a(20),
     >res3a(20),res4a(20),res5a(20),res6a(20),res7a(20),res8a(20)

      parameter(NMAX=7,MMAX=90,NB=20)

      common/fgs03/x1a,x2a,ya1,ya2,ya3,ya4,
     >ya5,ya6,ya7,ya8

      x2a(1)=1.
      x2a(2)=2.
      x2a(3)=3.
      x2a(4)=4.
      x2a(5)=5.
      x2a(6)=6.
      x2a(7)=7.  

      call GAUSAB(Bold,wBold,0.,50.,NB)
      do ib=1,NB
      b1a(ib)=Bold(ib)**0.5
c      print*,b1a(ib)
      enddo

      do ib=1,20
      do j=1,90

      do k=1,7

      yntmp1(k)=ya1(j,k,ib)
      yntmp2(k)=ya2(j,k,ib)
      yntmp3(k)=ya3(j,k,ib)
      yntmp4(k)=ya4(j,k,ib)
      yntmp5(k)=ya5(j,k,ib)
      yntmp6(k)=ya6(j,k,ib)
      yntmp7(k)=ya7(j,k,ib)
      yntmp8(k)=ya8(j,k,ib)

      enddo

      x2d=Log(x2)/Log(4.)

      call linear(NMAX,x2a,yntmp1,x2d,ymtmp1(j))
      call linear(NMAX,x2a,yntmp2,x2d,ymtmp2(j))
      call linear(NMAX,x2a,yntmp3,x2d,ymtmp3(j))
      call linear(NMAX,x2a,yntmp4,x2d,ymtmp4(j))
      call linear(NMAX,x2a,yntmp5,x2d,ymtmp5(j))
      call linear(NMAX,x2a,yntmp6,x2d,ymtmp6(j))
      call linear(NMAX,x2a,yntmp7,x2d,ymtmp7(j))
      call linear(NMAX,x2a,yntmp8,x2d,ymtmp8(j))

      enddo

      call linear(MMAX,x1a,ymtmp1,x1,res1a(ib))
      call linear(MMAX,x1a,ymtmp2,x1,res2a(ib))
      call linear(MMAX,x1a,ymtmp3,x1,res3a(ib))
      call linear(MMAX,x1a,ymtmp4,x1,res4a(ib))
      call linear(MMAX,x1a,ymtmp5,x1,res5a(ib))
      call linear(MMAX,x1a,ymtmp6,x1,res6a(ib))
      call linear(MMAX,x1a,ymtmp7,x1,res7a(ib))
      call linear(MMAX,x1a,ymtmp8,x1,res8a(ib))

c      print*,b1a(ib),bb,res3a(ib)

      enddo

      call linear(NB,b1a,res1a,bb,res(1))
      call linear(NB,b1a,res2a,bb,res(2))
      call linear(NB,b1a,res3a,bb,res(3))
      call linear(NB,b1a,res4a,bb,res(4))
      call linear(NB,b1a,res5a,bb,res(5))
      call linear(NB,b1a,res6a,bb,res(6))
      call linear(NB,b1a,res7a,bb,res(7))
      call linear(NB,b1a,res8a,bb,res(8))
      

       return
       end

C Linear interpolation
C Vadim Guzey, 06/10/2003

      subroutine linear(n,xa,ya,x,y)
      
      implicit double precision (a-h,o-z)
      integer n,ns,i

      dimension xa(n),ya(n)

C just in case there is some rounding error
      y=0.0
      
C Later added by T. Teckentrup      
      ns=0

      dif=abs(x-xa(1))

      do i=1,n-1

C Finding the closest index i
      dift=abs(x-xa(i))   
      if (dift.le.dif) then
      ns=i
      dif=dift
      endif
      enddo

C x could be greater or smaller than xa(i)

      if (x.gt.xa(ns)) then   

      y1=ya(ns)
      y2=ya(ns+1)
      x1=xa(ns)
      x2=xa(ns+1) 

      y=((y1-y2)*x+(y2*x1-y1*x2))/(x1-x2)

      else

      y1=ya(ns-1)
      y2=ya(ns)
      x1=xa(ns-1)
      x2=xa(ns) 

      y=((y1-y2)*x+(y2*x1-y1*x2))/(x1-x2)   

      endif

      return
      end

      
C Integration

      SUBROUTINE GAUSAB( X, W, A, B, NPT )

C=======================================================================
C                               /b
C     Gaussian quadratures for  | f(x)dx.   NPT must be even & <100.
C                              /a
C
C=======================================================================


      REAL*8 X(1), W(1)

      REAL*8 XX, WW, C1, C2
      COMMON / GAUSCM / XX(100), WW(100)
C-----------------------------------------------------------------------
      NPT2 = NPT/2

      IF( A .EQ. 0. .AND. B .GT. .5E20 ) THEN
      CALL GAUS(NPT,1)
      DO 5 N = 1, NPT
      X(N) = XX(N)
      W(N) = WW(N)
5     CONTINUE
      RETURN
      END IF

      CALL GAUS(NPT,0)

      C1 = (B-A)/2.
      C2 = (B+A)/2.

      DO 10 N = 1, NPT2
      M = NPT2 + 1 - N
      X(N) = - XX(M)*C1 + C2
      W(N) =   WW(M)*C1
10    CONTINUE

      DO 20 N = NPT2+1, NPT
      M = N - NPT2
      X(N) = XX(M)*C1 + C2
      W(N) = WW(M)*C1
20    CONTINUE

      RETURN

      END
C=======================================================================
        SUBROUTINE GAUS(NN,ITYPE)
C
C       STANDARD SETTINGS ARE: IALF=IBTA=0,0<NN<100,
C               ITYPE=0 : -1 TO +1 POINTS
C               ITYPE=1 :  0 TO INFINITY
C
      IMPLICIT REAL*8 (A-H,O-Z)



      DATA IALF, IBTA / 0, 0 /
C***********************************************************************
C
C     GAUJAC COMPUTES GAUSS - JACOBI INTEGRATION WEIGHTS AND NODES
C
C***********************************************************************
C
      DIMENSION X(100),A(100),B(100),C(100)
      DIMENSION XX(100), WW(100)
      DIMENSION G(100)
C
        COMMON/GAUSCM/XX,WW
C
      PI=3.141592653589793238D0
C
      ISUM=IALF+IBTA
      IDIF=IBTA-IALF
      FNUM=ISUM*IDIF
      ALF=IALF
      BTA=IBTA
C
C     CALCULATE COEFFICIENTS REQUIRED BY JACOBI
C
      DO 20 I=1,100
      I2=I+I
      IM=I-1
      DEN=(ISUM+I2)*(ISUM+I2-2)
      IF(ISUM.EQ.0) GO TO 67
      B(I)=FNUM/DEN
      GO TO 68
67    B(I)=0.D0
68    FNUMP=4*IM*(IALF+IM)*(IBTA+IM)*(ISUM+IM)
      DENP=(ISUM+I2-1)*(ISUM+I2-2)**2*(ISUM+I2-3)
      IF(((ISUM.EQ.0).OR.(ISUM.EQ.1)).AND.(I.EQ.1)) GO TO 69
      C(I)=FNUMP/DENP
      GO TO 20
69    C(I)=0.D0
20    CONTINUE

C-----------------------------------------------------------------------
C     For outputting quadratures
C-----------------------------------------------------------------------
C      TYPE 34
C34    FORMAT(1X,'FILENAME FOR INTERMEDIATE RESULTS =? (A10)')
C      ACCEPT 35,FNAMX
C35    FORMAT(A10)
C      OPEN(UNIT=3,FILE=FNAMX,STATUS='NEW')
C-----------------------------------------------------------------------

      CALL JACOBI(NN,X,G,ALF,BTA,B,C,EPS,CSX,CSA,TSX,TSA)

C-----------------------------------------------------------------------
C     For outputting quadratures
C-----------------------------------------------------------------------
C      TYPE 70
C70    FORMAT(1X,'OUTPUT FILENAME ? (A10)')
C      ACCEPT 71,FNAM
C71    FORMAT(A10)
C      OPEN(UNIT=7,FILE=FNAM,STATUS='NEW')
C-----------------------------------------------------------------------

C     OUTPUT DATA STATEMENT OF WEIGHTS AND NODES
C
      IF(ITYPE.EQ.0) GO TO 130
      DO 100 I=1,NN
      II=NN+1-I
C     XX(I)=(1.0D+00+X(II))/(1.0D+00-X(II))
C     WW(I)=2.0D+00*G(II)/((1.0D+00-X(II))**(IALF+2)*
C    1     (1.0D+00+X(II))**IBTA)
      ZZ=PI*(1.0D+00+X(II))/4.0D+00
      XX(I)=SIN(ZZ)/COS(ZZ)
      WW(I)=G(II)*PI*0.25D+00/(COS(0.25D+00*PI*(1.0D+00
     1      +X(II)))**2)

C-----------------------------------------------------------------------
C      WRITE(7,120)I,I,XX(I),WW(I)
120   FORMAT(6X,6HDATAX(,I2,6H,M),W(,I2,4H,M)/,D22.15,1H,,
     1      D22.15,1H/)
C-----------------------------------------------------------------------

100   CONTINUE
        RETURN
 130  CONTINUE
      NI=NN/2
      DO 140 I=1,NI
      II=NI+1-I
      XX(I)=X(II)
      WW(I)=G(II)

C-----------------------------------------------------------------------
C      WRITE(7,120)I,I,XX(I),WW(I)
C-----------------------------------------------------------------------

140   CONTINUE
      RETURN
      END


      SUBROUTINE JACOBI(NN,X,A,ALF,BTA,B,C,EPS,CSX,CSA,TSX,TSA)
      IMPLICIT REAL*8 (A-H,O-Z)
C


C***********************************************************************
C
C     CALCULATES THE ZEROS X(I) OF THE NN-TH ORDER JACOBI POLYNOMIAL
C     PN(ALF,BTA) FOR THE SEGMENT (-1,1)
C     THE LARGEST ZERO WILL BE STORED IN X(1).
C     ALSO CALCULATES THE CORRESPONDING COEFFICIENTS A(I) OF THE
C     NN-TH ORDER GAUSS-JACOBI QUADRATURE FORMULA OF DEGREE 2*NN-1.
C     THIS SUBROUTINE MUST BE GIVEN THE COEFFICIENTS:
C
C                (ALF+BTA)(BTA-ALF)
C        B(N) = -----------------------------
C               (ALF+BTA+2N)(ALF+BTA+2N-2)
C
C                4(N-1)(ALF+N-1)(BTA+N-1)(ALF+BTA+N-1)
C        C(N) = ---------------------------------------------
C               (ALF+BTA+2N-1)(ALF+BTA+2N-2)**2(ALF+BTA+2N-3)
C
C
C     IN THE RECURSION RELATION
C
C        P(N) = (X-B(N))*P(N-1)-C(N)*P(N-2)
C
C     FOR ALL N LESS THAN OR EQUAL TO THE HIGHEST DEGREE NN.
C
C        CSX = CALC SUM X(I)      TSX = TRUE SUM X(I)
C        CSA = CALC SUM A(I)      TSA = TRUE SUM A(I)
C
C***********************************************************************
C
      DIMENSION X(100),A(100),B(100),C(100)
      FN=NN
      CSX=0.D0
      CSA=0.D0
      BETA=EXP(FLGAMA(ALF+1.D0)+FLGAMA(BTA+1.D0)-FLGAMA(ALF+BTA+2.D0))
      CC=2.**(ALF+BTA+1.)*BETA
      TSX=FN*(BTA-ALF)/(ALF+BTA+2.*FN)
      TSA=CC
      DO 1 J=2,NN
1     CC=CC*C(J)
      DO 12 I=1,NN
      IF(I-1)12,2,3
C        LARGEST ZERO
2     AN=ALF/FN
      BN=BTA/FN
      R1=(1.D0+ALF)*(2.78/(4.+FN*FN)+.768*AN/FN)
      R2=1.D0+1.48*AN+.96*BN+.452*AN*AN+.83*AN*BN
      XT=1.D0-R1/R2
      GO TO 11
3     IF(I-2)12,4,5
C        SECOND ZERO
4     R1=(4.1+ALF)/((1.+ALF)*(1.+.156*ALF))
      R2=1.D0+.06*(FN-8.)*(1.D0+.12*ALF)/FN
      R3=1.D0+.012*BTA*(1.D0+.25*ABS(ALF))/FN
      RATIO=R1*R2*R3
      XT=XT-RATIO*(1.D0-XT)
      GO TO 11
5     IF(I-3) 12,6,7
C        THIRD ZERO
6     R1=(1.67+.28*ALF)/(1.D0+.37*ALF)
      R2=1.D0+.22*(FN-8.D0)/FN
      R3=1.D0+8.*BTA/((6.28D0+BTA)*FN*FN)
      RATIO=R1*R2*R3
      XT=XT-RATIO*(X(1)-XT)
      GO TO 11
7     IF(NN-I-1)10,9,8
C        MIDDLE ZEROS
8     XT=3.*X(I-1)-3.*X(I-2)+X(I-3)
      GO TO 11
C        SECOND LAST ZERO
9     R1=(1.D0+.235*BTA)/(.766D0+.119*BTA)
      R2=1./(1.D0+.639*(FN-4.D0)/(1.D0+.71*(FN-4.D0)))
      R3=1./(1.D0+20.*ALF/((7.5D0+ALF)*FN*FN))
      RATIO=R1*R2*R3
      XT=XT+RATIO*(XT-X(I-2))
      GO TO 11
C        LAST ZERO
10    R1=(1.D0+.37*BTA)/(1.67D0+.28*BTA)
      R2=1./(1.D0+.22*(FN-8.D0)/FN)
      R3=1./(1.D0+8.*ALF/((6.28D0+ALF)*FN*FN))
      RATIO=R1*R2*R3
      XT=XT+RATIO*(XT-X(I-2))
C
11    CALL GSROOT(XT,NN,ALF,BTA,DPN,PN1,B,C,EPS)
      X(I)=XT
      A(I)=CC/(DPN*PN1)
C-----------------------------------------------------------------------
C      WRITE(3,20)ALF,BTA,NN,I,XT,A(I)
C-----------------------------------------------------------------------

      CSX=CSX+XT
      CSA=CSA+A(I)
12    CONTINUE

C-----------------------------------------------------------------------
C      WRITE(3,22)CSX,CSA,TSX,TSA
20    FORMAT(1X,2F6.2,2I3,2(1X,D26.18),1X,(1X,D26.18))
22    FORMAT(1H0,/,5X,'CSX =',D25.18,5X,'CSA =',D25.18,/,
     1 5X,'TSX =',D25.18,5X,'TSA =',D25.18)
C-----------------------------------------------------------------------

      RETURN
      END

      SUBROUTINE GSROOT(X,NN,ALF,BTA,DPN,PN1,B,C,EPS)
      IMPLICIT REAL*8 (A-H,O-Z)



C***********************************************************************
C
C     IMPROVES THE APPROXIMATE ROOT X
C     IN ADDITION WE ALSO OBTAIN:
C      DPN = DERIVATIVE OF P(N) AT X
C      PN1 = VALUE OF P(N-1) AT X
C
C***********************************************************************
C
      DIMENSION B(100),C(100)
      ITER=0
1     ITER=ITER+1
      CALL RECR(P,DP,PN1,X,NN,ALF,BTA,B,C)
      D=P/DP
      X=X-D
      IF(ABS(D)-EPS)3,3,2
2     IF(ITER-10)1,3,3
3     DPN=DP
      RETURN
      END

      SUBROUTINE RECR(PN,DPN,PN1,X,NN,ALF,BTA,B,C)
      IMPLICIT REAL*8 (A-H,O-Z)



      DIMENSION B(100),C(100)
      P1=1.D0
      P=X+(ALF-BTA)/(ALF+BTA+2.D0)
      DP1=0.D0
      DP=1.D0
      DO 1 J=2,NN
      Q=(X-B(J))*P-C(J)*P1
      DQ=(X-B(J))*DP+P-C(J)*DP1
      P1=P
      P=Q
      DP1=DP
      DP=DQ
1     CONTINUE
      PN=P
      DPN=DP
      PN1=P1
      RETURN
      END

       FUNCTION FLGAMA(W)
      IMPLICIT REAL*8 (A-H,O-Z)




C
C***********************************************************************
C
C     CALCULATES LOG(BASE E)GAMMA(W) FOR W REAL AND GAMMA(W) POSITIVE.
C     USES STIRLING'S APPROXIMATION
C     ACCURATE TO ABOUT 12 SIGNIFICANT PLACES
C
C***********************************************************************
C
      PI=3.141592653589793238D0
      X=W
      M=0
      FK=-1.D0
      IF(X-.5D0)1,1,2
C        W LESS EQ .5
1     M=1
      XPI=X*PI
      X=1.D0-X
2     FK=FK+1.D0
      IF(X+FK-6.D0)2,2,3
3     Z=X+FK
      ZZ=Z*Z
C        LOG GAMMA(Z), Z GREATER 6.
      Y=(Z-.5)*LOG(Z)-Z+.9189385332047+(((((((-10861851./ZZ
     1 +2356200.)/ZZ-704820.)/ZZ+309400.)/ZZ-218790.)/ZZ
     2 +291720.)/ZZ-1021020.)/ZZ+30630600.)/Z/367567200.
      IF(FK)6,6,4
4     IK=FK
      DO 5 I=1,IK
      FK=FK-1.
5     Y=Y-LOG(X+FK)
6     IF(M)7,11,7
7     P=PI/SIN(XPI)
      IF(P)8,8,10
8     WRITE(6,9)W
9     FORMAT(2X,6HGAMMA(D11.4,13H) IS NEGATIVE)
      Y=0.
      GO TO 11
10    Y=LOG(P)-Y
11    FLGAMA=Y
      RETURN
      END
c***************************************************************************