*     ANDOYER.F : HAMILTONIANA DE ANDOYER DE DOS ARMONICOS VIA SYLVIO.
*
      IMPLICIT REAL*8 (A-H,K-Z)
      REAL*8 Y(4),DYDX(4)
      REAL*8 BIN(0:20,0:20),RPER(0:2,0:15,0:15,-15:15,-15:15)
      REAL*8 S1(0:5,0:5),T2(0:5,0:5),DII(0:5,0:5)
      COMPLEX*16 CI,EMT2(0:5,-5:5,0:5),ET1(0:5,-5:5,0:5)
      COMPLEX*16 EMT1(0:5,-5:5,0:5)
      COMPLEX*16 EJ1(0:5,-5:5,0:5),PROD(0:5,-5:5,0:5)
      COMPLEX*16 A_IN(0:5,-5:5,0:5),B_IN(0:5,-5:5,0:5)
      COMPLEX*16 F1_1(0:5,-5:5,0:5)
      COMPLEX*16 FUNC,APROX2,TERME,aprox3,ASD,ASD1,ASD2
      COMMON /CTE/ TWOPI,CERO,UNO,DOS,UNO2,TRES2,RAD,ERROR,BASE
      COMMON /ORB/ SEMI1,SEMI2,ALFA,ALFA0,E1,E2,TITA1,TITA2,L
      COMMON /FUN/ RPER
      COMMON /COE/ A,B,EPS,T,NU11,BB,MM,E20,R1M1_1,R1M1_2
      EXTERNAL DERIVS

C     ESPECIFICACION DE ARCHIVOS.
      OPEN (1,FILE='rper.coef',STATUS='OLD')
      OPEN (2,FILE='elem_medios.dat',STATUS='OLD')
      OPEN (3,FILE='andoyer.dat',STATUS='UNKNOWN')
      
C     INICIALIZACION DE VARIABLES GLOBALES.
      TWOPI = 8.0D0*DATAN(1.0D0)
      CERO  = 0.0D0
      UNO   = 1.0D0
      DOS   = 2.0D0
      UNO2  = UNO/DOS
      UNO3  = UNO/3.0D0
      TRES2 = 3.0D0/2.0D0
      RAD   = TWOPI/360.0
      PI    = UNO2*TWOPI
      GG    = 0.01720209895**2
      ERROR = 1.0D-13
      BASE  = 999.9D0
      CI    = DCMPLX(CERO,UNO)

C     CALCULO DE COEFICIENTES BINOMIALES.
      CALL BINOM (BIN)

C     MASAS Y PARAMETROS DE LA RESONANCIA.
      P = DOS
      Q = UNO
      ALFA0  = (P/(P+Q))**(0.66666666)
      SEMI2  = 5.2D0
      E2     = 0.05
      SEMI10 = SEMI2*ALFA0
      GJ     = CERO
      MASS0  = UNO
      MASS2  = 1.0D-3
      MU     = GG*(MASS0 + MASS2)
      ENE2   = DSQRT(MU)/(SEMI2**1.5)
      EPS    = GG*MASS2
            
C     LEE COEFICIENTES RESONANTES FINALES DEL DESARROLLO GLOBAL DE R.
      RPER = CERO
 1    READ (1,*,END=2) IA,J,IK,IU,IT,RPER(IA,J,IK,IU,IT)
       RPER(IA,J,IK,IU,IT) = -RPER(IA,J,IK,IU,IT)/SEMI2
       GOTO 1
 2    CONTINUE

C     LEE ELEMENTOS MEDIOS DEL ASTEROIDE.
      READ (2,*) SEMI1,E1,TITA1,TITA2
      TITA1 = TITA1*RAD
      TITA2 = TITA2*RAD
      ALFA = SEMI1/SEMI2

C     CALCULO DE LOS MOMENTOS CANONICOS.
      LRES = DSQRT(MU*SEMI10)
      L    = DSQRT(MU*SEMI1)
      G    = L*DSQRT(UNO - E1*E1)
      TAU  = UNO2*(ENE2 - GJ)*(LRES - 3.0*L) - GJ*G
      DJ1  = L - G
      DJ2  = (ENE2*G + TAU)/(ENE2-GJ)
      DJ3  = LRES
      NU11 = -12.0*(MU**2)/(DJ3**4)
      WRITE (*,*) 'DJ1,DJ2',DJ1,DJ2
      WRITE (*,*) 'NU11',NU11

C     CALCULAMOS COEFICIENTES PRINCIPALES.
      CALL COEFS (R10,R01,BB,MM,R1M1,R1M1_1,R1M1_2)
      BB   = BB/DSQRT(DJ3)
      BETA = MM/BB
      
C     CALCULAMOS SESSIN.
      K = DSQRT(DOS*DJ1)*COS(TITA1) + BETA*E2*COS(TITA2)
      H = DSQRT(DOS*DJ1)*SIN(TITA1) + BETA*E2*SIN(TITA2)
      S = DJ2 + UNO2*BETA*BETA*E2*E2
      S = S - BETA*E2*(K*COS(TITA2) + H*SIN(TITA2))

C     MAS PARAMETROS QUE SURGEN DE LA TRANSFORMACION DE SESSIN.
      A = NU11*S
      B = UNO2*NU11
      T = BB
      TETA1C = DATAN2(H,K)
      DJ = UNO2*(K*K + H*H)
	write(*,*)A, B, EPS, T
	STOP
C     CAMBIO DE SIGNO DE A,B & T. 
      A = -A
      B = -B
C     PARA ACERTAR SIGNOS, CAMBIAMOS THETA1C --> THETA1C + PI
      T = T
C     (AHORA TENEMOS B >0 Y EPS*TAU>0)

C     VALORES DEL KERNEL PARA LA CONDICION INICIAL.
      F430 = A*DJ + B*DJ*DJ + EPS*T*DSQRT(ABS(DOS*DJ))*DCOS(TETA1C)

      WRITE (*,*) 'A,B:',A,B
      WRITE (*,*) 'EPS,T:',EPS,T
      WRITE(*,*)'DJ:', DJ
      WRITE (*,*) 'F430:',F430
      WRITE (*,*) 'BB,MM:',BB,MM
      
C     FACTORES DE ESCALA.
      A1  = ABS(4.0*EPS*T/B)**(0.333333333)
      A1S = -3.0*B*A1*A1/4.0D0
      H1S =  0.25*B*(A1**4)
      WRITE (*,*) 'A1:',A1

C     OTROS PARAMETROS MAS (Y VAN...).
C      WRITE (*,*) 'J1+J2,S+J',DJ1+DJ2,S+DJ
C     SESSIN ESTA BIEN DEFINIDA, O POR LO MENOS DA LO MISMO!
      
      R  = A/A1S
      C1 = F430/H1S

C     POLINOMIO CARACTERISTICO.
      DJ = UNO2*(K*K + H*H)
      TETA1C = DATAN2(H,K)
      DI = DJ/(A1**2)
      POL = DOS*DI - (C1 + 3.0*R*DI - 4.0*DI*DI)**2
      G2 = UNO3*(8.0*C1 + 4.5*R*R)**2 - 12.0*R
      G3 = ((8.0*C1 + 4.5*R*R)**3)/27.0D0
      G3 = G3 - DOS*R*(8.0*C1 + 4.5*R*R) + 4.0D0
      WRITE (*,*) 'C1,R:',C1,R
      WRITE (*,*)'INVARIANTES  G2,G3:',G2,G3
     
C     RESOLUCION NUMERICA DEL POLINOMIO CARACTERISTICO POL.
      FUNC0 = BASE
      DIF = UNO
      DI0 = CERO
      DO WHILE (DIF.GT.CERO)
       DI0 = DI0 + 1.0d-4
       POL = DOS*DI0 - (C1 + 3.0*R*DI0 - 4.0*DI0*DI0)**2
       IF (FUNC0.EQ.BASE) THEN
        FUNC0 = POL
       ELSE
        IF (POL*FUNC0.LT.CERO) GOTO 12
       END IF
       FUNC0 = POL
      END DO
 12   CONTINUE
      WRITE (*,*) 'DI0,DJ0:',DI0,DI0*(A1**2)
      
C     CRITERIO DE LIBRACION (LIBRA SI DELTA < 0).
      DELTA = G2**3 - 27.0*G3*G3
      WRITE(*,*)'VALOR DE DELTA:', DELTA

      ee_1=0.
 9    ee_2= ee_1 - (4.0*ee_1**3-g2*ee_1-g3)/(12.*ee_1**2-g2)
      if(abs(4.0*ee_2**3-g2*ee_2-g3).gt.1.0d-9) then
         ee_1 = ee_2 
         print*,ee_2
         go to 9
      end if


      IF (DELTA.LT.CERO) THEN
C     RESOLUCION DE LA ECUACION RESOLVENTE.
       IF (G2.GT.CERO) THEN
        ARG = DSQRT(27.0*G3*G3/(G2**3))
        WRITE (*,*) 'ARG',ARG
        DACOSHARG = DLOG(ARG + DSQRT(ARG*ARG - UNO))
        GI  = UNO3*DACOSHARG
        EE1 = DSQRT(G2/3.0)*DCOSH(GI)
       ELSE
        ARG = DSQRT(-27.0*G3*G3/(G2**3))
        DASINHARG = DLOG(ARG + DSQRT(ARG*ARG + UNO))
        GI  = UNO3*DASINHARG
        EE1 = DSQRT(-G2/3.0)*DSINH(GI)
       END IF
       ETA = DSQRT(3.0*EE1*EE1 - G2/4.0)
       KA  = DSQRT(UNO2 - 3.0*EE1/ETA/4.0)
       WRITE (*,*) 'ETA,KA',ETA,KA
       WRITE (*,*) 'EE1, ee_2 ,GI',EE1,ee_2, GI
      ELSE
       WRITE (*,*) 'CIRCULA!!'
      END IF

C     AMPLITUD DE LIBRACION.
      AP = DOS - DOS*(C1+3.0*R*DI0-4.0*DI0*DI0)*(3.0*R-8.0*DI0)
      AP = AP/4.0D0/ETA
      BP = -DOS*(3.0*R-8.0*DI0)*(3.0*R-8.0*DI0) +
     *     DOS*(C1+3.0*R*DI0-4.0*DI0*DI0)*8.0
      BP = (BP + 24.0*(ETA-EE1))/24.0/ETA
      QQ = AP/(UNO - BP)
      WRITE (*,*) 'AP,BP',AP,BP
      WRITE (*,*) 'QQ,QQ [deg]',QQ,QQ/RAD

C     CENTRO DE LIBRACION.
      DIC = DI0 + UNO2*QQ

C     INTEGRAL ELIPTICA K(k).
      KA2 = KA*KA
      KELIP = UNO + 0.25*KA2 + 9.0*KA2*KA2/64.0 + 25.0*(KA2**3)/256.0
      KELIP = KELIP + 1225.0*(KA2**4)/16384.0 + 3969.0*(KA2**5)/65536.0
      KELIP = UNO2*PI*KELIP
      
C     PERIODO DE LIBRACION (ESCALEADO).
      PER = 4.0*KELIP/DOS/DSQRT(ETA)

C     PERIODO DE LIBRACION (EN TIEMPO FISICO).
      T0 = PER*(4.0/A1/A1/B)
      WRITE (*,*) 'PER,T0',PER,T0

C     ACCION.
      Q2 = QQ*QQ
      Q3 = QQ*Q2
      ACC = -A1*A1*Q2*(3.0*R-8.0*DI0)*PI/(4.0*PER)
      ACC = ACC - 4.0*Q3 + BP*BP*(3.0*R-8.0*DI0)*Q2*Q2/8.0/AP/AP
      ACC = ACC - BP*Q2*Q2/DOS/AP - BP*BP*Q2*Q3/4.0/AP/AP
      ACC = ACC + BP*BP*BP*(3.0*R-8.0*DI0)*Q2*Q3/8.0/AP/AP/AP
      WRITE (*,*) 'ACC',ACC
      WRITE (*,*) 'ERROR RELATIVO EN ACCION',
     *     BP*BP*BP*(3.0*R-8.0*DI0)*Q2*Q3/8.0/AP/AP/AP/ACC
      
C     DESARROLLO EN SERIE DE SENO DE TETA1C: 
C     SIN(TETA1C) = SUM S1(I,J) QQ*I SIN(J*w1).
C********************
C     SI USAMOS TIEMPO FISICO, PONER T0 EN VEZ DE PER
C********************
      COU = BP/AP
      S1 = CERO
      S1(1,1) = PI/T0/DSQRT(DOS*DIC)*(4.0/A1/A1/B)
      S1(2,2) = (PI/per/DSQRT(DOS*DIC))*(4.0/A1/A1/B)*DSQRT(DOS)/16.0
      S1(2,2) =  S1(2,2)*(UNO/DIC - 4.0*COU)
      
C     DESARROLLO EN SERIE DE TETA2C (PARTE NO CIRCULATORIA): 
      FAC = NU11*A1*A1*T0/TWOPI
      WRITE (*,*) 'FAC',FAC
      T2 = CERO
      T2(1,1) = -UNO2
      T2(2,2) =COU/4.0
      DO I = 0,5
      DO J = 0,5
       T2(I,J) = FAC*T2(I,J)
      END DO
      END DO
C     NOTA: RECORDAR QUE FALTA EL TERMINO SUMANDO XI.

C     DESARROLLO EN SERIE DE DI-DI0.
      DII = CERO
      DII(1,0) =  UNO2
      DII(1,1) = -UNO2
      DII(2,0) = -COU/8.0
      DII(2,2) =  COU/8.0

C     FRECUENCIA FUNDAMENTAL DE TITA2.
      APROX = DI0
      DO I = 0,5
       APROX = APROX + DII(I,0)*(QQ**I)
      END DO
      M = NU11*S + NU11*A1*A1*APROX
      WRITE (*,*) 'PER_TITA2',TWOPI/M,-TWOPI/(NU11*(DJ1+DJ2))


CCCC  LAZO SOBRE EL TIEMPO. CALCULO DE LA ORBITA ASOCIADA.
      E20 = E2
      TITA10 = TITA1
      TITA20 = TITA2
      Y(1) = DJ1
      Y(2) = DJ2
      Y(3) = TITA1
      Y(4) = TITA2
C     
      DELT = T0/100.0
      SIGMAX = CERO
c      TIEMPO = CERO
      DO 10 IN=0,1000
         TIEMPO = CERO +  DFLOAT(IN)*DELT
C         TIEMPO = CERO,10*T0,DELT
       CALL DERIVS (Y,DYDX)
       CALL RK4    (Y,DYDX,4,TIEMPO,DELT,Y,DERIVS)
       DJ1   = Y(1)
       DJ2   = Y(2)
       TITA1 = Y(3)
       TITA2 = Y(4)
       TITA1 = DMOD(TITA1,TWOPI)
       TITA2 = DMOD(TITA2,TWOPI)
       IF (TITA1.LT.CERO) TITA1 = TITA1 + TWOPI
       IF (TITA2.LT.CERO) TITA2 = TITA2 + TWOPI
       IF (TITA1.GT.PI)   TITA1 = TITA1 - TWOPI

       Knum = dsqrt(dos*dj1)*dcos(tita1)+ beta*e2*dcos(tita2)
       Hnum = dsqrt(dos*dj1)*dsin(tita1)+ beta*e2*dsin(tita2)

C     SOLUCION ANALITICA.
       W1 = TWOPI*TIEMPO/T0
       
C     CALCULO DE DI,TETA1C.
       DIA = CERO
       T1CA = CERO
       T2CA = CERO
       DO I = 0,1
       DO J = 0,1
        DIA  = DIA  + DII(I,J)*(QQ**I)*DCOS(J*W1)
        T1CA = T1CA +  S1(I,J)*(QQ**I)*DSIN(J*W1)
        T2CA = T2CA +  T2(I,J)*(QQ**I)*DSIN(J*W1)
       END DO
       END DO
       DIA = DIA + DI0
C
       TITA2A = TITA20 + M*TIEMPO + T2CA
       TITA2A = DMOD(TITA2A,TWOPI)
       IF (TITA2A.LT.CERO) TITA2A = TITA2A + TWOPI
       
C     LA SOLUTION ANALITICA.
       DJA  = DIA*A1*A1
       T1CA = DASIN(T1CA)
       IF (T1CA.LT.CERO) T1CA = T1CA + TWOPI
       KA = DSQRT(DOS*DJA)*DCOS(T1CA)
       HA = DSQRT(DOS*DJA)*DSIN(T1CA)
C
       TITA1A = DATAN2(HA-E2*BETA*DSIN(TITA2A),KA-E2*BETA*DCOS(TITA2A))
       TITA1A = DMOD(TITA1A,TWOPI) + TITA10
       IF (TITA1A.LT.CERO) TITA1A = TITA1A + TWOPI
       IF (TITA1A.GT.PI)   TITA1A = TITA1A - TWOPI
C
       DJ2A = S - UNO2*E2*E2*BETA**2
       DJ2A = DJ2A + BETA*E2*(KA*DCOS(TITA2A) + HA*DSIN(TITA2A))
       DJ1A = S + DJA - DJ2A
C
       XI = TITA20 + M*TIEMPO


C     VARIABLES REGULARES
       
       KK  = DSQRT(DOS*DJ1)*DCOS(TITA1)
       HH  = DSQRT(DOS*DJ1)*DSIN(TITA1)
       KKA = DSQRT(DOS*DJ1A)*DCOS(TITA1A)
       HHA = DSQRT(DOS*DJ1A)*DSIN(TITA1A)

C     SALIDA DE RESULTADOS.
c       WRITE (3,100) TIEMPO,KK, HH, KKA,HHA
       WRITE (3,100) TIEMPO,Knum,Hnum,KA,HA
c      WRITE (3,100) TIEMPO,DJ,DJA,TETA1C/RAD,T1CA/RAD
C      WRITE (3,100) TIEMPO,TITA1/RAD,TITA1A/RAD,TITA2/RAD,TITA2A/RAD
C      WRITE (3,100) TIEMPO,DJ1+DJ2-S,DJA
c      WRITE (3,100) TIEMPO,DJ1,DJ1A,DJ2,DJ2A
c      WRITE (3,100) TIEMPO,KK,HH,KKA,HHA
C     
 10   CONTINUE
C     
 100  FORMAT (1P9E15.6)
C
      END


C*********************************************************************
C     SUBRUTINAS


      SUBROUTINE DERIVS (Y,DYDX)
C     
      IMPLICIT REAL*8 (A-H,K-Z)
      REAL*8 Y(4),DYDX(4)
      COMMON /CTE/ TWOPI,CERO,UNO,DOS,UNO2,TRES2,RAD,ERROR,BASE
      COMMON /COE/ A,B,EPS,T,NU11,BB,MM,E2,R1M1_1,R1M1_2
C
c       bb = -bb
      DJ1   = Y(1)
      DJ2   = Y(2)
      TITA1 = Y(3)
      TITA2 = Y(4)

C     HAMILTONIANO.
      F = UNO2*NU11*(DJ1+DJ2)**2 + EPS*BB*DSQRT(DOS*DJ1)*DCOS(TITA1)
      F = F + EPS*MM*E2*DCOS(TITA2)

C     ECUACIONES DE MOVIMIENTO.
      DYDX(1) = -EPS*BB*DSQRT(DOS*DJ1)*DSIN(TITA1)
      DYDX(2) = -EPS*MM*E2*DSIN(TITA2)
      DYDX(3) =  NU11*(DJ1+DJ2) + EPS*BB*DCOS(TITA1)/DSQRT(DOS*DJ1)
      DYDX(4) =  NU11*(DJ1+DJ2)
      DYDX(3) = -DYDX(3)
      DYDX(4) = -DYDX(4)
C
      DO I = 1,4
       DYDX(I) = -DYDX(I)
      END DO
C
      RETURN
C
      END


      SUBROUTINE COEFS (R10,R01,BB,MM,R1M1,R1M1_1,R1M1_2)
C
      IMPLICIT REAL*8 (A-H,K-Z)
      REAL*8 R(0:2,0:15,0:15,-15:15,-15:15)
      COMMON /CTE/ TWOPI,CERO,UNO,DOS,UNO2,TRES2,RAD,ERROR,BASE
      COMMON /ORB/ SEMI1,SEMI2,ALFA,ALFA0,E1,E2,TITA1,TITA2,L
      COMMON /FUN/ R
C
      IMAX = 15

C     FUNCION PERTURBADORA.
      DELA = ALFA - ALFA0
      IF (DELA.LT.1.0E-15) DELA = 1.0E-15
      r00 = cero
      R01 = CERO
      r10 = cero
      r20 = cero
      r02 = cero
      r11 = cero
      r1m1 = cero
      difb = cero
      difm = cero
      r20_1 = cero
      r20_2 = cero
      db_10 = cero
      r11_1 = cero
      r11_2 = cero
      r1m1_1 = cero
      r1m1_2 = cero
      r00_0 = cero
      r00_1 = cero
      r00_2 = cero
      dm_00 = cero
      dm_10 = cero
      dm_20 = cero
      dr10ds = cero
      s10_1 = cero
      s10_2 = cero
      dr01ds = cero
      s01_1 = cero
      s01_2 = cero
      BB    = CERO
      MM    = CERO
C      
      DO IA = 0,2
      DO J  = 0,IMAX
      DO IK = 0,IMAX
      DO IU = -IMAX,IMAX
      DO IT = -IMAX,IMAX
       FAC  = (E1**J)*(E2**IK)*(DELA**IA)
       
       if (abs(it).eq.0.and.abs(iu).eq.0) then
        if (abs(j).ne.0.or.abs(ik).ne.0) then
         r00 = r00 + FAC*R(IA,J,IK,IU,IT)
         if (j.eq.0) then
          r00_0 = r00_0 + FAC*R(IA,J,IK,IU,IT)
         end if
         if (j.eq.2) then
          r00_1 = r00_1 + FAC*R(IA,J,IK,IU,IT)
         end if
         if (j.eq.4) then
          r00_2 = r00_2 + FAC*R(IA,J,IK,IU,IT)
         end if
        end if
       end if
       
       if (abs(it).eq.0.and.abs(iu).eq.1) then
        r10 = r10 + FAC*R(IA,J,IK,IU,IT)
        dr10ds = dr10ds + J*FAC*R(IA,J,IK,IU,IT)/sqrt(l)
        if (j.ne.1) then
         s10_1 = s10_1 + J*FAC*R(IA,J,IK,IU,IT)/sqrt(l)
        end if
        if (j.ne.3) then
         s10_2 = s10_2 + J*FAC*R(IA,J,IK,IU,IT)/sqrt(l)
        end if
       end if
       
       if (abs(it).eq.0.and.abs(iu).eq.1) then
        if (j.ne.1) then
         difb = difb + FAC*R(IA,J,IK,IU,IT)
         if (j.eq.3) then
          db_10 = db_10 + FAC*R(IA,J,IK,IU,IT)
         end if
        end if
       end if
       
       if (it.eq.iu.and.abs(it).eq.1) then
        if (ik.ne.1) then
         difm = difm + FAC*R(IA,J,IK,IU,IT)
         if (j.eq.0) then
          dm_00 = dm_00 + FAC*R(IA,J,IK,IU,IT)
         end if
         if (j.eq.2) then
          dm_10 = dm_10 + FAC*R(IA,J,IK,IU,IT)
         end if
         if (j.eq.4) then
          dm_20 = dm_20 + FAC*R(IA,J,IK,IU,IT)
         end if
        end if
       end if
       
       if (abs(it).eq.0.and.abs(iu).eq.2) then
        r20 = r20 + FAC*R(IA,J,IK,IU,IT)
        if (j.eq.2) then
         r20_1 = r20_1 + FAC*R(IA,J,IK,IU,IT)
        end if
        if (j.eq.4) then
         r20_2 = r20_2 + FAC*R(IA,J,IK,IU,IT)
        end if
       end if
       
       if (abs(it).eq.1.and.abs(iu).eq.0) then
        r1m1 = r1m1 + FAC*R(IA,J,IK,IU,IT)
        if (j.eq.1) then
         r1m1_1 = r1m1_1 + FAC*R(IA,J,IK,IU,IT)
        end if
        if (j.eq.3) then
         r1m1_2 = r1m1_2 + FAC*R(IA,J,IK,IU,IT)
        end if
       end if
       
       if (it.eq.1.and.iu.eq.2) then
        r11 = r11 + FAC*R(IA,J,IK,IU,IT)
        if (j.eq.1) then
         r11_1 = r11_1 + FAC*R(IA,J,IK,IU,IT)
        end if
        if (j.eq.3) then
         r11_2 = r11_2 + FAC*R(IA,J,IK,IU,IT)
        end if
       end if
       if (it.eq.-1.and.iu.eq.-2) then
        r11 = r11 + FAC*R(IA,J,IK,IU,IT)
        if (j.eq.1) then
         r11_1 = r11_1 + FAC*R(IA,J,IK,IU,IT)
        end if
        if (j.eq.3) then
         r11_2 = r11_2 + FAC*R(IA,J,IK,IU,IT)
        end if
       end if
       
       if (it.eq.iu.and.abs(it).eq.1) then
        r01 = r01 + FAC*R(IA,J,IK,IU,IT)
        dr01ds = dr01ds + J*FAC*R(IA,J,IK,IU,IT)/sqrt(l)
        if (j.ne.2) then
         s01_1 = s01_1 + J*FAC*R(IA,J,IK,IU,IT)/sqrt(l)
        end if
        if (j.ne.4) then
         s01_2 = s01_2 + J*FAC*R(IA,J,IK,IU,IT)/sqrt(l)
        end if
        
       end if
       
       if (it.eq.iu.and.abs(it).eq.2) then
        r02 = r02 + FAC*R(IA,J,IK,IU,IT)
       end if
      END DO
      END DO
      END DO
      END DO
      END DO
C
      DO IA = 0,2
      DO J  = 0,IMAX
      DO IK = 0,IMAX
      DO IU = -IMAX,IMAX
      DO IT = -IMAX,IMAX
       FAC  = (E1**J)*(E2**IK)*(DELA**IA)
       if (abs(it).eq.0.and.abs(iu).eq.1.and.j.eq.1) then
        BB = BB + FAC*R(IA,J,IK,IU,IT)
       end if
       if (it.eq.iu.and.abs(it).eq.1.and.ik.eq.1.and.j.eq.0) then
        MM = MM + FAC*R(IA,J,IK,IU,IT)
       end if
      END DO
      END DO
      END DO
      END DO
      END DO
C
      bb = bb/e1
      mm = mm/e2
C
      RETURN
C
      END


      SUBROUTINE RK4 (Y,DYDX,IN,X,H,YOUT,DERIVS)
C     
      IMPLICIT REAL*8 (A-H,K-Z)
      PARAMETER (INMAX=10)
      REAL*8 Y(IN),DYDX(IN),YOUT(IN),YT(INMAX),DYT(INMAX),DYM(INMAX)
C     
      HH = 0.5D0*H
      H6 = H/6.0D0
      XH = X + HH
      DO 11 I = 1,IN
       YT(I) = Y(I) + HH*DYDX(I)
 11   CONTINUE
      CALL DERIVS (YT,DYT)
      DO 12 I = 1,IN
       YT(I) = Y(I) + HH*DYT(I)
 12   CONTINUE
      CALL DERIVS (YT,DYM)
      DO 13 I = 1,IN
       YT(I)  = Y(I) + H*DYM(I)
       DYM(I) = DYT(I) + DYM(I)
 13   CONTINUE
      CALL DERIVS (YT,DYT)
      DO 14 I = 1,IN
       YOUT(I) = Y(I) + H6*(DYDX(I) + DYT(I) + 2.0*DYM(I))
 14   CONTINUE
C     
      RETURN
C     
      END


      SUBROUTINE BINOM (BIN)
C     
      IMPLICIT REAL*8 (A-H,K-Z)
      REAL*8 BIN(0:20,0:20)
      COMMON /CTE/ TWOPI,CERO,UNO,DOS,UNO2,TRES2,RAD,ERROR,BASE

C     CALCULO DE LA MATRIZ DE COEFICIENTES BINOMIALES.
      DO I = 0,20
       BIN(I,0) = UNO
       BIN(I,I) = UNO
      END DO
      DO I = 1,19
       DO J = 1,I-1
        BIN(I,J) = BIN(I,J-1)*DFLOAT(I + 1 - J)/DFLOAT(J)
       END DO
      END DO
C
      RETURN
C
      END


      SUBROUTINE PRODUCTO (A,B,PROD)
C     
      IMPLICIT REAL*8 (A-H,K-Z)
      COMPLEX*16 A(0:5,-5:5,0:5),B(0:5,-5:5,0:5),PROD(0:5,-5:5,0:5)
      COMPLEX*16 TERM
      COMMON /CTE/ TWOPI,CERO,UNO,DOS,UNO2,TRES2,RAD,ERROR,BASE
C
      PROD = CERO
      DO I1  =  0,5
      DO J1  = -5,5
      DO IK1 =  0,5
       DO I2  =  0,5
       DO J2  = -5,5
       DO IK2 =  0,5
        I  = I1  + I2
        J  = J1  + J2
        IK = IK1 + IK2
        IMAX = MAX(I,ABS(J),ABS(IK))
        IF (IMAX.LE.5) THEN
         TERM = A(I1,J1,IK1)*B(I2,J2,IK2)
         PROD(I,J,IK) = PROD(I,J,IK) + TERM
        END IF
       END DO
       END DO
       END DO
      END DO
      END DO
      END DO
C
      RETURN
C
      END
