*     ANDOYER_PELADO.F : RESOLUCION ANALITICA DEL ANDOYER 
*
      IMPLICIT REAL*8 (A-H,K-Z)
      REAL*8 Y(2),DYDX(2)
      COMMON /CTE/ TWOPI,CERO,UNO,DOS,UNO2,TRES2,RAD,ERROR,BASE
      COMMON /COE/ A,B,EPS,T
      EXTERNAL DERIVS

C     ESPECIFICACION DE ARCHIVOS.

      OPEN (2,FILE='polin.dat',STATUS='UNKNOWN')
      OPEN (3,FILE='cubica.dat',STATUS='UNKNOWN')
      OPEN (4,FILE='soluciones.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
      ERROR = 1.0D-13

C     CONDICIONES INICIALES EN J Y TETA
      DJ    = 0.01
      TETA  = 20.0*cero
      TETA  = TETA*RAD 
       write(*,*)'ondiciones iniciales', dj,teta 
C     COEFICIENTES DEL HAMILTONIANO DE ANDOYER
c      A   = -0.003
c      B   =  0.04
c      EPS =  0.004

      A    = -0.00001
      B    = 0.01
      EPS  = 0.001
      T    = 0.01
      

C     VALORES DEL KERNEL PARA LA CONDICION INICIAL.
      F430 = A*DJ + B*DJ*DJ + EPS*T*DSQRT(ABS(DOS*DJ))*DCOS(TETA)
      WRITE (*,*) 'F430:',F430
     
C     FACTORES DE ESCALA.
      A1  = (4.0*EPS*T/B)**(0.333333333)
      A1ES = -3.0*B*A1*A1/4.0D0
      H1ES =  0.25*B*(A1**4)
      WRITE(*,*)'H1ES', H1ES
      ALFA  = A/A1ES
      C1 = F430/H1ES
      write(*,*)'A1ES', A1ES
      WRITE (*,*) 'a1, C1,ALFA:',a1,C1,ALFA

C     POLINOMIO CARACTERISTICO.
      DI = DJ/(A1**2)
C      DI = 0.287381654170384
      POL  =  DOS*DI - (C1 + 3.0*ALFA*DI - 4.0*DI*DI)**2
C      WRITE(*,*)'POL EN LA RAIZ',POL
C      STOP
      POL1 =  (DOS-6.0*C1*ALFA) + DOS*(8.0*C1-9.0*ALFA**2)*DI 
      POL1 = POL1 + 72.0*ALFA*DI**2 - 64.0*DI**3
      POL2 = (16.0*C1 - 18*ALFA**2)+ 144.0*ALFA*DI -192*DI**2
     
C     INVARIANTES DE LA ELIPTICA DE WEIERSTRASS
      G2 = UNO3*(8.0*C1 + 4.5*ALFA**2)**2 - 12.0*ALFA
      G3 = ((8.0*C1 + 4.5*ALFA**2)**3)/27.0D0
      G3 = G3 - DOS*ALFA*(8.0*C1 + 4.5*ALFA**2) + 4.0D0
      WRITE (*,*)'INVARIANTES  G2,G3:',G2,G3
     
C     CRITERIO DE LIBRACION (LIBRA SI DELTA < 0).
      DELTA = G2**3 - 27.0*G3*G3
      WRITE(*,*)'VALOR DE DELTA:', DELTA

C     GRAFICO DEL POLINOMIO CARACTERISTICO
c      DIG = -1
      PASO = 0.001
      DO I = 0,1000
         DIG =  PASO*DFLOAT(I)
         POL  =  DOS*DIG - (C1 + 3.0*ALFA*DIG - 4.0*DIG*DIG)**2
         WRITE(2,*) DIG, POL
      END DO
      
C     RESOLUCION NUMERICA DEL POLINOMIO CARACTERISTICO POL.
      DIN1= 0.15
 11   POL   =  DOS*DIN1 - (C1 + 3.0*ALFA*DIN1 - 4.0*DIN1**2)**2
      POL1  =  (DOS-6.0*C1*ALFA) + DOS*(8.0*C1-9.0*ALFA**2)*DIN1
      POL1  = POL1 + 72.0*ALFA*DIN1**2 - 64.0*DIN1**3
      DIN2  = DIN1 - POL/POL1
      POLI2 = DOS*DIN2 - (C1 + 3.0*ALFA*DIN2 - 4.0*DIN2**2)**2
      if(Dabs(POLI2).gt.1.0E-10) then
         DIN1 = DIN2 
         go to 11
      end if
      WRITE(*,*)'raiz del polinomio',DIN2
      DI0 = DIN1

C     GRAFICO DE LA ECUACION CUBICA 
      PASO = 0.005
      DO J = 0,1000  
         EE = -2.0+ PASO*DFLOAT(J)
         ECU  =  4.0*ee**3-g2*ee-g3
         WRITE(3,*) EE, ECU
      END DO
      
C     RESOLUCION DE LA RESOLVENTE CUBICA 4E^3 - G2E -G3 =0
      EE_1=2.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 
         go to 9
      end if
      WRITE(*,*) 'raiz de la cubica resolvente', EE_2
      EE1 = EE_2

      IF (DELTA.LT.CERO) THEN
         ETA = DSQRT(3.0*EE1*EE1 - G2/4.0)
         KA  = DSQRT(UNO2 - 3.0*EE1/ETA/4.0)
         WRITE (*,*) 'ETA,KA',ETA,KA
      ELSE
         WRITE (*,*) 'CIRCULA!!'
      END IF
      
C     AMPLITUD DE LIBRACION.
      AP = DOS - DOS*(C1+3.0*ALFA*DI0-4.0*DI0*DI0)*(3.0*ALFA-8.0*DI0)
      AP = AP/4.0D0/ETA
      BP = -DOS*(3.0*ALFA-8.0*DI0)*(3.0*ALFA-8.0*DI0) +
     *     DOS*(C1+3.0*ALFA*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
      write(*,*)'DIC', DIC
      
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
      stop

CCCC  LAZO SOBRE EL TIEMPO. CALCULO DE LA ORBITA ASOCIADA.
c      TETA0 = TETA
      Y(1) = DJ
      Y(2) = TETA 
C     
      DELT = T0/100.0
      SIGMAX = CERO
      DO 10 IN=0,1000
         TIEMPO = CERO +  DFLOAT(IN)*DELT
         CALL DERIVS (Y,DYDX)
         CALL RK4    (Y,DYDX,2,TIEMPO,DELT,Y,DERIVS)
c     CALL RK4    (Y,DYDX,2,TIEMPO,pass,Y,DERIVS)
         DJ    = Y(1)
         TETA  = Y(2)
         TETA = DMOD(TETA,TWOPI)
         IF (TETA.LT.CERO) TETA = TETA + TWOPI
         IF (TETA.GT.PI)   TETA = TETA - TWOPI
         
         KNUM = DSQRT(ABS(DOS*DJ))*DCOS(TETA)
         HNUM = DSQRT(ABS(DOS*DJ))*DSIN(TETA)
         
         
C     SOLUCION ANALITICA.
         W1    = TWOPI*TIEMPO/T0
         COU   = BP/AP       
         DIA   = DI0 + UNO2*QQ -UNO2*QQ*DCOS(W1)
         DIA   = DIA - UNO2*COU *Q2 + UNO2*COU *Q2*DCOS(DOS*W1)
         STETA = QQ*DSIN(W1)+ (UNO/8.0/DI0- UNO2*COU)*Q2*DSIN(DOS*W1)
         STETA = UNO2*TWOPI/PER/DSQRT(DOS*DIC)*STETA
         DIA1   = DI0 + UNO2*QQ -UNO2*QQ*DCOS(W1)
         STETA1 = QQ*DSIN(W1)
         STETA1 = UNO2*TWOPI/PER/DSQRT(DOS*DIC)*STETA1
         
         
C     LA SOLUTION ANALITICA.
         DJA  = DIA*A1*A1
         T1CA = DASIN(STETA)
         IF (T1CA.LT.CERO) T1CA = T1CA + TWOPI
         KA = DSQRT(DOS*DJA)*DCOS(T1CA)
         HA = DSQRT(DOS*DJA)*DSIN(T1CA)
         
         DJA1  = DIA1*A1*A1
         T1CA1 = DASIN(STETA1)
         IF (T1CA1.LT.CERO) T1CA1 = T1CA1 + TWOPI
         KA1 = DSQRT(DOS*DJA1)*DCOS(T1CA1)
         HA1 = DSQRT(DOS*DJA1)*DSIN(T1CA1)
         
         WRITE (4,100)  KA1, HA1, KNUM, HNUM
         WRITE (*,*)  KA1, HA1, KNUM, HNUM
     
 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(2),DYDX(2)
      COMMON /CTE/ TWOPI,CERO,UNO,DOS,UNO2,TRES2,RAD,ERROR,BASE
      COMMON /COE/ A,B,EPS,T
C
      DJ   = Y(1)
      TETA   = Y(2)

C     HAMILTONIANO.
      F = A* DJ + B*DJ**2 + EPS*T*DSQRT(DOS*DJ)*DCOS(TETA)

C     ECUACIONES DE MOVIMIENTO.
      DYDX(1) = - EPS*T*DSQRT(abs(DOS*DJ))*DSIN(TETA)
      DYDX(2) = A + DOS*B*DJ + EPS*T*DCOS(TETA)/DSQRT(abs(DOS*DJ))
C      DYDX(2) =  DYDX(2) 
      
      RETURN
C
      END

C*******************************************************************
      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
