!===========================================================
!Discrete Logistic Map:
!Liapunov exponent from sum_i ln|f'(x_i)|
!Calculation for r in [rmin,rmax] with RSTEPS steps
! rsteps: values or r studied: r=rmin+(rmax-rmin)/RSTEPS
! ntrans: number of discarted iteration in order to discart
!         transient behaviour
! nsteps: number of terms in the sum
! xstart: value of initial x0 for every r
!===========================================================
program logistic_map
 implicit none
 integer , parameter :: dp     = 8
 real(dp), parameter :: rmin   = 2.5_dp
 real(dp), parameter :: rmax   = 4.0_dp
 real(dp), parameter :: xstart = 0.2_dp
 integer , parameter :: rsteps = 1000
 integer , parameter :: nsteps = 60000
 integer , parameter :: ntrans = 2000
 integer             :: i,ir
 real(dp)            :: r,x0,x1,sum,dr
 real(dp), parameter :: ONE    = 1.0_dp , TWO = 2.0_dp
 open(unit=33,file='lia.dat')
 dr    = (rmax-rmin)/(rsteps-1)
 do ir = 0,RSTEPS-1
  r    = rmin+ir*dr
  x0   = xstart
  do i = 1,ntrans
   x1  = r   * x0        * (ONE -     x0)
   x0  =       x1
  enddo
  sum  =       log(abs(r * (ONE - TWO*x0)))
! ----- Calculate:
  do i = 2,nsteps
   x1  = r   * x0        * (ONE -     x0)
   sum = sum + log(abs(r * (ONE - TWO*x1)))
   x0  = x1
  enddo
  write(33,*)r,sum/nsteps
 enddo !do ir=0,rsteps-1
 close(33)
end program logistic_map


!  ---------------------------------------------------------------------
!  Copyright by Konstantinos N. Anagnostopoulos (2004-2022)
!  Physics Dept., National Technical University,
!  konstant@mail.ntua.gr, www.physics.ntua.gr/konstant
!  
!  This program is free software: you can redistribute it and/or modify
!  it under the terms of the GNU General Public License as published by
!  the Free Software Foundation, version 3 of the License.
!  
!  This program is distributed in the hope that it will be useful, but
!  WITHOUT ANY WARRANTY; without even the implied warranty of
!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
!  General Public License for more details.
!  
!  You should have received a copy of the GNU General Public Liense along
!  with this program.  If not, see <http://www.gnu.org/licenses/>.
!  -----------------------------------------------------------------------