C***********************************************************************
      PROGRAM freq
C Calcula la frecuencia del pendulo en funcion de la energia 
C y la aproxmación a la frecuencia en la vecindad de la separatriz
C***********************************************************************
      IMPLICIT REAL*8(a-h,o-z)
      CHARACTER arch*40
      WRITE(*,'(a,$)')' Archivo: '
      READ(*,'(a)')arch
      OPEN(50,FILE=arch)
      pi=4.d0*datan(1.d0)
C rango de energías y se toman 10000 valores
      x1=-1.d0
      x2=3.d0
      dx=(x2-x1)/10000. 
C se define la frecuencia de pequeñas oscilaciones
      om0=1.d0
C se calcula la frecuencia para cada energía, x=1 es la separatriz
      DO i=1,10000
      x=x1+(i-1)*dx
C s es k^2 para oscilaciones
      s=(x+1.d0)/2.d0
C frecuencia aproximada en la vecindad de la separatriz
      w=x-1.d0
      IF(w.ne.0.d0) omap=pi*om0/log(32.d0/(abs(w)))
      IF(s.LT.1.d0)THEN
C xk(s) es la integral elíptica en función de k^2
      z=xk(s)
      om=pi*om0/z/2.d0
      ELSE
C k debe ser menor a 1
       IF(s.NE.1.d0)THEN
C para rotaciones k = 1/k_osc
       z=xk(1/s)
       wr=DSQRT(s)*om0
       om= pi*wr/z/2.d0
       ENDIF
      ENDIF 
      WRITE(50,*)x,om,omap
      ENDDO 
      END 

C**********************************************************************
      FUNCTION xk(xm)
C Error en la estimacion: 2.d-8
C**********************************************************************
       IMPLICIT REAL*8(a-h,o-z)
       a0=1.38629436112
       a1=0.09666344259
       a2=0.03590092383
       a3=0.03742563713
       a4=0.01451196212
       b0=0.5
       b1=0.12498593597
       b2=0.06880248576
       b3=0.03328355346
       b4=0.00441787012
       xm1=1.d0-xm
       xk=a0+a1*xm1+a2*xm1**2+a3*xm1**3+a4*xm1**4+(b0+
     % b1*xm1+b2*xm1**2+b3*xm1**3+b4*xm1**4)*DLOG(1/xm1)
       RETURN
       END



