CC****************************************************************
      PROGRAM pen
CC--- Calcula las trayectorias en el espacio de fases del péndulo
CC****************************************************************
      IMPLICIT REAL*8(a-h,o-z)

C-- Se usa un vector de dos dimensiones y=(q,p) ----    
      DIMENSION y(2)

C-- Rutina externa donde se escriben las ecuaciones de Hamilton --
      EXTERNAL fcn

C-- Archivo de datos y su lectura ---------
      OPEN(21,FILE='pen.dat')

C-- Energía inicial, final, valor inicial del ángulo y el número de
C-- valores de energía a considerar
      READ(21,*)ei,ef,xi,ne

C-- Tiempo inicial y final a integrar y el paso de integración ---
      READ(21,*)t0,tf,paso

C-- Parámetro para reducir la cantidad de salidas a imprimir, a cada
C-- paso de tiempo, dos pasos, etc.
      READ(21,*)icada
      CLOSE(21)

C-- Archivo de salida de datos ------------
      OPEN(3,FILE='pen.out')

      pi=4.d0*datan(1.d0)

C-- Se define el paso en energías -----------
      IF(ei.NE.ef)THEN
      epaso=(ef-ei)/ne
      ELSE
      epaso=1.d10
      ENDIF

C-- Lazo que barre las energias ----------
      DO k0=1,ne

C-- Cond. iniciales de la orbita --------------
       x0=xi*pi/180.d0 
       e0=ei+(k0-1)*epaso
       potencial=-DCOS(x0)
       empot=e0-potencial
       IF(empot.LT.0.d0)GOTO 1

C-- Impulso inicial a partir del ángulo inicial y la energía
       p0=DSQRT(2.d0*empot)

C-- y(1), y(2) son las variables a integrar, y(1) el ángulo, y(2) el
C-- impulso

       y(1)=x0
       y(2)=p0
      
C-- Indice que avanza el tiempo      
      i=0

C-- x es el tiempo inicial y xend el tiempo final entre los cuales se
C-- integra
            xend=t0

C-- Integracion ---------------------------------------
            DO WHILE(xend.LT.tf)
C-- Parámetros de la rutina de integración
               n=2
               eps=1.d-15
               hmax=1.d0
               h=.1d0

               x=t0+i*paso
               xend=x+paso

c-- Se llama a la rutina de integración, devuelve y(1), y(2) luego de
C-- integrar un intervalo xend-x = paso
          CALL dopri8(n,fcn,x,y,xend,eps,hmax,h)

C-- El icada es  para hacer que el archivo de salida
C-- contenga menos datos, si icada = 1, imprime todos

                IF(MOD(i,icada).EQ.0)THEN
                 phi=DMOD(y(1),2.d0*pi)
                 IF(phi.LT.-pi)phi=phi+2.d0*pi 
                 IF(phi.GT.pi)phi=phi-2.d0*pi 
         WRITE(3,198)phi,y(2),-y(2)
 198     FORMAT(F14.8,1x,F14.8,1x,F14.8)
                  ENDIF
               i=i+1
           ENDDO
         
    1    ENDDO           !lazo en e 
       close(3)
       END
C**********************************************************************
C-- Ecuaciones de Hamilton a integrar: dy = dy/dt
      SUBROUTINE fcn(n,x,y,dy)
C**********************************************************************
      IMPLICIT REAL*8(a-h,o-z)
      DIMENSION y(n),dy(2)

      dy(1)=y(2)
      dy(2)=-DSIN(y(1))


      RETURN
      END
C**********************************************************************
C Rutina de integración, RK orden 8 Prince & Dormand
      SUBROUTINE DOPRI8(N,FCN,X,Y,XEND,EPS,HMAX,H)
C**********************************************************************
      IMPLICIT REAL*8(A-H,O-Z)
      EXTERNAL fcn
      REAL*8 K1(51),K2(51),K3(51),K4(51),K5(51),K6(51),K7(51)
     *,Y(N),Y1(51)
      LOGICAL REJECT
      COMMON/STAT/NFCN,NSTEP,NACCPT,NREJCT
      COMMON/COEF/C2,C3,C4,C5,C6,C7,C8,C9,C10,C11,C12,C13,
     *A21,A31,A32,A41,A43,A51,A53,A54,A61,A64,A65,A71,A74,A75,A76,
     *A81,A84,A85,A86,A87,A91,A94,A95,A96,A97,A98,A101,A104,A105,A106,
     *A107,A108,A109,A111,A114,A115,A116,A117,A118,A119,A1110,A121,
     *A124,A125,A126,A127,A128,A129,A1210,A1211,A131,A134,A135,A136,
     *A137,A138,A139,A1310,A1311,B1,B6,B7,B8,B9,B10,B11,B12,B13,
     *BH1,BH6,BH7,BH8,BH9,BH10,BH11,BH12
C     UROUND:REL. MASCHINENGENAUIGKEIT
C     ANGEPASST AUF PE8/32 UND SIEMENS
      DATA INI/0/,NMAX/3000/,UROUND/.2220446049D-15/
      IF(INI.EQ.0)THEN
         INI=1
         CALL COEFST
      ENDIF
      POSNEG=DSIGN(1.D0,XEND-X)
      HMAX=DABS(HMAX)
      H=DMIN1(DMAX1(1.D-10,DABS(H)),HMAX)
      H=DSIGN(H,POSNEG)
      EPS=DMAX1(EPS,13.D0*UROUND)
      REJECT=.FALSE.
      NACCPT=0
      NREJCT=0
      NFCN=0
      NSTEP=0
 1    IF(NSTEP.GT.NMAX.OR.X+.03D0*H.EQ.X)GO TO 79
      IF((X-XEND)*POSNEG+UROUND.GT.0.D0)RETURN
      IF((X+H-XEND)*POSNEG.GT.0.D0)H=XEND-X
      CALL FCN(N,X,Y,K1)
 2    CONTINUE
      NSTEP=NSTEP+1
      DO 22 I=1,N
 22   Y1(I)=Y(I)+H*A21*K1(I)
      CALL FCN(N,X+C2*H,Y1,K2)
      DO 23 I=1,N
 23   Y1(I)=Y(I)+H*(A31*K1(I)+A32*K2(I))
      CALL FCN(N,X+C3*H,Y1,K3)
      DO 24 I=1,N
 24   Y1(I)=Y(I)+H*(A41*K1(I)+A43*K3(I))
      CALL FCN(N,X+C4*H,Y1,K4)
      DO 25 I=1,N
 25   Y1(I)=Y(I)+H*(A51*K1(I)+A53*K3(I)+A54*K4(I))
      CALL FCN(N,X+C5*H,Y1,K5)
      DO 26 I=1,N
 26   Y1(I)=Y(I)+H*(A61*K1(I)+A64*K4(I)+A65*K5(I))
      CALL FCN(N,X+C6*H,Y1,K6)
      DO 27 I=1,N
 27   Y1(I)=Y(I)+H*(A71*K1(I)+A74*K4(I)+A75*K5(I)+A76*K6(I))
      CALL FCN(N,X+C7*H,Y1,K7)
      DO 28 I=1,N
 28   Y1(I)=Y(I)+H*(A81*K1(I)+A84*K4(I)+A85*K5(I)+A86*K6(I)+A87*K7(I))
      CALL FCN(N,X+C8*H,Y1,K2)
      DO 29 I=1,N
 29   Y1(I)=Y(I)+H*(A91*K1(I)+A94*K4(I)+A95*K5(I)+A96*K6(I)+A97*K7(I)
     *      +A98*K2(I))
      CALL FCN(N,X+C9*H,Y1,K3)
      DO 30 I=1,N
 30   Y1(I)=Y(I)+H*(A101*K1(I)+A104*K4(I)+A105*K5(I)+A106*K6(I)
     *       +A107*K7(I)+A108*K2(I)+A109*K3(I))
      DO 61 I=1,N
      Y11S=A111*K1(I)+A114*K4(I)+A115*K5(I)+A116*K6(I)+A117*K7(I)
     *     +A118*K2(I)+A119*K3(I)
      Y12S=A121*K1(I)+A124*K4(I)+A125*K5(I)+A126*K6(I)+A127*K7(I)
     *     +A128*K2(I)+A129*K3(I)
      K4(I)=A131*K1(I)+A134*K4(I)+A135*K5(I)+A136*K6(I)+A137*K7(I)
     *     +A138*K2(I)+A139*K3(I)
      K5(I)=B1*K1(I)+B6*K6(I)+B7*K7(I)+B8*K2(I)+B9*K3(I)
      K6(I)=BH1*K1(I)+BH6*K6(I)+BH7*K7(I)+BH8*K2(I)+BH9*K3(I)
      K2(I)=Y11S
 61   K3(I)=Y12S
      CALL FCN(N,X+C10*H,Y1,K7)
      DO 31 I=1,N
 31   Y1(I)=Y(I)+H*(K2(I)+A1110*K7(I))
      CALL FCN(N,X+C11*H,Y1,K2)
      XPH=X+H
      DO 32 I=1,N
 32   Y1(I)=Y(I)+H*(K3(I)+A1210*K7(I)+A1211*K2(I))
      CALL FCN(N,XPH,Y1,K3)
      DO 33 I=1,N
 33   Y1(I)=Y(I)+H*(K4(I)+A1310*K7(I)+A1311*K2(I))
      CALL FCN(N,XPH,Y1,K4)
      NFCN=NFCN+13
      DO 35 I=1,N
      K5(I)=Y(I)+H*(K5(I)+B10*K7(I)+B11*K2(I)+B12*K3(I)+B13*K4(I))
 35   K6(I)=Y(I)+H*(K6(I)+BH10*K7(I)+BH11*K2(I)+BH12*K3(I))
      ERR=0.D0
      DO 41 I=1,N
      DENOM=DMAX1(1.D-6,DABS(K5(I)),DABS(Y(I)),2.D0*UROUND/EPS)
 41   ERR=ERR+((K5(I)-K6(I))/DENOM)**2
      ERR=DSQRT(ERR/DFLOAT(N))
      FAC=DMAX1((1.D0/6.D0),DMIN1(3.D0,(ERR/EPS)**(1.D0/8.D0)/.9D0))
      HNEW=H/FAC
      IF(ERR.GT.EPS)GO TO 51
      NACCPT=NACCPT+1
      DO 44 I=1,N
 44   Y(I)=K5(I)
      X=XPH
      IF(DABS(HNEW).GT.HMAX)HNEW=POSNEG*HMAX
      IF(REJECT)HNEW=POSNEG*DMIN1(DABS(HNEW),DABS(H))
      REJECT=.FALSE.
      H=HNEW
      GO TO 1
 51   REJECT=.TRUE.
      H=HNEW
      IF(NACCPT.GE.1)NREJCT=NREJCT+1
      NFCN=NFCN-1
      GO TO 2
C 79   WRITE(6,979)X
C 979  FORMAT(' EXIT OF DOPRI8 AT X=',D16.7/)
 79   RETURN
      END
C**********************************************************************
      SUBROUTINE COEFST
C**********************************************************************
      IMPLICIT REAL*8(A-H,O-Z)
      COMMON/COEF/C2,C3,C4,C5,C6,C7,C8,C9,C10,C11,C12,C13,
     *A21,A31,A32,A41,A43,A51,A53,A54,A61,A64,A65,A71,A74,A75,A76,
     *A81,A84,A85,A86,A87,A91,A94,A95,A96,A97,A98,A101,A104,A105,A106,
     *A107,A108,A109,A111,A114,A115,A116,A117,A118,A119,A1110,A121,
     *A124,A125,A126,A127,A128,A129,A1210,A1211,A131,A134,A135,A136,
     *A137,A138,A139,A1310,A1311,B1,B6,B7,B8,B9,B10,B11,B12,B13,
     *BH1,BH6,BH7,BH8,BH9,BH10,BH11,BH12
      C2=1.D0/18.D0
      C3=1.D0/12.D0
      C4=1.D0/8.D0
      C5=5.D0/16.D0
      C6=3.D0/8.D0
      C7=59.D0/400.D0
      C8=93.D0/200.D0
      C9=5490023248.D0/9719169821.D0
      C10=13.D0/20.D0
      C11=1201146811.D0/1299019798.D0
      C12=1.D0
      C13=1.D0
      A21=C2
      A31=1.D0/48.D0
      A32=1.D0/16.D0
      A41=1.D0/32.D0
      A43=3.D0/32.D0
      A51=5.D0/16.D0
      A53=-75.D0/64.D0
      A54=-A53
      A61=3.D0/80.D0
      A64=3.D0/16.D0
      A65=3.D0/20.D0
      A71=29443841.D0/614563906.D0
      A74=77736538.D0/692538347.D0
      A75=-28693883.D0/1125.D6
      A76=23124283.D0/18.D8
      A81=16016141.D0/946692911.D0
      A84=61564180.D0/158732637.D0
      A85=22789713.D0/633445777.D0
      A86=545815736.D0/2771057229.D0
      A87=-180193667.D0/1043307555.D0
      A91=39632708.D0/573591083.D0
      A94=-433636366.D0/683701615.D0
      A95=-421739975.D0/2616292301.D0
      A96=100302831.D0/723423059.D0
      A97=790204164.D0/839813087.D0
      A98=800635310.D0/3783071287.D0
      A101=246121993.D0/1340847787.D0
      A104=-37695042795.D0/15268766246.D0
      A105=-309121744.D0/1061227803.D0
      A106=-12992083.D0/490766935.D0
      A107=6005943493.D0/2108947869.D0
      A108=393006217.D0/1396673457.D0
      A109=123872331.D0/1001029789.D0
      A111=-1028468189.D0/846180014.D0
      A114=8478235783.D0/508512852.D0
      A115=1311729495.D0/1432422823.D0
      A116=-10304129995.D0/1701304382.D0
      A117=-48777925059.D0/3047939560.D0
      A118=15336726248.D0/1032824649.D0
      A119=-45442868181.D0/3398467696.D0
      A1110=3065993473.D0/597172653.D0
      A121=185892177.D0/718116043.D0
      A124=-3185094517.D0/667107341.D0
      A125=-477755414.D0/1098053517.D0
      A126=-703635378.D0/230739211.D0
      A127=5731566787.D0/1027545527.D0
      A128=5232866602.D0/850066563.D0
      A129=-4093664535.D0/808688257.D0
      A1210=3962137247.D0/1805957418.D0
      A1211=65686358.D0/487910083.D0
      A131=403863854.D0/491063109.D0
      A134=-5068492393.D0/434740067.D0
      A135=-411421997.D0/543043805.D0
      A136=652783627.D0/914296604.D0
      A137=11173962825.D0/925320556.D0
      A138=-13158990841.D0/6184727034.D0
      A139=3936647629.D0/1978049680.D0
      A1310=-160528059.D0/685178525.D0
      A1311=248638103.D0/1413531060.D0
      B1=14005451.D0/335480064.D0
      B6=-59238493.D0/1068277825.D0
      B7=181606767.D0/758867731.D0
      B8=561292985.D0/797845732.D0
      B9=-1041891430.D0/1371343529.D0
      B10=760417239.D0/1151165299.D0
      B11=118820643.D0/751138087.D0
      B12=-528747749.D0/2220607170.D0
      B13=1.D0/4.D0
      BH1=13451932.D0/455176623.D0
      BH6=-808719846.D0/976000145.D0
      BH7=1757004468.D0/5645159321.D0
      BH8=656045339.D0/265891186.D0
      BH9=-3867574721.D0/1518517206.D0
      BH10=465885868.D0/322736535.D0
      BH11=53011238.D0/667516719.D0
      BH12=2.D0/45.D0
      RETURN
      END

