*     HANSEN.F  :  DETERMINES MODIFIED HANSEN COEFFICIENS.
*     
      IMPLICIT REAL*8 (A-H,K-Z)
      INTEGER NMAX,NMAX1,NMAX2
      PARAMETER (NMAX=30,NMAX1=50,NMAX2=40)
      REAL*8 Y(0:NMAX2,0:NMAX2,-NMAX1:NMAX1,-NMAX:NMAX)
      REAL*8 B(-NMAX1:NMAX1,0:NMAX,0:NMAX2,0:NMAX2)
      REAL*8 C(-NMAX1:NMAX1,0:NMAX,0:NMAX2,0:NMAX2)
      COMMON /GEN/ TWOPI,CERO,UNO,DOS,TRES2,G,UNO2,RAD,ERROR
C     
      OPEN (4,FILE='hansen.coef',STATUS='UNKNOWN',FORM='UNFORMATTED')
      
C     DEFINITION OF GLOBAL CONSTANTS.
      TWOPI = 8.0D0*DATAN(1.0D0)
      CERO  = 0.0D0 
      UNO   = 1.0D0
      DOS   = 2.0D0
      TRES  = 3.0D0
      UNO2  = UNO/DOS
      UNO3  = UNO/TRES
      DOS3  = DOS/TRES
      TRES2 = TRES/DOS
      PI    = UNO2*TWOPI
      CIEN  = 100.0D0
      RAD   = PI/180.0D0
      ERROR = 1.0D-13
      CI    = DCMPLX(CERO,UNO)
      
C     INICIALIZATION OF COEFFICIENTS.
      DO IN = -NMAX1,NMAX1
      DO IM =      0,NMAX
      DO J  =      0,NMAX2
      DO IK =      0,NMAX2
       B(IN,IM,J,IK) = CERO
       C(IN,IM,J,IK) = CERO
      END DO
      END DO
      END DO
      END DO
      
C     MODIFIED NEWCOMB OPERATORS.
      CALL NEWCOMB (Y)

C     FINAL COEFFICIENTS.
      DO IN = -NMAX1,NMAX1
      DO IM =      0,NMAX
      DO IK =      0,NMAX2
       IA = MAX(0,IK-IM)
       IB = MAX(0,IM-IK)
       DO IS = 0,NMAX2
        J = ABS(IM-IK) + 2*IS
        IF (J.LE.NMAX2.AND.(IS+IA).LT.NMAX2.AND.(IS+IB).LT.NMAX2) THEN
         B(IN,IM,J,IK) = B(IN,IM,J,IK) + Y(IS+IA,IS+IB,IN,IM)
         C(IN,IM,J,IK) = C(IN,IM,J,IK) + Y(IS+IA,IS+IB,IN,IM)
        END IF
       END DO
       IF (IK.NE.0) THEN
        IA = MAX(0,-IK-IM)
        IB = MAX(0,IM+IK)
        DO IS = 0,NMAX2
         J = ABS(IM+IK) + 2*IS
         IF (J.LE.NMAX2.AND.(IS+IA).LT.NMAX2.AND.(IS+IB).LT.NMAX2) THEN
          B(IN,IM,J,IK) = B(IN,IM,J,IK) + Y(IS+IA,IS+IB,IN,IM)
          C(IN,IM,J,IK) = C(IN,IM,J,IK) - Y(IS+IA,IS+IB,IN,IM)
         END IF
        END DO
       END IF
      END DO
      END DO
      END DO

C     WRITES FINAL COEFFICIENTS IN FILE.
      DO IN = -NMAX1,NMAX1
      DO IM =      0,NMAX
      DO J  =      0,NMAX2
      DO IK =      0,NMAX2
       IF (B(IN,IM,J,IK).NE.CERO.OR.C(IN,IM,J,IK).NE.CERO) THEN
        WRITE (4) IN,IM,J,IK,B(IN,IM,J,IK),C(IN,IM,J,IK)
       END IF
      END DO
      END DO
      END DO
      END DO
C     
      END

      
      SUBROUTINE NEWCOMB (Y)
C     
      IMPLICIT REAL*8 (A-H,K-Z)
      INTEGER NMAX,NMAX1,NMAX2
      PARAMETER (NMAX=30,NMAX1=50,NMAX2=40)
      REAL*8 GAMLN(0:NMAX1)
      REAL*8 Y(0:NMAX2,0:NMAX2,-NMAX1:NMAX1,-NMAX:NMAX)
      COMMON /GEN/ TWOPI,CERO,UNO,DOS,TRES2,G,UNO2,RAD,ERROR
C     
      EXTERNAL BETA
C
      do i=0,nmax2
      do j=0,nmax2
      do ik=-nmax1,nmax1      
      do il=-nmax,nmax
      Y(i,j,ik,il) = CERO
      end do
      end do
      end do
      end do

C     GAMLN.
      DO J = 0,NMAX1
       GAMLN(J) = BETA(1.5D0,DFLOAT(J))
      END DO
      
C     OPERADORES DE NEWCOMB Y(I,J,IN,IM).
      DO IN = -NMAX1,NMAX1
      DO IM =  -NMAX,NMAX
       Y(0,0,IN,IM) = UNO
       Y(1,0,IN,IM) = DFLOAT(IM) - 0.5*DFLOAT(IN)
      END DO
      END DO
C
      DO I  =      1,NMAX2
      DO IN = -NMAX1,NMAX1
      DO IM =  -NMAX,NMAX
       SUMA = CERO
       IF (ABS(IM+1).LE.NMAX) THEN
        SUMA = SUMA + DOS*(2*IM-IN)*Y(I-1,0,IN,IM+1)
       END IF
       IF (I.GT.1.AND.ABS(IM+2).LT.NMAX) THEN
        SUMA = SUMA + (IM-IN)*Y(I-2,0,IN,IM+2)
       END IF
       Y(I,0,IN,IM) = SUMA/(4.0*DFLOAT(I))
      END DO
      END DO
      END DO
C     
      DO J  =      1,NMAX2
      DO I  =      0,NMAX2
      DO IN = -NMAX1,NMAX1
      DO IM =  -NMAX,NMAX
       SUMA = CERO
       IF (ABS(IM-1).LE.NMAX) THEN
        SUMA = SUMA - DOS*(2*IM+IN)*Y(I,J-1,IN,IM-1)
       END IF
       IF (J.GT.1.AND.ABS(IM-2).LE.NMAX) THEN
        SUMA = SUMA - (IM+IN)*Y(I,J-2,IN,IM-2)
       END IF
       IF (I.GT.0) SUMA = SUMA - (I-5*J+4+4*IM+IN)*Y(I-1,J-1,IN,IM)
       DO IK = 2,MIN(I,J)
        SUMA = SUMA + DOS*(I-J+IM)*((-UNO)**IK)*GAMLN(IK)
     *       *Y(I-IK,J-IK,IN,IM)
       END DO
       Y(I,J,IN,IM) = SUMA/(4.0*DFLOAT(J))
      END DO
      END DO
      END DO
      END DO
C     
      RETURN
C     
      END
      
      
      FUNCTION BETA (A,B)
C     
      IMPLICIT REAL*8 (A-H,K-Z)
      COMMON /GEN/ TWOPI,CERO,UNO,DOS,TRES2,G,UNO2,RAD,ERROR
      EXTERNAL GAMMLN
C     
      BETA = CERO
      AB = A - B
      Z  = UNO + AB
      IF (Z.GT.CERO) THEN
       IF (Z.GT.UNO) THEN
        GAMAMB = GAMMLN(Z)
       ELSE
        GAMAMB = LOG(EXP(GAMMLN(UNO+Z))/Z)
       END IF
      ELSE
       EGAMMA = UNO2*TWOPI/EXP(GAMMLN(UNO-Z))/DSIN(UNO2*TWOPI*Z)
       SIGNO = EGAMMA/ABS(EGAMMA)
       GAMAMB = LOG(ABS(EGAMMA))
      END IF
C     
      BETA = SIGNO*DEXP(GAMMLN(A+UNO)-GAMMLN(B+UNO)-GAMAMB)
C     
      RETURN
C     
      END
      
      
      FUNCTION GAMMLN (XX)
C     
      IMPLICIT REAL*8 (A-H,K-Z)
      REAL*8 COF(6)
      COMMON /GEN/ TWOPI,CERO,UNO,DOS,TRES2,G,UNO2,RAD,ERROR
C     
      DATA COF,STP /76.18009172947146,-86.50532032941677,
     *     24.01409824083091,-1.231739572450155,
     *     0.001208650973866179,-0.5395239384953D-5,
     *     2.5066282746310005/
C     
      X = XX
      Y = X
      TMP = X + 5.5
      TMP = (X+0.5)*DLOG(TMP) - TMP
      SER = 1.0000000000190015
      DO J = 1,6
       Y = Y + 1.0
       SER = SER + COF(J)/Y
      END DO
      GAMMLN = TMP + DLOG(STP*SER/X)
C     
      RETURN
C     
      END
      
