CC****************************************************************
      PROGRAM stan_lyap 
CC  Calcula el exponente de Lyapunov para el mapa standard
CC  normalizando el vector desviación
CC  Variables del mapa standard en (0,1)
CC****************************************************************
      IMPLICIT REAL*8(a-h,o-z)
      OPEN(21,FILE='stan_lyap.dat')
      READ(21,*)xk
      READ(21,*)N,icada
      READ(21,*)s0,t0
      READ(21,*)ds0,dt0
      CLOSE(21)
      OPEN(31,FILE='stan_lyap.out')
C-----NORMA INICIAL DEL VECTOR DESVIACION
      v0=dsqrt(ds0*ds0+dt0*dt0)
C-----SE NORMALIZA EL VECTOR DESVIACION (INICIAL)
      ds0=ds0/v0;dt0=dt0/v0

      dospi=8.d0*DATAN(1.d0)
      xk1=xk/dospi
      xlyapsum0=0.d0

      DO 4 k=1,N
C-----MAPA
      co=DCOS(dospi*t0)
      se=DSIN(dospi*t0)
      s=s0+xk1*se
      t=t0+s
       t=DMOD(t,1.d0)
       IF (t.lt.0.d0) t=t+1.d0
        s=DMOD(s,1.d0) 
       IF(s.LT.0.d0)s=s+1.d0 
       s0=s
       t0=t
C-----MAPA TANGENTE       
       ds=ds0+dospi*xk1*co*dt0
       dt=dt0+ds
C------SE CALCULA LA NORMA DEL VECTOR DESVIACION
       v=dsqrt(ds*ds+dt*dt)
C------SE CALCULA LA SUMA DEL LOG(NORMA) EN CADA ITERACION
       xlyapsum=xlyapsum0+dlog(v)
C------SE CALCULA EL EXPONENTE DE LYAPUNOV
       xlyap=xlyapsum/k
C------SE NORMALIZA EL VECTOR DESVIACION CADA VEZ 
       dt=dt/v;ds=ds/v
       xlyapsum0=xlyapsum
       ds0=ds
       dt0=dt
      IF (MOD(k,icada).eq.0) WRITE(31,*)k,xlyap
    4 CONTINUE
      END

