program Vibsymbol_0728
      implicit none
      real*8 :: Gamma, velocity, F, G, e_coeff, A, pi, phi, H, H_ks
      real*8, dimension(0:50, 0:50) :: horse_dustX, horse_dustY
      integer, dimension(0:50, 0:50) :: omega_ten
      integer, dimension(0:10, 0:50, 0:50) :: omega
      integer, dimension(0:2015) :: P
      real*8, dimension(0:2015) :: Abs_P
      integer :: Count, Num, I, J, KS_n

// initial
      print *, "H_n ? (n=1 to 5) : "
      read *, KS_n

      Gamma = 0.25
      pi = acos(-1.0d0)
      e_coeff = 1.0
      A = 9.8 / Gamma
      Count = 10

      do I = 0, 50
        do J = 0, 50
          horse_dustX(I, J) = I * 2.0 * pi / 50.0
          horse_dustY(I, J) = J * 2.0 * pi / 50.0 - horse_dustX(I, J)
          if(horse_dustX(I, J) < pi) then
            omega(0, I, J) = 0
          else
            omega(0, I, J) = 1
          end if
        end do
      end do

// main
      do Num = 1, Count
      do I = 0, 50
        do J = 0, 50
          phi = horse_dustX(I, J)
          horse_dustX(I, J) = F(horse_dustX(I, J), horse_dustY(I, J), pi)
          if(horse_dustX(I, J) < 0.0) then
            horse_dustX(I, J) = horse_dustX(I, J) + 2.0 * pi
          end if
          horse_dustY(I, J) = G(horse_dustX(I, J), horse_dustY(I, J), phi, e_coeff, A, pi)
          
          if(horse_dustX(I, J) < pi) then
            omega(Num, I, J) = 0
          else
            omega(Num, I, J) = 1
          end if

        end do
      end do
      end do

// 2 to 10
      do I = 0, 50
        do J = 0, 50
          do Num = 0, 10
            if(Num == 5) then
              
            else if(Num == 0 .OR. Num == 10) then
              if(KS_n <= 4) then
                
              else
                omega_ten(I, J) = omega_ten(I, J) + omega(Num, I, J) * (2 ** Num)
              end if
            else if(Num == 1 .OR. Num == 9) then
              if(KS_n <= 3) then
                
              else
                omega_ten(I, J) = omega_ten(I, J) + omega(Num, I, J) * (2 ** Num)
              end if
            else if(Num == 2 .OR. Num == 8) then
              if(KS_n <= 2) then
                
              else
                omega_ten(I, J) = omega_ten(I, J) + omega(Num, I, J) * (2 ** Num)
              end if
            else if(Num == 3 .OR. Num == 7) then
              if(KS_n <= 1) then
                
              else
                omega_ten(I, J) = omega_ten(I, J) + omega(Num, I, J) * (2 ** Num)
              end if
            else 
              omega_ten(I, J) = omega_ten(I, J) + omega(Num, I, J) * (2 ** Num)
            end if
          end do
        end do
      end do

// calculate P for I,J
      do I = 0, 50
        do J = 0, 50
          do Num = 0, 2015
            if(omega_ten(I, J) == Num) then
              P(Num) = P(Num) + 1
            end if
          end do
        end do
      end do

// average
      do Num = 0, 2015
        Abs_P(Num) = real(P(Num) / 2601)
        H = H + Abs_P(Num) * log(Abs_P(Num))
      end do
      H_ks = - H / (2.0 * real(KS_n))
      print *, "H_ks (for n = ", n, ") = ", H_ks

      end program Vibsymbol_0728

// function
      function F(dustX, dustY, pi)
        real*8 :: F, dustX, dustY, pi
        F = DMOD(dustX + dustY, 2.0 * pi)
      end function F

      function G(dustX, dustY, phi, e_coeff, A, pi)
        real*8 :: G, dustX, dustY, phi, e_coeff, A, pi
        G = e_coeff * (dustX - phi) + (1.0 + e_coeff) * A * COS(dustX) / 4.9
      end function G

--------------------------------------------
      program VibLyapnov_0729
      implicit none
      real*8 :: Gamma, velocity, F, G, e_coeff, A, pi, phi, phi2, Prelambda, Sum
      real*8, dimension(1:1000) :: Lyapnov
      real*8, dimension(0:50, 0:50) :: horse_dustX, horse_dustY
      real*8, dimension(0:999, 0:50, 0:50) :: U, V, norm_En, Mu_n, Lambda
      integer :: Count, Num, I, J

// initial
      open(8, STATUS="unknown", FILE="result_lyapnov.txt", POSITION="append")
      write(8, '(1X, A)') "#PROGRAM Ball_Laypnov"
      write(8, '(1X, A, F5.3)') "#gamma = ", 0.25

      Gamma = 0.25
      pi = acos(-1.0d0)
      e_coeff = 1.0
      A = 9.8 / Gamma
      Count = 1000

      do I = 0, 50
        do J = 0, 50
          horse_dustX(I, J) = I * 2.0 * pi / 50.0
          horse_dustY(I, J) = J * 2.0 * pi / 50.0 - horse_dustX(I, J)
          norm_En(0, I, J) = sqrt(horse_dustX(I, J) ** 2.0 + horse_dustY(I, J) ** 2.0)
        end do
      end do

// main
      do Num = 1, Count
      do I = 0, 50
        do J = 0, 50
          phi = horse_dustX(I, J)
          phi2 = horse_dustY(I, J)
          horse_dustX(I, J) = F(horse_dustX(I, J), horse_dustY(I, J), pi)
          if(horse_dustX(I, J) < 0.0) then
            horse_dustX(I, J) = horse_dustX(I, J) + 2.0 * pi
          end if
          horse_dustY(I, J) = G(horse_dustX(I, J), horse_dustY(I, J), phi, e_coeff, A, pi)
          
          U(Num, I, J) = horse_dustX(I, J) - phi
          V(Num, I, J) = horse_dustY(I, J) - phi2
          norm_En(Num, I, J) = sqrt(U(Num, I, J) ** 2.0 + V(Num, I, J) ** 2.0)
          Mu_n(Num, I, J) = norm_En(Num, I, J) / norm_En(Num -1, I, J)
          
        end do
      end do
      end do

// calculate lambda
      do I = 0, 50
        do J = 0, 50
          do Num = 1, Count
            Prelambda = Prelambda + log(Mu_n(Num, I, J))
            Lambda(Num, I, J) = Prelambda / real(Num)
          end do
        end do
      end do

//average
      do Num = 1, Count
        do I = 0, 50
          do J = 0, 50
            Sum = Sum + Lambda(Num, I, J)
          end do
        end do
        Lyapnov(Num) = Sum / 2601.0
        write(8, *) Lyapnov(Num)
        Sum = 0.0
      end do

      end program VibLyapnov_0728

// function
      function F(dustX, dustY, pi)
        real*8 :: F, dustX, dustY, pi
        F = DMOD(dustX + dustY, 2.0 * pi)
      end function F

      function G(dustX, dustY, phi, e_coeff, A, pi)
        real*8 :: G, dustX, dustY, phi, e_coeff, A, pi
        G = e_coeff * (dustX - phi) + (1.0 + e_coeff) * A * COS(dustX) / 4.9
      end function G