*      	 RPER.F  :  CALCULA LOS COEFICIENTES RESONANTES DE <R> PARA 
*                   RESONANCIA DE MOV. MEDIOS (P+Q)/P. UTILIZA HANSEN.COEF 
*                   COMO DATO DE ENTRADA. INCLUYE TERMINOS DIRECTOS E 
*                   INDIRECTOS PARA EL PROBLEMA GENERAL DE TRES CUERPOS.
*
        IMPLICIT REAL*8 (A-H,K-Z)
        INTEGER NMAX,NMAX1,NMAX2
	PARAMETER (NMAX=30,NMAX1=50,NMAX2=40)
        REAL*8 AKI(0:100,0:100)
	REAL*8 B(-NMAX1:NMAX1,0:NMAX,0:NMAX2,0:NMAX2)
	REAL*8 C(-NMAX1:NMAX1,0:NMAX,0:NMAX2,0:NMAX2)
	REAL*8 D(0:2,0:15,0:15,-15:15,-15:15)
        COMMON /GEN/ TWOPI,CERO,UNO,DOS,TRES2,G,UNO2,RAD,ERR
	COMMON /RES/ P,Q
C
        OPEN (1,FILE='hansen.coef',STATUS='OLD',FORM='UNFORMATTED')
        OPEN (3,FILE='rper.coef',STATUS='UNKNOWN')

C          DEFINICION DE VARIABLES GLOBALES.
        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
        ERR   = 1.0D-13

C	   ESPECIFICAMOS LA RESONANCIA.
c	WRITE (*,*) '  Resonancia (p+q)/p '
c	WRITE (*,*) '  Ingrese valores de p,q '
c	READ  (*,*) P,Q
        P = 2.0D0
        Q = 1.0D0

C	   CALCULA LOS COEFICIENTES AKI(K,I).
	CALL COEFAKI (AKI)

C          LA PARTE INDIRECTA.
        AKI(1,0) = AKI(1,0) - UNO2

C	   INICIALIZACION DE COEFICIENTES.
	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          LEE COEFICIENTES DE HANSEN.
 1	READ (1,END=2) IN,IM,J,IK,B(IN,IM,J,IK),C(IN,IM,J,IK)
	 GOTO 1
 2	CONTINUE

C	   COEFICIENTES DIRECTOS RESONANTES FINALES.
	ALFA0 = (P/(P+Q))**(0.66666666D0)
	DO IA =  0,2
	WRITE (*,*) IA
        DO IU = -15,15
        DO IL =  0,15
c        DO IL =  0,NMAX
        IJMIN = ABS(IU*P - IL)*cero
	IKMIN = ABS(IU*(P+Q) - IL)*cero
	IT = IU*(P+Q) - IL
	IF (ABS(IT).LE.15) THEN
         DO IJ = IJMIN,15
         DO IK = IKMIN,15
	  D(IA,IJ,IK,IU,IT) = CERO
	  DO I = 0,NMAX - IL
	   J = 2*I + IL
	   IF (J.LT.NMAX1) THEN
	    FAC = UNO
	    DO IO = 1,IA
	     FAC = FAC*DFLOAT(J+1-IO)/DFLOAT(IO)
	    END DO
	    PROD = UNO2*FAC*AKI(IL,I)*(ALFA0**(J-IA))
	    IM = IU*P
	    IN = IU*(P+Q)
            SIGNM = UNO
            SIGNN = UNO
            IF (IM.NE.0) SIGNM = DABS(DFLOAT(IM))/DFLOAT(IM)
	    IF (IN.NE.0) SIGNN = DABS(DFLOAT(IN))/DFLOAT(IN)
	    IF (ABS(IM).LE.NMAX2.AND.ABS(IN).LE.NMAX2) THEN
             PROD1 =    B(J,IL,IJ,ABS(IM)) + SIGNM*C(J,IL,IJ,ABS(IM))
             PROD2 = B(-J-1,IL,IK,ABS(IN)) + SIGNN*C(-J-1,IL,IK,ABS(IN))
	     D(IA,IJ,IK,IU,IT) = D(IA,IJ,IK,IU,IT) + PROD*PROD1*PROD2
	    END IF
	   END IF
          END DO
	 END DO
	 END DO
	END IF
	END DO
	END DO
	END DO

C	   SALIDA DE RESULTADOS (DIRECTOS+INDIRECTOS).
	DO IA = 0,2
        DO IJ = 0,15
        DO IK = 0,15
	DO IU = -15,15
        DO IT = -15,15
         if (abs(D(IA,IJ,IK,IU,IT)).gt.1.0d-10) then
          WRITE (3,100) IA,IJ,IK,IU,IT,D(IA,IJ,IK,IU,IT)
         end if
	END DO
	END DO
	END DO
	END DO
	END DO
C
 100	FORMAT (5I6,1PE28.18)
C
        END


	SUBROUTINE COEFAKI (AKI)
C
        IMPLICIT REAL*8 (A-H,K-Z)
	INTEGER NMAX
        REAL*8 A(0:100),B(0:100),C(0:100),BIN(0:200,0:200)
        REAL*8 AKI(0:100,0:100),COEF(90)
        COMMON /GEN/ TWOPI,CERO,UNO,DOS,TRES2,G,UNO2,RAD,ERR
C
C          CALCULA COEFICIENTES A(J) Y DETERMINA NMAX.
	CALL UNOMX (COEF,NMAX)
	DO I = 1,90
	 A(I-1) = COEF(I)
	END DO

C          PASO A COEFICIENTES INTERMEDIOS B(J).
        DOS9 = DSQRT(DOS/9.0D0)
        SIE9 = -7.0D0/9.0D0
        DO J = 0,NMAX
         B(J) = DOS9*A(J)*(SIE9**J)
	END DO

C          CALCULO DE LA MATRIZ DE COEFICIENTES BINOMIALES.
        FAC = 1.0D-20
        DO I = 0,200
         BIN(I,0) = FAC
         BIN(I,I) = FAC
	END DO
        DO I = 1,199
        DO J = 1,I-1
         BIN(I,J) = BIN(I,J-1)*DFLOAT(I + 1 - J)/DFLOAT(J)
	END DO
	END DO

C          PASO A COEFICIENTES FINALES C(J).
        DOS7 = -DOS/7.0D0
        DO J = 0,NMAX
         C(J) = CERO
         DO IN = J,NMAX
          C(J) = C(J) + B(IN)*BIN(IN,J)*(DOS7**J)/FAC
	 END DO
	END DO

C          CALCULO DE LOS COEFICIENTES HIPER-FINALES AKI(K,I).
        DO IK = 0,NMAX
        GAMAK = (-UNO)**IK
        IF (IK.EQ.0) GAMAK = UNO2*GAMAK
        DO I = 0,NMAX-IK
         AKI(IK,I) = CERO
         JMAX = MIN(2*I,NMAX-IK)
         DO J = I,JMAX
          JIK = 2*J - 2*I + IK
          AKI(IK,I)=AKI(IK,I)+C(J+IK)*BIN(J+IK,JIK)*BIN(JIK,J-I)/FAC/FAC
         END DO
         AKI(IK,I) = AKI(IK,I)*GAMAK
        END DO
	END DO
C
	RETURN
C
        END


	SUBROUTINE UNOMX (COEF,NMAX)
C
        IMPLICIT REAL*8 (A-H,O-Z)
        REAL*8 X(300),Y(300),COEF(90),LISTA(90),COVAR(90,90),AFUNC(90)
        COMMON /GEN/ TWOPI,CERO,UNO,DOS,TRES2,G,UNO2,RAD,ERR
        EXTERNAL FUNCS

C	   INGRESO DE DATOS DE LA EXPANSION.
	MA = 40
	MA = MA + 1
c	DMIN = 0.9
	DMIN = 0.99
c	DMIN = 0.92
	NMAX = MA
C      
        NDATA = 100

C          LAZO DE CALCULO DE LA FUNCION EXACTA.
        DO 10 I = 1,NDATA
         X(I) = -DMIN + (UNO+DMIN-0.01)*DFLOAT(I-1)/99.0D0
         Y(I) = UNO/DSQRT(UNO + X(I))
 10     CONTINUE
 
        MFIT = MA
        NCVM = 60
        DO I = 1,MFIT
         LISTA(I) = I
        END DO

        CALL LFIT (X,Y,NDATA,COEF,MA,LISTA,MFIT,COVAR,NCVM,CHISQ,FUNCS)

        DO 30 I = 1,NDATA
         CALL FUNCS (X(I),AFUNC,MA)
         FA = 0.0D0
         DO 20 J = 1,MA
          FA = FA + COEF(J)*AFUNC(J)
 20      CONTINUE
 30     CONTINUE
C
	RETURN
C
        END


      SUBROUTINE LFIT (X,Y,NDATA,A,MA,LISTA,MFIT,COVAR,NCVM,
     *                                                 CHISQ,FUNCS)
C
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER (MMAX=90)
      REAL*8 X(NDATA),Y(NDATA),A(MA),LISTA(MA)
      REAL*8 COVAR(NCVM,NCVM),BETA(MMAX),AFUNC(MMAX)
      EXTERNAL FUNCS
C
      KK = MFIT + 1
      DO 12 J = 1,MA
      IHIT = 0
      DO 11 K = 1,MFIT
      IF (LISTA(K).EQ.J) IHIT = IHIT + 1
 11   CONTINUE
      IF (IHIT.EQ.0) THEN
      LISTA(KK) = J
      KK = KK + 1
      ELSE IF (IHIT.GT.1) THEN
      WRITE (*,*) 'IMPROPER SET IN LISTA'
      END IF
 12   CONTINUE
      IF (KK.NE.(MA+1)) WRITE (*,*) 'IMPROPER SET IN LISTA'
      DO 14 J = 1,MFIT
      DO 13 K = 1,MFIT
      COVAR(J,K) = 0.0D0
 13   CONTINUE
      BETA(J) = 0.0D0
 14   CONTINUE
      DO 18 I = 1,NDATA
      CALL FUNCS (X(I),AFUNC,MA)
      YM = Y(I)
      IF (MFIT.LT.MA) THEN
      DO 15 J = MFIT+1,MA
      YM = YM - A(LISTA(J))*AFUNC(LISTA(J))
 15   CONTINUE
      END IF
      DO 17 J = 1,MFIT
      WT = AFUNC(LISTA(J))
      DO 16 K = 1,J
      COVAR(J,K) = COVAR(J,K) + WT*AFUNC(LISTA(K))
 16   CONTINUE
      BETA(J) = BETA(J) + YM*WT
 17   CONTINUE
 18   CONTINUE
      IF (MFIT.GT.1) THEN
      DO 21 J = 2,MFIT
      DO 19 K = 1,J-1
      COVAR(K,J) = COVAR(J,K)
 19   CONTINUE
 21   CONTINUE
      END IF
      CALL GAUSSJ (COVAR,MFIT,NCVM,BETA,1,1)
      DO 22 J = 1,MFIT
      A(LISTA(J)) = BETA(J)
 22   CONTINUE
      CHISQ = 0.0D0
      DO 24 I = 1,NDATA
      CALL FUNCS (X(I),AFUNC,MA)
      SUM = 0.0D0
      DO 23 J = 1,MA
      SUM = SUM + A(J)*AFUNC(J)
 23   CONTINUE
      CHISQ = CHISQ + (Y(I) - SUM)**2
 24   CONTINUE
      CALL COVSRT (COVAR,NCVM,MA,LISTA,MFIT)
C
      RETURN
C
      END


      SUBROUTINE COVSRT (COVAR,NCVM,MA,LISTA,MFIT)
C
      IMPLICIT REAL*8 (A-H,O-Z)
      REAL*8 COVAR(NCVM,NCVM),LISTA(MFIT)
C
      DO 12 J = 1,MA-1
      DO 11 I = J+1,MA
      COVAR(I,J) = 0.0D0
 11   CONTINUE
 12   CONTINUE
      DO 14 I = 1,MFIT-1
      DO 13 J = I+1,MFIT
      IF (LISTA(J).GT.LISTA(I)) THEN
      COVAR(LISTA(J),LISTA(I)) = COVAR(I,J)
      ELSE 
      COVAR(LISTA(I),LISTA(J)) = COVAR(I,J)
      END IF
 13   CONTINUE
 14   CONTINUE
      SWAP = COVAR(1,1)
      DO 15 J = 1,MA
      COVAR(1,J) = COVAR(J,J)
      COVAR(J,J) = 0.0D0
 15   CONTINUE
      COVAR(LISTA(1),LISTA(1)) = SWAP
      DO 16 J = 2,MFIT
      COVAR(LISTA(J),LISTA(J)) = COVAR(1,J)
 16   CONTINUE
      DO 18 J = 2,MA
      DO 17 I = 1,J-1
      COVAR(I,J) = COVAR(J,I)
 17   CONTINUE
 18   CONTINUE
C
      RETURN
C
      END


      SUBROUTINE GAUSSJ (A,N,NP,B,M,MP)
C
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER (NMAX=90)
      REAL*8 A(NP,NP),B(NP,MP),IPIV(NMAX),INDXR(NMAX),INDXC(NMAX)
C
      DO 11 J = 1,N
      IPIV(J) = 0
 11   CONTINUE
      DO 22 I = 1,N
      BIG = 0.0D0
      DO 13 J = 1,N
      IF (IPIV(J).NE.1) THEN
      DO 12 K = 1,N
      IF (IPIV(K).EQ.0) THEN
      IF (DABS(A(J,K)).GE.BIG) THEN
      BIG = DABS(A(J,K))
      IROW = J
      ICOL = K
      END IF
      ELSE IF (IPIV(K).GT.1) THEN
      WRITE (*,*) 'SINGULAR MATRIX'
      END IF
 12   CONTINUE
      END IF
 13   CONTINUE
      IPIV(ICOL) = IPIV(ICOL) + 1
C
      IF (IROW.NE.ICOL) THEN
      DO 14 L = 1,N
      DUM = A(IROW,L)
      A(IROW,L) = A(ICOL,L)
      A(ICOL,L) = DUM
 14   CONTINUE
      DO 15 L = 1,M
      DUM = B(IROW,L)
      B(IROW,L) = B(ICOL,L)
      B(ICOL,L) = DUM
 15   CONTINUE
      END IF
      INDXR(I) = IROW
      INDXC(I) = ICOL
      IF (A(ICOL,ICOL).EQ.0.) WRITE (*,*) 'SINGULAR MATRIX'
      PIVINV = 1.0D0/A(ICOL,ICOL)
      A(ICOL,ICOL) = 1.0D0
      DO 16 L = 1,N
      A(ICOL,L) = A(ICOL,L)*PIVINV
 16   CONTINUE
      DO 17 L = 1,M
      B(ICOL,L) = B(ICOL,L)*PIVINV
 17   CONTINUE
      DO 21 LL = 1,N
      IF (LL.NE.ICOL) THEN
      DUM = A(LL,ICOL)
      A(LL,ICOL) = 0.0D0
      DO 18 L = 1,N
      A(LL,L) = A(LL,L) - A(ICOL,L)*DUM
 18   CONTINUE
      DO 19 L = 1,M
      B(LL,L) = B(LL,L) - B(ICOL,L)*DUM
 19   CONTINUE
      END IF
 21   CONTINUE
 22   CONTINUE
      DO 24 L = N,1,-1
      IF (INDXR(L).NE.INDXC(L)) THEN
      DO 23 K = 1,N
      DUM = A(K,INDXR(L))
      A(K,INDXR(L)) = A(K,INDXC(L))
      A(K,INDXC(L)) = DUM
 23   CONTINUE
      END IF
 24   CONTINUE
C
      RETURN
C
      END


      SUBROUTINE FUNCS (X,P,NP)
C
      IMPLICIT REAL*8 (A-H,O-Z)
      REAL*8 P(NP)
C
      P(1) = 1.0D0
      DO 11 J = 2,NP
      P(J) = P(J-1)*X
 11   CONTINUE
C
      RETURN
C
      END




