program LyapunovSpectrum

! -----------------------------------------------------
! Program calculates the Lyapunov spectrum of the given
! dynamical system ie its Lyapunov exponents calculated
! on the section of the systems' parameter space.
! -----------------------------------------------------

! ban Fortran's implicit initialization of the variables 
implicit none


! ----------------------------------------------------
! 0 ADOPTED SHORTHAND REFERENCES OF THE COMMON NOTIONS
! ----------------------------------------------------

! LCES := 'dims' Lyapunov characteristic exponents
! parameter := parameter of the system
! parameters := 'params' parameters of the system
! when fixed := when they (parameters) are fixed
! when changing := when they (parameters) are changing
! # := number of ..
! sim := simulation := integration of the continuous system
! 					or iteration of the discrete system
! Lyapunov spectrum := LCES calculated for many values of
!					   the parameters (outside this code
!					   it denotes the set of all LCES)
! sub or SUB := subroutine returning scalar or vector value
! fun or FUN := function returning scalar or vector value


! -------------------------
! 1 CONSTANTS OF THE SYSTEM
! -------------------------

! LCES base --> (*)
real(8), parameter :: base = 2.0

! ----------------
! 1.1 Henon System
! ----------------

! number of dimensions of the system
!integer, parameter :: dims = 2
! number of parameters of the system
!integer, parameter :: params = 2

! --------------------------------------
! 1.2 Lorenz System & 1.3 Rossler System
! --------------------------------------

!integer, parameter :: dims = 3
!integer, parameter :: params = 3

! -------------------------------
! 1.4 Rabinovich-Fabrikant System
! -------------------------------

integer, parameter :: dims = 3
integer, parameter :: params = 2

! (*) "Please note that base 'e' and NOT base '2'
!	   is the one commonly used in the literature."


! -------------------------
! 2 VARIABLES OF THE SYSTEM
! -------------------------

! name of the system (up to 20 characters)
 character(len = 20) :: system_name_20
! "system is continuous in time" - .True. or .False.
logical :: continuous
! names of the parameters (up to 5 characters)
 character(len = 5), dimension(params) :: PARAMS_NAMES_5
! currently changing parameter + its shorthand notation
integer :: var_idx, i
! minimum values of the parameters when changing
real(8), dimension(params) :: VAR_MIN
! maximum values of the parameters when changing
real(8), dimension(params) :: VAR_MAX
! values of parameters when fixed
real(8), dimension(params) :: FIXED
! formats of the parameters' names (up to 4 characters) +
! formats of the parameters' values (up to 4 characters)
 character(len = 4), dimension(2 * params) :: PARAMS_FORMATS_4
! formats of the LCES' names (up to 4 characters) +
! formats of the LCES' values (up to 4 characters)
 character(len = 4), dimension(2 * dims) :: EXPONENTS_FORMATS_4
! # the LCES' calculations for the parameters when changing
integer, dimension(params) :: CALCS
! 'dims' initial conditions i.e. initial state of the system
real(8), dimension(dims) :: ICS
! time interval integrated during one simulation
! for the parameters when changing --> (**)
real(8), dimension(params) :: SIM_TIME
! # time steps taken during one simulation
! for the parameters when changing --> (**)
integer, dimension(params) :: STEPS
! ASSUMED # the simulations needed for transients of the system
! to vanish for the parameters when changing.
integer, dimension(params) :: TRANS_SIMS
! maximum allowed # the simulations
! for the parameters when changing
integer, dimension(params) :: MAX_SIMS
! relative error used during LCES convergence test
! for the parameters when changing
real(8), dimension(params) :: REL_ERROR
! absolute error used during LCES convergence test
! for the parameters when changing
real(8), dimension(params) :: ABS_ERROR
! maximum allowed absolute value of the Lyapunov exponent
! still counted as zero for the parameters when changing
real(8), dimension(params) :: ZERO_TOL

! (**) "Please initiate to -1.d0 for every parameter IF the system
!  	   is DISCRETE. It's just due to poor implementation reasons."


! ------------------
! 3 LCES CALCULATION
! ------------------

! ----------------
! 3.1 Henon system
! ----------------

! "Please remember to change # of dimensions of the system
! to 2 and # of parameters of the system to 2 (up in 1.1).
! We apologize for inconvenience."

!system_name_20 = "Henon"
! continuous = .False.
!PARAMS_NAMES_5 = (/ "a", "b" /)
!VAR_MIN = (/ 5.d-3, 1.d-3 /)
!VAR_MAX = (/ 1.4d0, 3.d-1 /)
!FIXED = (/ 1.4d0, 3.d-1 /)
!PARAMS_FORMATS_4 = (/ "A7  ", "A7  ", "F7.3", "F7.3" /)
!EXPONENTS_FORMATS_4 = (/ "A8  ", "A8  ", "F8.3", "F8.3" /)
!CALCS = (/ 281, 300 /)
!ICS = (/ 1, 0 /)
!SIM_TIM = (/ -1.d0, -1.d0 /) ! unused
!STEPS = (/ -1.d0, -1.d0 /)	  ! unused
!TRANS_SIMS = (/ 50, 50 /)
!MAX_SIMS = (/ 10**7, 10**7 /)
!REL_ERROR = (/ 1.d-7, 1.d-7 /)
!ABS_ERROR = (/ 1.d-7, 1.d-7 /)
!ZERO_TOL = (/ 1.d-3, 1.d-3 /)

! -----------------
! 3.2 Lorenz system
! -----------------

! "Please remember to change # of dimensions of the system
! to 3 and # of parameters of the system to 3 (up in 1.2).
! Also please change the function to be integrated to
! 'LORENZ' (down in 5.5).
! We apologize for inconvenience."

!system_name_20 = "Lorenz"
! continuous = .True.
!PARAMS_NAMES_5 = (/ "s", "r", "b" /)
!VAR_MIN = (/ 0, 0, 0 /)
!VAR_MAX = (/ 60, 60, 60 /)
!FIXED = (/ 1.6d1, 4.592d1, 4.d0 /)
!PARAMS_FORMATS_4 = (/ "A7  ", "A7  ", "A7  ", "F7.2", "F7.2", "F7.2" /)
!EXPONENTS_FORMATS_4 = (/ "A8  ", "A8  ", "A8  ", "F8.3", "F8.3", "F8.3" /)
!CALCS = (/ 601, 601, 601 /)
!ICS = (/ 10, 1, 0 /)
!SIM_TIME = (/ 5.d-1, 5.d-1, 5.d-1 /)
!STEPS = (/ 50, 50, 50 /)
!TRANS_SIMS = (/ 50, 50, 50 /)
!MAX_SIMS = (/ 55000, 55000, 55000 /)
!REL_ERROR = (/ 1.d-6, 1.d-6, 1.d-6 /)
!ABS_ERROR = (/ 1.d-6, 1.d-6, 1.d-6 /)
!ZERO_TOL = (/ 1.49d-3, 1.49d-3, 1.49d-3 /)

! ------------------
! 3.3 Rossler system
! ------------------

! "Please remember to change # of dimensions of the system
! to 3 and # of parameters of the system to 3 (up in 1.3).
! Also please change the function to be integrated to
! 'ROSSLER' (down in 5.5)."
! We apologize for inconvenience."

!system_name_20 = "Rossler"
! continuous = .True.
!PARAMS_NAMES_5 = (/ "a", "b", "c" /)
!VAR_MIN = (/ 0.d0, 1.d-1, 4.d-1 /)
!VAR_MAX = (/ 3.75d-1, 2.d1, 3.7d1 /)
!FIXED = (/ 1.5d-1, 2.d-1, 10.d0 /)
!PARAMS_FORMATS_4 = (/ "A7  ", "A7  ", "A7  ", "F7.2", "F7.2", "F7.2" /)
!EXPONENTS_FORMATS_4 = (/ "A8  ", "A8  ", "A8  ", "F8.3", "F8.3", "F8.3" /)
!CALCS = (/ 376, 501, 367 /)
!ICS = (/ 10, 1, 0 /)
!SIM_TIME = (/ 2.d0, 5.d-1, 5.d-1 /)
!STEPS = (/ 200, 50, 50 /)
!TRANS_SIMS = (/ 20, 50, 50 /)
!MAX_SIMS = (/ 10**5, 6*10**4, 6*10**4 /)
!REL_ERROR = (/ 1.d-6, 1.d-5, 1.d-5 /)
!ABS_ERROR = (/ 1.d-6, 1.d-5, 1.d-5 /)
!ZERO_TOL = (/ 1.1d-3, 5.d-3, 5.d-3 /)

! -------------------------------
! 3.4 Rabinovich-Fabrikant System
! -------------------------------

! "Please remember to change # of dimensions of the system
! to 3 and # of parameters of the system to 2 (up in 1.4).
! Also please change the function to be integrated to
! 'RABFAB' (down in 5.5)."
! We apologize for inconvenience."

system_name_20 = "RabFab"
 continuous = .True.
PARAMS_NAMES_5 = (/ "a", "g" /)
VAR_MIN = (/ 1.39, 0.0 /)
VAR_MAX = (/ 1.61, 1.375 /)
FIXED = (/ 1.59, 1.26 /)
PARAMS_FORMATS_4 = (/ "A7  ", "A7  ", "F7.3", "F7.3" /)
EXPONENTS_FORMATS_4 = (/ "A8  ", "A8  ", "A8  ", "F8.3", "F8.3", "F8.3" /)
CALCS = (/ 221, 276 /)
ICS = (/ -1.0, 1.0, 0.5 /)
SIM_TIME = (/ 0.5, 0.5 /)
STEPS = (/ 50, 50 /)
TRANS_SIMS = (/ 50, 50 /)
MAX_SIMS = (/ 5*10**4, 5*10**4 /)
REL_ERROR = (/ 1.d-6, 1.d-6 /)
ABS_ERROR = (/ 1.d-6, 1.d-6 /)
ZERO_TOL = (/ 1.d-3, 1.d-3 /)


! -------------------------------
! 4 LYAPUNOV SPECTRUM CALCULATION
! -------------------------------

do var_idx = 1, params

	i = var_idx

	call CALC_LYAPUNOV_SPECTRUM(system_name_20, dims, continuous, params, &
								PARAMS_NAMES_5, var_idx, VAR_MIN(i), VAR_MAX(i), &
							    FIXED, PARAMS_FORMATS_4, EXPONENTS_FORMATS_4, &
							    CALCS(i), ICS, SIM_TIME(i), STEPS(i), TRANS_SIMS(i), &
							    MAX_SIMS(i), REL_ERROR(i), ABS_ERROR(i), base, ZERO_TOL(i))

end do


! --------------------------------------------------
! 5 SUBROUTINES AND FUNCTIONS DEFINITIONS AND BODIES
! --------------------------------------------------

contains

! --------------------------
! 5.1 'CREATE_EVOLUTION' SUB
! --------------------------

subroutine CREATE_EVOLUTION(dims, ICS, EVOLUTION)

	! creates 'EVOLUTION' array containing 'ICS' state
	! and 'dims' ortonormal axes of the initial state ball

	implicit none

	integer, intent(in) :: dims
	real(8), dimension(dims), intent(in) :: ICS
	real(8), dimension(dims * (1 + dims)), intent(out) :: EVOLUTION

	integer :: idx

	EVOLUTION = 0.0
	EVOLUTION(:dims) = ICS

	do idx = 1, dims

		EVOLUTION(idx * (1 + dims)) = 1.0

	end do

end subroutine CREATE_EVOLUTION

! --------------
! 5.2 'norm' fun
! --------------

function norm(dims, VECTOR)

	! calculates the norm of array 'VECTOR'
	! containing 'dims' elements

	implicit none

	integer, intent(in) :: dims
	real(8), dimension(dims), intent(in) :: VECTOR
	real(8) :: norm

	intrinsic sqrt, dot_product
	
	norm = sqrt(dot_product(VECTOR, VECTOR))

end function norm

! --------------------
! 5.3 'PROJECTION' FUN
! --------------------

function PROJECTION(dims, VECTOR, VERSOR)

	! determines the projection of array 'VECTOR'
	! on array 'VERSOR' both containing 'dims' elements

	implicit none

	integer, intent(in) :: dims
	real(8), dimension(dims), intent(in) :: VECTOR
	real(8), dimension(dims), intent(in) :: VERSOR
	real(8), dimension(dims) :: PROJECTION

	intrinsic dot_product

	PROJECTION = dot_product(VECTOR, VERSOR) * VERSOR
	
end function PROJECTION

! ------------------------------
! 5.4 'GS_ORTONORMALIZATION' SUB
! ------------------------------

subroutine GS_ORTONORMALIZATION(dims, EVOLUTION, NORMS)

	! performs Gram-Schmidt ortonormalization process on
	! 'dims' axes contained in the 'EVOLUTION' array +
	! returns 'NORMS' array containing norms of axes after
	! their ortogonalization BUT before their normalization

	implicit none

	integer, intent(in) :: dims
	real(8), dimension(dims * (1 + dims)) :: EVOLUTION
	real(8), dimension(dims), intent(out) :: NORMS

	integer :: idx, jdx
	real(8), dimension(dims) :: VECTOR, VERSOR

	do idx = 1, dims

		VECTOR = EVOLUTION(1 + idx * dims : (1 + idx) * dims)

		do jdx = 1, idx - 1

			VERSOR = EVOLUTION(1 + (idx - jdx) * dims : (1 + idx - jdx) * dims)
			VECTOR = VECTOR - PROJECTION(dims, VECTOR, VERSOR)

		end do

		NORMS(idx) = norm(dims, VECTOR)
		VECTOR = VECTOR / NORMS(idx)
		EVOLUTION(1 + idx * dims : (1 + idx) * dims) = VECTOR

	end do

end subroutine GS_ORTONORMALIZATION

! --------------------
! 5.5 'RK4_SOLVER' SUB
! --------------------

subroutine RK4_SOLVER(dims, EVOLUTION, sim_time, steps, PARAMETERS)

	! Integrates chosen 'dims' dimensional system depending on
	! parameters whose values are given in 'PARAMETERS.'
	! Unfortunately, despite substantial efforts to change it,
	! on has to specify manualy system to be integrated in the
	! marked body of the subroutine.

	implicit none

	integer, intent(in) :: dims
	real(8), dimension(dims * (1 + dims)) :: EVOLUTION
	real(8), intent(in) :: sim_time
	integer, intent(in) :: steps
	real(8), dimension(10), intent(in) :: PARAMETERS

	real(8) :: time_step
	integer :: step
	real(8), dimension(dims * (1 + dims)) :: RK4_K1, RK4_K2
	real(8), dimension(dims * (1 + dims)) :: RK4_K3, RK4_K4
	real(8), dimension(dims * (1 + dims)) :: RK4_K

	intrinsic dble

	time_step = sim_time / dble(steps)

	do step = 1, steps

		! "Please choose a continuous system to be integrated from a list below
		! or paste analogous 4 lines with the name of the system to be integrated,
		! if you have written a function calculating values of the system yourself."

		!RK4_K1 = LORENZ(dims, EVOLUTION,			     PARAMETERS) * time_step
		!RK4_K2 = LORENZ(dims, EVOLUTION + RK4_K1 / 2.d0, PARAMETERS) * time_step
		!RK4_K3 = LORENZ(dims, EVOLUTION + RK4_K2 / 2.d0, PARAMETERS) * time_step
		!RK4_K4 = LORENZ(dims, EVOLUTION + RK4_K3,	     PARAMETERS) * time_step

		!RK4_K1 = ROSSLER(dims, EVOLUTION,			       PARAMETERS) * time_step
		!RK4_K2 = ROSSLER(dims, EVOLUTION + RK4_K1 / 2.d0, PARAMETERS) * time_step
		!RK4_K3 = ROSSLER(dims, EVOLUTION + RK4_K2 / 2.d0, PARAMETERS) * time_step
		!RK4_K4 = ROSSLER(dims, EVOLUTION + RK4_K3,	       PARAMETERS) * time_step

		RK4_K1 = RABFAB(dims, EVOLUTION,			      PARAMETERS) * time_step
		RK4_K2 = RABFAB(dims, EVOLUTION + RK4_K1 / 2.d0, PARAMETERS) * time_step
		RK4_K3 = RABFAB(dims, EVOLUTION + RK4_K2 / 2.d0, PARAMETERS) * time_step
		RK4_K4 = RABFAB(dims, EVOLUTION + RK4_K3,	      PARAMETERS) * time_step

		RK4_K = (RK4_K1 + 2.d0 * RK4_K2 + 2.d0 * RK4_K3 + RK4_K4) / 6.d0

		EVOLUTION = EVOLUTION + RK4_K

	end do

end subroutine RK4_SOLVER

! ---------------
! 5.6 'HENON' FUN
! ---------------

function HENON(dims, EVOLUTION, PARAMETERS)

	! Calculates the next value of the Henon map for the point
	! evolving in the state space + calculates 'dims' next values
	! of the Henon map Jacobian for the 'dims' axes of the
	! sphere of initial states.

	implicit none

	integer, intent(in) :: dims
	real(8), dimension(dims * (1 + dims)), intent(in) :: EVOLUTION
	real(8), dimension(10), intent(in) :: PARAMETERS
	real(8), dimension(dims * (1 + dims)) :: HENON

	real(8) :: a, b
	real(8) :: x, y
	real(8), dimension(dims, dims) :: JACOBIAN
	integer :: idx
	real(8), dimension(dims) :: VECTOR

	intrinsic TRANSPOSE, RESHAPE, shape, MATMUL

	HENON = EVOLUTION

	a = PARAMETERS(1)
	b = PARAMETERS(2)

	x = HENON(1)
	y = HENON(2)

	HENON(1) = 1.d0 - a * x * x + y
	HENON(2) = b * x

	JACOBIAN = TRANSPOSE(RESHAPE((/ &
			   -2.d0 * a * x, 1.d0, &
			   			   b, 0.d0  &
			   /), shape(JACOBIAN)))

	do idx = 1, dims

		VECTOR = HENON(1 + idx * dims : (1 + idx) * dims)
		VECTOR = MATMUL(JACOBIAN, VECTOR)
		HENON(1 + idx * dims : (1 + idx) * dims) = VECTOR

	end do

end function HENON

! ----------------
! 5.7 'LORENZ' FUN
! ----------------

function LORENZ(dims, EVOLUTION, PARAMETERS)

	! Calculates the next value of the Lorenz system for the point
	! evolving in the state space + calculates 'dims' next values
	! of the Lorenz system Jacobian for the 'dims' axes of the
	! sphere of initial states.

	implicit none

	integer, intent(in) :: dims
	real(8), dimension(dims * (1 + dims)), intent(in) :: EVOLUTION
	real(8), dimension(10), intent(in) :: PARAMETERS
	real(8), dimension(dims * (1 + dims)) :: LORENZ

	real(8) :: s, r, b
	real(8) :: x, y, z
	real(8), dimension(dims, dims) :: JACOBIAN
	integer :: idx
	real(8), dimension(dims) :: VECTOR

	intrinsic TRANSPOSE, RESHAPE, shape, MATMUL

	LORENZ = EVOLUTION

	s = PARAMETERS(1)
	r = PARAMETERS(2)
	b = PARAMETERS(3)

	x = LORENZ(1)
	y = LORENZ(2)
	z = LORENZ(3)

	LORENZ(1) = s * (y - x)
	LORENZ(2) = (r - z) * x - y
	LORENZ(3) = x * y - b * z

	JACOBIAN = TRANSPOSE(RESHAPE((/ &
			      -s,     s, 0.d0, &
			   r - z, -1.d0,   -x, &
			   	   y,     x,   -b  &
			   /), shape(JACOBIAN)))

	do idx = 1, dims

		VECTOR = LORENZ(1 + idx * dims : (1 + idx) * dims)
		VECTOR = MATMUL(JACOBIAN, VECTOR)
		LORENZ(1 + idx * dims : (1 + idx) * dims) = VECTOR

	end do

end function LORENZ

! -----------------
! 5.8 'ROSSLER' FUN
! -----------------

function ROSSLER(dims, EVOLUTION, PARAMETERS)

	! Calculates the next value of the Rossler system for the point
	! evolving in the state space + calculates 'dims' next values
	! of the Rossler system Jacobian for the 'dims' axes of the
	! sphere of initial states.

	implicit none

	integer, intent(in) :: dims
	real(8), dimension(dims * (1 + dims)), intent(in) :: EVOLUTION
	real(8), dimension(10), intent(in) :: PARAMETERS
	real(8), dimension(dims * (1 + dims)) :: ROSSLER

	real(8) :: a, b, c
	real(8) :: x, y, z
	real(8), dimension(dims, dims) :: JACOBIAN
	integer :: idx
	real(8), dimension(dims) :: VECTOR

	intrinsic TRANSPOSE, RESHAPE, shape, MATMUL

	ROSSLER = EVOLUTION

	a = PARAMETERS(1)
	b = PARAMETERS(2)
	c = PARAMETERS(3)

	x = ROSSLER(1)
	y = ROSSLER(2)
	z = ROSSLER(3)

	ROSSLER(1) = -(y + z)
	ROSSLER(2) = x + a * y
	ROSSLER(3) = b + z * (x - c)

	JACOBIAN = TRANSPOSE(RESHAPE((/ &
			   0.d0, -1.d0, -1.d0, &
			   1.d0,     a,  0.d0, &
			   	  z,  0.d0, x - c  &
			   /), shape(JACOBIAN)))

	do idx = 1, dims

		VECTOR = ROSSLER(1 + idx * dims : (1 + idx) * dims)
		VECTOR = MATMUL(JACOBIAN, VECTOR)
		ROSSLER(1 + idx * dims : (1 + idx) * dims) = VECTOR

	end do

end function ROSSLER

! ----------------
! 5.9 'RABFAB' FUN
! ----------------

function RABFAB(dims, EVOLUTION, PARAMETERS)

	! Calculates the next value of the Rabinovich-Fabrikant system
	! for the point evolving in the state space + calculates 'dims'
	! next values of the Rabinovich-Fabrikant system Jacobian for
	! the 'dims' axes of the sphere of initial states.

	implicit none

	integer, intent(in) :: dims
	real(8), dimension(dims * (1 + dims)), intent(in) :: EVOLUTION
	real(8), dimension(10), intent(in) :: PARAMETERS
	real(8), dimension(dims * (1 + dims)) :: RABFAB

	real(8) :: a, g
	real(8) :: x, y, z
	real(8), dimension(dims, dims) :: JACOBIAN
	integer :: idx
	real(8), dimension(dims) :: VECTOR

	intrinsic TRANSPOSE, RESHAPE, shape, MATMUL

	RABFAB = EVOLUTION

	a = PARAMETERS(1)
	g = PARAMETERS(2)

	x = RABFAB(1)
	y = RABFAB(2)
	z = RABFAB(3)

	RABFAB(1) = y * (z - 1.d0 + x**2) + g * x
	RABFAB(2) = x * (3.d0 * z + 1.d0 - x**2) + g * y
	RABFAB(3) = -2.d0 * z * (a + x * y)

	JACOBIAN = TRANSPOSE(RESHAPE((/ &
	    					2.d0*x*y + g, z - 1.d0 + x**2,               y, &
			   3.d0*z + 1.d0 - 3.d0*x**2,               g,          3.d0*x, &
		                       -2.d0*y*z,       -2.d0*x*z, -2.d0*(a + x*y)  &
			   /), shape(JACOBIAN)))

	do idx = 1, dims

		VECTOR = RABFAB(1 + idx * dims : (1 + idx) * dims)
		VECTOR = MATMUL(JACOBIAN, VECTOR)
		RABFAB(1 + idx * dims : (1 + idx) * dims) = VECTOR

	end do

end function RABFAB

! ----------------------------------
! 5.10 'CALC_LYAPUNOV_EXPONENTS' SUB
! ----------------------------------

subroutine CALC_LYAPUNOV_EXPONENTS(dims, continuous, ICS, sim_time, steps, &
								   trans_sims, max_sims, rel_error, abs_error, &
								   base, PARAMETERS, sims, LYAPUNOV_EXPONENTS)

	! Calculates Lyapunov exponents in the 'base' base for the 'dims' dimensional
	! continuous ('continuous' = .True.) or discrete ('continuous' = .False.)
	! system depending on parameters with 'PARAMETERS' values and starting from
	! 'ICS' iniial state. Subroutine performs up to 'max_sims' simulations
	! each of 'sim_time' time length devided into 'steps' time steps unless
	! the relative-absolute test with 'rel_error' relative error value and
	! 'abs_error' absolute error value is satisfied: in that case calculation
	! is terminated and Lyapunov exponents are returned. Before actual simulation
	! 'trans_sims' simulations are performed, in order to ged rid of the transients.

	implicit none

	integer, intent(in) :: dims
	logical, intent(in) :: continuous
	real(8), dimension(dims), intent(in) :: ICS
	real(8), intent(in) :: sim_time
	integer, intent(in) :: steps
	integer, intent(in) :: trans_sims
	integer, intent(in) :: max_sims
	real(8), intent(in) :: rel_error
	real(8), intent(in) :: abs_error
	real(8), intent(in) :: base
	real(8), dimension(10), intent(in) :: PARAMETERS

	real(8), dimension(dims * (1 + dims)) :: EVOLUTION
	integer :: trans_sim
	real(8), dimension(dims) :: LYAPUNOV_SUMS
	real(8), dimension(dims), intent(out) :: LYAPUNOV_EXPONENTS
	integer, intent(out) :: sims
	logical :: converged
	real(8), dimension(dims) :: PREV_LYAPUNOV_EXPONENTS
	real(8), dimension(dims) :: NORMS

	intrinsic log

	call CREATE_EVOLUTION(dims, ICS, EVOLUTION)

	do trans_sim = 1, trans_sims

		if (continuous) then

			call RK4_SOLVER(dims, EVOLUTION, sim_time, steps, PARAMETERS)

		else

			EVOLUTION = HENON(dims, EVOLUTION, PARAMETERS)

		end if

		call GS_ORTONORMALIZATION(dims, EVOLUTION, NORMS)

	end do

	LYAPUNOV_SUMS = 0.0
	LYAPUNOV_EXPONENTS = 0.0
	sims = 0
	converged = .false.

	do while (.NOT. converged) 

		sims = sims + 1

		if (sims >= max_sims) then

			write(*,*) "No convergence."
			exit

		end if

		PREV_LYAPUNOV_EXPONENTS = LYAPUNOV_EXPONENTS

		if (continuous) then

			call RK4_SOLVER(dims, EVOLUTION, sim_time, steps, PARAMETERS)

		else

			EVOLUTION = HENON(dims, EVOLUTION, PARAMETERS)

		end if

		call GS_ORTONORMALIZATION(dims, EVOLUTION, NORMS)

		LYAPUNOV_SUMS = LYAPUNOV_SUMS + LOG(NORMS) / log(base)

		if (continuous) then

			LYAPUNOV_EXPONENTS = LYAPUNOV_SUMS / (sims * sim_time)

		else

			LYAPUNOV_EXPONENTS = LYAPUNOV_SUMS / sims

		end if

		converged = norm(dims, PREV_LYAPUNOV_EXPONENTS - LYAPUNOV_EXPONENTS) &
				  < rel_error * norm(dims, LYAPUNOV_EXPONENTS) + abs_error

	end do

end subroutine CALC_LYAPUNOV_EXPONENTS

! ----------------------
! 5.11 'INSERT_SORT' SUB
! ----------------------

subroutine INSERT_SORT(dims, ARRAY)

	! Sorts the given array descendantly.

	implicit none

	integer, intent(in) :: dims
	real(8), dimension(dims) :: ARRAY

	integer :: idx, jdx
	real(8) :: key

	do idx = 2, dims

		key = ARRAY(idx)
		jdx = idx - 1

		do while (jdx > 0 .and. key > ARRAY(jdx))

			ARRAY(jdx + 1) = ARRAY(jdx)
			ARRAY(jdx) = key
			jdx = jdx - 1

		end do

	end do

end subroutine INSERT_SORT

! -----------------------------
! 5.12 'lyapunov_dimension' fun
! -----------------------------

function lyapunov_dimension(dims, LYAPUNOV_EXPONENTS, zero_tol)

	! Calculates Lyapunov dimension given the array
	! 'LYAPUNOV_EXPONENTS' containing Lyapunov exponents.
	! Since obtainment of zero-exponenent is unreachable
	! numerically, values whose absolute value is less
	! than 'zero_tol' value will be counted as zero-exponents.

	implicit none

	integer, intent(in) :: dims
	real(8), dimension(dims), intent(in) :: LYAPUNOV_EXPONENTS
	real(8), intent(in) :: zero_tol
	real(8) :: lyapunov_dimension

	real(8), dimension(dims) :: LYAPUNOV_EXPS_COPY
	real(8) :: exponents_sum
	real(8) :: prev_exponents_sum
	integer :: idx

	intrinsic abs

	LYAPUNOV_EXPS_COPY = LYAPUNOV_EXPONENTS

	do idx = 1, dims
	
		if (abs(LYAPUNOV_EXPS_COPY(idx)) <= zero_tol) then

			LYAPUNOV_EXPS_COPY(idx) = 0.d0

		end if

	end do
	
	call INSERT_SORT(dims, LYAPUNOV_EXPS_COPY)

	exponents_sum = LYAPUNOV_EXPS_COPY(1)
	idx = 1

	do while (exponents_sum >= 0 .and. idx < dims)

		idx = idx + 1
		prev_exponents_sum = exponents_sum
		exponents_sum = exponents_sum + LYAPUNOV_EXPS_COPY(idx)
		lyapunov_dimension = (idx - 1) + prev_exponents_sum &
						      / abs(LYAPUNOV_EXPS_COPY(idx))

	end do

	if (abs(lyapunov_dimension) < zero_tol) then

		lyapunov_dimension = 0.d0

	end if 

end function lyapunov_dimension

! -------------------------
! 5.13 'attractor_type' fun
! -------------------------

function attractor_type(lyapunov_dimension)

	! Returns human-readable name of the attractor type
	! corresponding to the Lyapunov dimension value.

	implicit none

	real(8), intent(in) :: lyapunov_dimension
	character(len = 10) :: attractor_type

	if (lyapunov_dimension == 0.0) then

		attractor_type = "fixed"

	else if (lyapunov_dimension == 1.0) then

		attractor_type = "periodic"

	else if (lyapunov_dimension == 2.0) then

		attractor_type = "2-periodic"

	else

		attractor_type = "chaotic"

	end if

end function attractor_type

! ---------------------------------
! 5.14 'CALC_LYAPUNOV_SPECTRUM' SUB
! ---------------------------------

subroutine CALC_LYAPUNOV_SPECTRUM(system_name_20, dims, continuous, params, &
								  PARAMS_NAMES_5, var_idx, var_min, var_max, &
								  FIXED, PARAMS_FORMATS_4, EXPONENTS_FORMATS_4, &
								  calcs, ICS, sim_time, steps, trans_sims, &
								  max_sims, rel_error, abs_error, base, zero_tol)

	! Calculates Lyapunov spectrum of the system calling 'CALC_LYAPUNOV_EXPONENTS'
	! numerous times. Each parameter is scanned from 'var_min' value to the
	! 'var_max' value using 'calcs' steps. Beside arguments explained in the
	! 'CALC_LYAPUNOV_EXPONENTS' subroutine description one shoul also specify
	! the name 'system_name_20' of the system and 'PARAMS_NAMES_5' parameters,
	! formats of the parameters 'PARAMS_FORMATS_4' and exponents 'EXPONENTS_FORMATS_4'
	! (first specifying names' formats of all parameters and then specifying values'
	! formats of all parameters) and values 'FIXED' of the parameters when they are
	! fixed.

	implicit none

	character(len = 20), intent(in) :: system_name_20
	integer, intent(in) :: dims
	logical, intent(in) :: continuous
	integer, intent(in) :: params
	character(len = 5), dimension(params), intent(in) :: PARAMS_NAMES_5
	integer, intent(in) :: var_idx
	real(8), intent(in) :: var_min, var_max
	real(8), dimension(params), intent(in) :: FIXED
	character(len = 4), dimension(2 * params), intent(in) :: PARAMS_FORMATS_4
	character(len = 4), dimension(2 * dims), intent(in) :: EXPONENTS_FORMATS_4
	integer, intent(in) :: calcs
	real(8), dimension(dims), intent(in) :: ICS
	real(8), intent(in) :: sim_time
	integer, intent(in) :: steps
	integer, intent(in) :: trans_sims
	integer, intent(in) :: max_sims
	real(8), intent(in) :: rel_error
	real(8), intent(in) :: abs_error
	real(8), intent(in) :: base
	real(8), intent(in) :: zero_tol

	integer :: open_c
	character(len = 2) :: char_idx
	real(8) :: total_time
	real(8), dimension(10) :: PARAMETERS
	integer :: idx
	integer :: calc
	real(8) :: start, finish
	real(8) :: var_step
	real(8) :: var_value
	integer :: sims
	real(8), dimension(dims) :: LYAPUNOV_EXPONENTS
	real(8) :: lyapunov_dim
	character(len = 10) :: attr_type
	
	intrinsic trim, adjustl, cpu_time, int

	open(unit = 2, file = "./LCE_" // trim(system_name_20) // "_" // &
		 trim(PARAMS_NAMES_5(var_idx)) // ".dat", action = "write", &
		 status = "replace", access = "append", iostat = open_c)

	print *,
	write(*,*) "LCE_" // trim(system_name_20) // "_" // &
			   trim(PARAMS_NAMES_5(var_idx))
	print *,

	write(*,"(A6)", advance = "no") "calc"

	do idx = 1, params
	
		write(*, "(" // trim(PARAMS_FORMATS_4(idx)) // ")", &
			  advance = "no") trim(PARAMS_NAMES_5(idx))
	
	end do

	write(*, 5, advance = "no") "sim_t", "tot_t", "sims"
	5 format(A7, A6, A8)

	do idx = 1, dims

		write(char_idx, "(I2)") idx
		write(*, "(" // trim(EXPONENTS_FORMATS_4(idx)) // ")", &
			  advance = "no") "L" // trim(adjustl(char_idx))

	end do

	write(*, 6) "I_dim", "attr_type"
	6 format(A7, A12)

	write(2,"(A6)", advance = "no") "calc"

	do idx = 1, params
	
		write(2, "("// trim(PARAMS_FORMATS_4(idx)) // ")", &
			  advance = "no") trim(PARAMS_NAMES_5(idx))

	end do

	write(2, 5, advance = "no") "sim_t", "tot_t", "sims"

	do idx = 1, dims
	
		write(char_idx, "(I2)") idx
		write(2, "("// trim(EXPONENTS_FORMATS_4(idx)) // ")", &
			  advance = "no") "L" // trim(adjustl(char_idx))
	
	end do

	write(2, 6) "I_dim", "attr_type"

	total_time = 0.0
	PARAMETERS = 0.0

	do idx = 1, params

		if (idx /= var_idx) then

			PARAMETERS(idx) = FIXED(idx)

		end if

	end do

	var_step = (var_max - var_min) / (calcs - 1)

	do calc = 1, calcs

		call cpu_time(start)

		var_value = var_min + (calc - 1) * var_step
		PARAMETERS(var_idx) = var_value

		call CALC_LYAPUNOV_EXPONENTS(dims, continuous, ICS, sim_time, steps, &
									 trans_sims, max_sims, rel_error, abs_error, &
									 base, PARAMETERS, sims, LYAPUNOV_EXPONENTS)

		call cpu_time(finish)

		total_time = total_time + (finish - start)

		lyapunov_dim = lyapunov_dimension(dims, LYAPUNOV_EXPONENTS, zero_tol)
		attr_type = attractor_type(lyapunov_dim)

		write(*,"(I6)", advance = "no") calc

		do idx = 1, params

			write(*, "("// trim(PARAMS_FORMATS_4(idx + params)) // ")", &
				  advance = "no") PARAMETERS(idx)

		end do

		write(*, 7, advance = "no") finish - start, int(total_time), sims
		7 format(F7.2, I6, I8)

		do idx = 1, dims

			write(*, "("// trim(EXPONENTS_FORMATS_4(idx + dims)) // ")", &
				  advance = "no") LYAPUNOV_EXPONENTS(idx)

		end do

		write(*, 8) lyapunov_dim, attr_type
		8 format(F7.3, A12)

		write(2,"(I6)", advance = "no") calc

		do idx = 1, params

			write(2, "("// trim(PARAMS_FORMATS_4(idx + params)) // ")", &
				  advance = "no") PARAMETERS(idx)

		end do

		write(2, 7, advance = "no") finish - start, int(total_time), sims

		do idx = 1, dims

			write(2, "("// trim(EXPONENTS_FORMATS_4(idx + dims)) // ")", &
				  advance = "no") LYAPUNOV_EXPONENTS(idx)

		end do

		write(2, 8) lyapunov_dim, attr_type

	end do

	print *,

	close(2)

end subroutine CALC_LYAPUNOV_SPECTRUM

end program LyapunovSpectrum
