!============================================================================
!                       General Documentation Section
! Name:
!   LAMBDA_LOG_DISP
!
! Purpose:
!   Calculate the time-dependent exponent (lambda) and logarithmic dis-
!   placement curves, as defined by Gao (1997).
!
! Calling Sequence:
!   call LAMBDA_LOG_DISP( x0, k, m, L, ie_value  &
!                       , tan_motion_lim  &
!                       , normalize_input  &
!                       , out_lambda, out_log_disp  &
!                       , out_npts, out_shells )
!
!   This subroutine must have an explicit interface where called.
!
! Example:
!   None.
!
! File I/O:
!   None.
!
! Entry and Exit States:
!   N/A.
!
! RCS Information:
!   $Id: lambda_log_disp.f90,v 1.1 2003/04/04 04:39:27 jlin Exp jlin $
!
! Revision History:
! - 1993:  Original version by Jianbo Gao, UCLA Dept. of Electrical Engin-
!   eering.
! - 22 May 2002:  Fortran 90 commenting and formatting by Johnny Lin,
!   CIRES/University of Colorado, Boulder.  Email:  air_jlin@yahoo.com.
!   Passed a few minimal tests.
!
! Notes:
! - References:
!   * Gao, J. B. (1997), "Recognizing randomness in a time series," 
!     Physica D, Vol. 106, pp. 49-56.
!   * Gao J. B., and Zheng Z. M. (1993), "Local exponential divergence 
!     plot and optimal embedding of a chaotic time-series," Phys. Lett. A,
!     Vol. 181, No. 2, pp. 153-158.
!   * Gao, J. and Z. Zheng (1994), "Direct dynamical test for deterministic
!     chaos and optimal embedding of a chaotic time series," Phys. Rev. E,
!     Vol. 49, No. 5, pp. 3807-3814.
!   * Theiler, J. (1990), "Estimating fractal dimension," J. Opt. Soc. Am.
!     A, Vol. 7, No. 6, pp. 1055-1073.
! - See http://www.johnny-lin.com/fort_refs/lambda_log_disp.notes.html for
!   important notes regarding parameter choice.
! - In computations, normalized time series (i.e. in interval [0,1]) is
!   assumed.  If the input is not normalized, set parameter normalize_input
!   to 1 to make timeseries normalized before beginning computations.
! - Prescribed rulers (squared):  Given an array rl(9), of base 2 and 
!   exponent (-i) {i=ie,ie+1,...}, to define shells.  Note that to properly 
!   define the ensemble averages at least several hundreds pairs of points 
!   are needed.  See out_npts for the numbers used in averaging.
! - The number of points actually used in calculation (local variable nn)
!   is not the same as the number of points in the input timeseries,
!   because one does not necessarily need all those points to adequately 
!   span the phase space.  Thus, "nn" is set to a value based on input
!   timeseries length, embedding dimension, etc. or 8000, whichever is 
!   smaller.  If this is unacceptable, goto the "Main Parameters Setting" 
!   sub-section and change "nn" accordingly.
! - Note:  This routine does not use "implicit none," because when it
!   was first written, there were a large number of variables that were 
!   implicitly declared.
! - No subprograms required by this routine besides those included in
!   the Fortran 90 standard.
!
! Copyright (c) 1993-2002 by Jianbo Gao and Johnny Lin.  For  licensing 
! and contact information see http://www.johnny-lin.com/flib.html.
!----------------------------------------------------------------------------
!                Variable Declaration and Description Section

! Overall Routine Declaration:
    subroutine lambda_log_disp( x0, k, m, L, ie_value  &
                              , tan_motion_lim  &
                              , normalize_input  &
                              , out_lambda, out_log_disp  &
                              , out_npts, out_shells )

! Include Files:
!   None.

! Modules Used and Module Variables Changed by Subroutine:
!   None.

! Input Arguments:
    real, intent(in)  &         !- Input scalar timeseries.  Real vector.
       :: x0(:)

    integer, intent(in)  &      !- Evolution times at which to calculate
       :: k(:)                  !  time-dependent exponent.  Integer
                                !  vector of equally spaced values of k,
                                !  starting with 1.

    integer, intent(in)  &      !- Embedding dimension.  Scalar integer.
       :: m

    integer, intent(in)  &      !- Delay time.  Scalar integer.
       :: L

    integer, intent(in)  &      !- Ruler initialization (to give rl(1)).
       :: ie_value              !  This variable defines shells used in
                                !  calculations and output in out_shells;
                                !  each shell is defined as the interval:
                                !     ( 2**[-(i+1)/2], 2**[-i/2] )
                                !  where i is an integer from ie_value to
                                !  ie_value+8.

    integer, intent(in)  &      !- Value of the criteria to remove tangen-
       :: tan_motion_lim        !  tial motion effects (i.e. ABS(i-j) must
                                !  be greater than this value to be consid-
                                !  ered for calculations).  If parameter 
                                !  is set to < 0, the subroutine uses a 
                                !  value of 10 for the tangential effects
                                !  removal test.

    integer, intent(in)  &      !- If value is 1, before using the input
       :: normalize_input       !  time0 series x0, the subroutine normal-
                                !  izes it by: (x0 - min)/(max - min), 
                                !  where min and max are the maximum and 
                                !  minimum of x0, respectively (i.e. normal-
                                !  ize to unit interval [0,1]).

! Output Arguments:
    real, intent(out)  &           !- Each row is the time-dependent expon-
       :: out_lambda(SIZE(k),9)    !  ent (lambda) at each evolution time k,
                                   !  for shells given by the shell interval
                                   !  (out_shells(i+1), out_shells(i)).

    real, intent(out)  &           !- Each row is the logarithmic displace-
       :: out_log_disp(SIZE(k),9)  !  ment at each evolution time k, for
                                   !  shells given by the shell interval
                                   !  (out_shells(i+1), out_shells(i)).

    integer, intent(out)  &        !- Number of point pairs used for the en-
       :: out_npts(SIZE(k),9)      !  semble averaging in each shell at each 
                                   !  evolution time k, for shells given
                                   !  by the shell interval (out_shells(i+1),
                                   !  out_shells(i)).

    real, intent(out)  &           !- Shell upper boundaries for each row
       :: out_shells(9)            !  in out_lambda and out_log_disp.  The
                                   !  bounds for each i row in those output 
                                   !  variables is the interval:
                                   !     ( out_shells(i+1), out_shells(i) ).
                                   !  The final row (i=9) corresponds to
                                   !  the interval:
                                   !     ( 0, out_shells(9) )
                                   !  which is a "shell" that is actually a
                                   !  "ball."

! Local Variables:

    real, dimension(SIZE(k)+1,50) :: a
    real, dimension(SIZE(k)+1) :: d
    real, dimension(SIZE(x0)) :: xy
    real, dimension(SIZE(x0)) :: x           !- normalized ver. of x0
    real    :: r(SIZE(k)+1)
    real    :: ssz(SIZE(k),9)
    real    :: rl(9)
    integer :: mapped(SIZE(k))
    integer :: nsz(SIZE(k),9)

    integer :: nxp                     !- num pts. in x0 time series
    integer :: map                     !- num pts. in evol. time vector
    integer :: nxp_p, nxp1, map1
    real    :: cmax, cmin, ci
    integer :: tan_motion_lim0
    integer :: nn                      !- num pts. actually used in calcs.
    integer :: ie, in, md, mt
    integer :: mstep, mmap
    integer :: i
!============================================================================




! ------------------------ Initial Parameter Setting ------------------------

if (tan_motion_lim .ge. 0) then      !- set test to remove tangential motion
   tan_motion_lim0 = tan_motion_lim  !  effects
else
   tan_motion_lim0 = 10
endif 

nxp   = SIZE(x0)     !- set some array bounds etc.
nxp_p = SIZE(x0)
nxp1  = SIZE(x0)
map   = SIZE(k)
map1  = map + 1

x = x0               !- set x to x0 (normalize below, if chosen)




! ----------------------- Normalize Input (If Chosen) -----------------------

if (normalize_input .eq. 1) then

	cmax=x(1)
	cmin=x(1)
	do 301 i=2,nxp
	ci=x(i)
	if(cmax.gt.ci) goto 302
	cmax=ci
302	if(cmin.lt.ci) goto 301
	cmin=ci
301	continue
	do 303 i=1,nxp
303	x(i)=(x(i)-cmin)/(cmax-cmin)

endif




! ------------------------- Main Parameters Setting -------------------------

!- set Jianbo's parameter ie using variables that are controlled via 
!  calling arguments, and calculate rulers rl:

ie = ie_value

	write(*,*) 'ie=             ',ie
	if (ie.lt.0) goto 4001
	rl(1)=1./2**ie
	goto 4002
4001	rl(1)=2**(-ie)
4002	write(*,*) '                 rl(1)=',  rl(1)
	do 113 i=2,9
113	rl(i)=rl(i-1)/2


!- set a few of Jianbo's parameters, using variables that are controlled
!  via calling arguments:

in = tan_motion_lim0   !- points for uncorrelation of tangential motion effects
md = m                 !- embedding dimension
mt = L                 !- delay time
mstep = k(2)-k(1)      !- step of k
mmap = map             !- number of k values

if (md .le. 0) then
   write(*,*)'error--bad md value'
   stop
endif


!- calculate nn, the number of points for calculation.  note that this is
!  not the same as the number of points in a timeseries, because you don't
!  necessarily need all those points to span the phase space.  nn is set
!  to the calculated value or 8000, whichever is smaller:

nn = nxp - ( in + ((md-1) * mt) + (mstep * mmap) )
nn = MIN(nn, 8000)

if (nn .le. 0) then
   write(*,*)'error--bad nn value'
   stop
endif


!- diagnostic messages:

      write(*,*) 'num pts in input series  :  nxp', nxp
	write(*,*) 'points for uncorrelation :  in', in
	write(*,*) 'total points for calc.   :  nn', nn
	write(*,*) 'embedding dimension      :  md', md
	write(*,*) 'delay time               :  mt', mt
	write(*,*) 'step of k                :  mstep', mstep
	write(*,*) 'num of k values          :  mmap', mmap




! -------------------------- Main Computation Block -------------------------

	do 3101 i=1,mmap
3101	mapped(i)=mstep*i
	do 3102 i=1,mmap
	do 3102 j=1,9
	ssz(i,j)=0.0
	nsz(i,j)=0
3102	continue

	nn1=nn+in
	do 1 i0=1,mt
	do 4 j0=i0+in,nn1
	i=i0
	j=j0
	do 104 kk=1,md
	i1=i+(kk-1)*mt
	j1=j+(kk-1)*mt
104	a(map1,kk)=(x(i1)-x(j1))*(x(i1)-x(j1))
339	d(map1)=0.0
	do 105 kk=1,md
105	d(map1)=d(map1)+a(map1,kk)
	dd=d(map1)
	if(dd.le.0) goto 41
	if(dd.ge.rl(1)) goto 41

	il=1
	if(dd.ge.rl(2)) goto 1000
	il=2
	if(dd.ge.rl(3)) goto 1000
	il=3
	if(dd.ge.rl(4)) goto 1000
	il=4
	if(dd.ge.rl(5)) goto 1000
	il=5
	if(dd.ge.rl(6)) goto 1000
	il=6
	if(dd.ge.rl(7)) goto 1000
	il=7
	if(dd.ge.rl(8)) goto 1000
	il=8
	if(dd.ge.rl(9)) goto 1000
	il=9

1000	r(map1)=ALOG(dd)
	imp=1
9001	net=mapped(imp)
	do 1041 kk=1,md
	i1=i+(kk-1)*mt
	j1=j+(kk-1)*mt
1041	a(imp,kk)=(x(i1+net)-x(j1+net))*(x(i1+net)-x(j1+net))
	d(imp)=0.0
	do 1042 kk=1,md
1042	d(imp)=d(imp)+a(imp,kk)
	d2=d(imp)
	if(d2.le.0) goto 411
	r(imp)=ALOG(d2)
	er=r(imp)-r(map1)
	ssz(imp,il)=ssz(imp,il)+er/2
	nsz(imp,il)=nsz(imp,il)+1
411	imp=imp+1
	if(imp.le.mmap) goto 9001

41	if(j.ge.nn1) goto 4
	i=i+mt
	j=j+mt
	do 106 kk=1,md-1
106	a(map1,kk)=a(map1,kk+1)
	ldm=mt*(md-1)
	ii1=i+ldm
	jj1=j+ldm
	a(map1,md)=(x(ii1)-x(jj1))*(x(ii1)-x(jj1))
	goto 339
4	continue
1	continue
	do 7001 il=1,9
	do 7001 i=1,mmap
	mm=nsz(i,il)
	if(mm.eq.0) goto 7001
	ssz(i,il)=ssz(i,il)/nsz(i,il)
7001	continue




! ------------------------ Prepare Output Variables -------------------------

do i=1,mmap
   out_log_disp(i,1) = ssz(i,1)-(ie+0)*ALOG(2.)/2.
   out_log_disp(i,2) = ssz(i,2)-(ie+1)*ALOG(2.)/2. 
   out_log_disp(i,3) = ssz(i,3)-(ie+2)*ALOG(2.)/2.
   out_log_disp(i,4) = ssz(i,4)-(ie+3)*ALOG(2.)/2.
   out_log_disp(i,5) = ssz(i,5)-(ie+4)*ALOG(2.)/2.
   out_log_disp(i,6) = ssz(i,6)-(ie+5)*ALOG(2.)/2.
   out_log_disp(i,7) = ssz(i,7)-(ie+6)*ALOG(2.)/2.
   out_log_disp(i,8) = ssz(i,8)-(ie+7)*ALOG(2.)/2.
   out_log_disp(i,9) = ssz(i,9)-(ie+8)*ALOG(2.)/2.
enddo

out_lambda = ssz
out_npts   = nsz
out_shells = SQRT(rl)




end subroutine lambda_log_disp

!======== end of file ========
