*     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 (4,FILE='ando_numerico.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.1
C      DJ     = 0.1
      TETA  = 30.0*cero
      TETA  = TETA*RAD 

      A    = 0.01
      B    = -0.005
      EPS  = 0.1
      T    = -0.1
      
C     VALORES DEL KERNEL PARA LA CONDICION INICIAL.
      F430 = A*DJ + B*DJ*DJ +EPS*T*DSQRT(DOS*DJ)*DCOS(TETA)
      WRITE (*,*) 'F430:',F430
     

CCCC  LAZO SOBRE EL TIEMPO. CALCULO DE LA ORBITA ASOCIADA.
      Y(1) = DJ
      Y(2) = TETA 
C     
      DELT = 1.
      SIGMAX = CERO
      DO 10 IN=0,10000
         TIEMPO = CERO +  DFLOAT(IN)*DELT
       CALL DERIVS (Y,DYDX)
       CALL RK4    (Y,DYDX,2,TIEMPO,DELT,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     SALIDA DE RESULTADOS.

       WRITE (4,*)tiempo,  KNUM, HNUM
C       WRITE (4,*)tiempo, TETA/rad, DJ
C     
 10   CONTINUE
 20   CONTINUE
C     

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     ECUACIONES DE MOVIMIENTO.
      DYDX(1) = EPS*T*DSQRT(DOS*DJ)*DSIN(TETA)
      DYDX(2) = A + DOS*B*DJ + EPS*T*DCOS(TETA)/DSQRT(DOS*DJ)
      
      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
