!====================================================================

 module global
   integer, parameter :: kdp=kind(0.d0)
   real (kind=kdp) dt, dt2, dt6
   real (kind=kdp) a                                ! some parameter
   real (kind=kdp), parameter :: pi = 4.d0*atan(1.d0)
 end module

!====================================================================
!
program main
!    
!====================================================================
!
!  This code integrates a 2D nonlinear system of the kind
!   
!  \dot(x1) = F1(x1,x2)
!  \dot(x2) = F2(x1,x2)
!
!  The code uses a 4th order Runge-Kutta algorithm.
!
!  This example illustrates the case of a non-hyperbolic fixed
!  point, a center which is changed into a spiral by arbitrarily 
!  small nonlinear terms (proportional to the parameter a, whose 
!  sign controls the stability of the spiral). This example is 
!  example 6.3.2 in S. Strogatz's book "Nonlinear dynamics and chaos"
!
!  The code can be easily generalized to higher dimension and to 
!  different nonlinear systems.
! 
!  Copyright 3/2017. F. Piazza, University of Orléans,
!                    Francesco.Piazza@gmail.com 
!
!====================================================================
!
use global
integer is, N
real (kind=kdp)  time
real (kind=kdp), dimension(2) :: x    ! the state vector
!
!====================================================================
!  Parameters

 dt   = 0.01d0            ! integration step 
 dt2  = dt/2.d0           !
 dt6  = dt/6.d0           !
 time = 0.d0              !
 is   = 10                ! sampling interval
 N    = 10000             ! Number of integration steps 
 a    = -0.02             ! nonlinearity parameter                      
 open (2,file='traj.dat') ! Trajectory file
! 
!====================================================================
! Initial condition 
!
x = (/ 2.d0, 0.01d0 /)
!
!====================================================================
! Time integration
!
 do i=1, N
    call one_RK4_step(x)
    time = i * dt
    if ( mod(i,is) == 0) then 
        print*, time, "   ", x(1), x(2)
        write (2,100) time, x(1), x(2) 
       end if
 end do
!
!====================================================================
 100 format(20(F16.8,3x)) 
!====================================================================
!
contains
!
!====================================================================

    subroutine one_RK4_step(x)   ! advances the time map one timestep
    use global 
    real (kind=kdp), dimension(2) :: x
    real (kind=kdp), dimension(2) :: k1, k2, k3, k4

    call Fvec(k1,x)
    call Fvec(k2,x + k1 * dt2)
    call Fvec(k3,x + k2 * dt2)
    call Fvec(k4,x + k3 * dt)
    x = x + dt6*(k1 + 2.d0*k2 + 2.d0*k3 + k4)

    end subroutine	
!
!--------------------------------------------------------------------
!    
    subroutine Fvec(F,x)   
    use global 
    real (kind=kdp), dimension(2) :: x, F

    F(1) = -x(2) + a*x(1)*(x(1)*x(1) + x(2)*x(2))                         
    F(2) =  x(1) + a*x(2)*(x(1)*x(1) + x(2)*x(2)) 

    end subroutine
!
!====================================================================
!
end program  
