c	By Liang Huang, Arizona State University, USA; Lanzhou University, China

c	The paper: L. Huang, Q. Chen, Y.-C. Lai, L. M. Pecora, ``Generic behavior of master-stability 
c	functions in coupled nonlinear dynamical systems,'' Physical Review E 80, 036204 (2009). 

c	This program calculates all the three Lyapunov exponents of coupled Rossler ossilators, 
c	where the largest is the master stability function.
c	The coupling is from ij1 to ij2 (see below).


c	include 'ddqrdcmp.f'
c	include 'ddshell.f'
c	USE numerical_libraries
	parameter(n=1,ni=3,h=0.001,a=0.2,b=0.2,c=9.,nd=n*ni,nt=1000000,
     *	n1=3000,nl=13000,ih=100)
      real*8 t,u(nd),k1(nd),k2(nd),k3(nd),k4(nd),v(nd),v1(nd),vd,r0
	real*8 tv(nd),A0(ni,ni),Od(ni,ni),Ot(ni,ni),Os(ni,ni),Ost(ni,ni)
	real*8 Ot1(ni,ni),Ot2(ni,ni),Ot3(ni,ni),Ot4(ni,ni)
	real*8 ly(nd),lyt(nd),lyma,Abk(ni,ni)
	Complex*16 eval(ni)
	LOGICAL sing
	real*8 cc(ni),dd(ni),QQ(ni,ni),RR(ni,ni),Rt(ni,ni),Qt(ni,ni)
	real*8 aph0(1000)
	character fn(7)

c	ij1 to ij2 coupling, coupled at A0(ij2,ij1)
	ij1=1
	ij2=1
c	This is Fig. 1(a).

c	do ij1=1,ni
c	do ij2=1,ni
	
	fn(1)=char(ij1+48)
	fn(2)=char(ij2+48)


	open (2,file='MSFRosslerAll.dat')


c	aph is "K" in the paper.
	it=0
	it1=0
	do i=1,1000
		aph=(i-1)*0.02
		if (it.eq.0 .and. aph.gt.1) it=i
		if (it.gt.0) aph=1+(i-it+1)*0.1
	
		if (it1.eq.0 .and. aph.gt.10) it1=i
		if (it1.gt.0) aph=10+(i-it1+1)*2

		aph0(i)=aph

		write (*,*) i,aph
		if (aph.gt.10) exit
	end do

	nss=i-1



	do iss=1,nss
		
		aph=aph0(iss)
	


	it0=50000*100
	ijk=0

c	write (*,*) a,b,c,aph
c	pause



	
	itt=0
	ly(1)=0
	ly(2)=0
	ly(3)=0

	

ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c	The constant values of the Jacobian matrix
	Abk(1,1)=0
	Abk(1,2)=-1
	Abk(1,3)=-1
	Abk(2,1)=1
	Abk(2,2)=a
	Abk(2,3)=0
	Abk(3,2)=0

	do i=1,ni
		do j=1,ni
			Od(j,i)=0
		end do
		Od(i,i)=1
	end do


	
	do i=1,ni
		do j=1,ni
			Os(j,i)=0
			RR(j,i)=0
			QQ(j,i)=0
		end do
		Os(i,i)=1
		RR(i,i)=0
		QQ(i,i)=1
		lyt(i)=0
	end do


	
ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccccccccccccccccccccccccccccccccccccccccccccccccccccccc



	vd=0
	do i=1,nd
	  u(i)=rand()
	  v(i)=0.1*rand()
		vd=vd+v(i)**2
	end do
	v(nd)=sqrt(1-vd)


	ip=0
	t=0
	do i=1,100000000
	 u0=u1
	 u1=u(1)
	 t=i*h

c	RK4 integration for both the original system and dO/dt=DF*O, see Ott's book, p143
	 call rk(u,k1,nd,n,ni,h,a,b,c)	 

c		 Evolve the system for n1 oscillations, then integrate the matrix equation dO/dt=DF*O	 
		 if (ip.ge.n1) then		
			
c			The elements of the Jacobian matrix which depend on the dynamical values			
			Abk(3,1)=u(3)
			Abk(3,3)=u(1)-c

			do j1=1,ni
				do j=1,ni
					A0(j1,j)=Abk(j1,j)
				end do
					A0(j1,j1)=Abk(j1,j1)-aph
			end do

c			A0(ij2,ij1)=Abk(ij2,ij1)-aph

			do j1=1,ni
				do j=1,ni
					Ot1(j1,j)=0
					do j2=1,ni
						Ot1(j1,j)=Ot1(j1,j)+A0(j1,j2)*Od(j2,j)
					end do
					Ot1(j1,j)=Ot1(j1,j)*h
				end do
			end do

			do j=1,ni
				do j1=1,ni
					Ot(j1,j)=Od(j1,j)+0.5*Ot1(j1,j)
				end do
			end do
		end if

	 do j=1,nd
	   tv(j)=u(j)+0.5*k1(j)
	 end do



	 call rk(tv,k2,nd,n,ni,h,a,b,c)

		 if (ip.ge.n1) then	
			
			
			Abk(3,1)=tv(3)
			Abk(3,3)=tv(1)-c

			do j1=1,ni
				do j=1,ni
					A0(j1,j)=Abk(j1,j)
				end do
					A0(j1,j1)=Abk(j1,j1)-aph
			end do


c			A0(ij2,ij1)=Abk(ij2,ij1)-aph

			do j1=1,ni
				do j=1,ni
					Ot2(j1,j)=0
					do j2=1,ni
						Ot2(j1,j)=Ot2(j1,j)+A0(j1,j2)*Ot(j2,j)
					end do
					Ot2(j1,j)=Ot2(j1,j)*h
				end do
			end do

			do j=1,ni
				do j1=1,ni
					Ot(j1,j)=Od(j1,j)+0.5*Ot2(j1,j)
				end do
			end do
		end if

	 do j=1,nd
	   tv(j)=u(j)+0.5*k2(j)
	 end do

	


	 call rk(tv,k3,nd,n,ni,h,a,b,c)

	
		 if (ip.ge.n1) then	
			
			Abk(3,1)=tv(3)
			Abk(3,3)=tv(1)-c

			do j1=1,ni
				do j=1,ni
					A0(j1,j)=Abk(j1,j)
				end do
					A0(j1,j1)=Abk(j1,j1)-aph
			end do


c			A0(ij2,ij1)=Abk(ij2,ij1)-aph

			do j1=1,ni
				do j=1,ni
					Ot3(j1,j)=0
					do j2=1,ni
						Ot3(j1,j)=Ot3(j1,j)+A0(j1,j2)*Ot(j2,j)
					end do
					Ot3(j1,j)=Ot3(j1,j)*h
				end do
			end do

			do j=1,ni
				do j1=1,ni
					Ot(j1,j)=Od(j1,j)+Ot3(j1,j)
				end do
			end do
		end if


	 do j=1,nd
	   tv(j)=u(j)+k3(j)
	 end do

	 call rk(tv,k4,nd,n,ni,h,a,b,c)


		 if (ip.ge.n1) then
			
			Abk(3,1)=tv(3)
			Abk(3,3)=tv(1)-c

			do j1=1,ni
				do j=1,ni
					A0(j1,j)=Abk(j1,j)
				end do
					A0(j1,j1)=Abk(j1,j1)-aph
			end do


c			A0(ij2,ij1)=Abk(ij2,ij1)-aph

			do j1=1,ni
				do j=1,ni
					Ot4(j1,j)=0
					do j2=1,ni
						Ot4(j1,j)=Ot4(j1,j)+A0(j1,j2)*Ot(j2,j)
					end do
					Ot4(j1,j)=Ot4(j1,j)*h
				end do
			end do

			do j=1,ni
				do j1=1,ni
					Od(j1,j)=Od(j1,j)+
     *				(Ot1(j1,j)+2*Ot2(j1,j)+2*Ot3(j1,j)+Ot4(j1,j))/6
				end do
			end do
		end if


	 do j=1,nd
	   u(j)=u(j)+(k1(j)+2*k2(j)+2*k3(j)+k4(j))/6.
	 end do

c	 ip: number of oscillations
	
	 if (u1.gt.u0 .and. u1.gt.u(1)) ip=ip+1
	 
c	 calculate the Lyapunov exponent for nl oscillations, then exit
	 if (ip.ge.n1+nl) exit

	 if (ip.ge.n1) then	
c			A0(3,1)=u(3)
c			A0(3,3)=u(1)-c

c	Eular integration is not good, leads to systematic errors			
c			do j1=1,ni
c				do j=1,ni
c					Ot(j1,j)=0
c					do j2=1,ni
c						Ot(j1,j)=Ot(j1,j)+A0(j1,j2)*Od(j2,j)
c					end do
c					Ot(j1,j)=Ot(j1,j)*h+Od(j1,j)
c				end do
c			end do

c			do j=1,ni
c				do j1=1,ni
c					Od(j1,j)=Ot(j1,j)
c				end do
c			end do


			
c	do QR decomposition	 
	   if (mod(i,ih).eq.0) then
			ijk=ijk+1

			do j1=1,ni
				do j=1,ni
					Ot(j1,j)=0
					do j2=1,ni
						Ot(j1,j)=Ot(j1,j)+Od(j1,j2)*QQ(j2,j)
					end do
				end do
			end do



c	        QR decomposition, Rev. Mod. Phy. 57, No. 3, Part 1, p. 617, 1985
			call ddqrdcmp(Ot,ni,ni,cc,dd,sing)
			
			do j=1,ni
				do j1=1,ni
					Rt(j,j1)=0
				end do
			end do
			do j=1,ni-1
				Rt(j,j)=dd(j)
				do j1=j+1,ni
					Rt(j,j1)=Ot(j,j1)
					Ot(j,j1)=0
				end do
			end do
			j=ni
			Rt(j,j)=dd(j)
cc	obtained R
			
			do j=1,ni
				do j1=1,ni
					QQ(j1,j)=0
				end do
				QQ(j,j)=1
			end do

			do j=1,ni-1
				do j1=1,ni
					do j2=1,ni
						Qt(j1,j2)=-Ot(j1,j)*Ot(j2,j)/cc(j)
						if (j1.eq.j2) Qt(j1,j2)=Qt(j1,j2)+1
					end do
				end do
					
				do j0=1,ni
				do j1=1,ni
					Ost(j0,j1)=0
					do j2=1,ni
						Ost(j0,j1)=Ost(j0,j1)+QQ(j0,j2)*Qt(j2,j1)
					end do
				end do
				end do	

				do j0=1,ni
				do j1=1,ni
					QQ(j0,j1)=Ost(j0,j1)
				end do
				end do
			end do
cc	obtained Q

c			test QR			
c			do j0=1,ni
c				do j1=1,ni
c					Ost(j0,j1)=0
c					do j2=1,ni
c						Ost(j0,j1)=Ost(j0,j1)+QQ(j0,j2)*Rt(j2,j1)
c					end do
c					write (*,*) j0,j1,Ost(j0,j1)-Od(j0,j1)
c				end do
c			end do	
ccccccccccccccccccccc
			
c	matrix product	
c			do j1=1,ni
c				do j=1,ni
c					Ot(j1,j)=0
c					do j2=1,ni
c						Ot(j1,j)=Ot(j1,j)+Rt(j1,j2)*RR(j2,j)
c					end do
c				end do
c			end do

	   	

c	Since R is uptriangular, the diagnals of the product equals the product of 
c	the corresponding diagnals, then log(product)=sum(log)
				do j=1,ni
					RR(j,j)=RR(j,j)+log(abs(Rt(j,j)))
				end do

c	reset
		do it=1,ni
		do j=1,ni
			Od(j,it)=0
		end do
		Od(it,it)=1
		end do


	   do j=1,nd
		  ly(j)=RR(j,j)
	   end do

c	sort in ascending order
c		call shell(nd,ly)

c	   if (mod(ijk,100).eq.0)write(2,100) t,(ly(j)/ijk/ih/h,j=1,ni)

c			if (mod(ijk,100).eq.0) then			
c				write (*,*) t,(ly(j)/ijk/ih/h,j=1,ni)
c			end if





	 end if

	end if
	 
c	  or if t is larger than nt, exit
	  if (t.ge.nt) exit	
	end do
	

	write (2,'(99e20.12)') aph,(ly(j)/ijk/ih/h,j=1,ni)

	write (*,*) ij1,ij2,aph,t,(ly(j)/ijk/ih/h,j=1,ni)
c	close (1)
c	pause
	end do

c	end do
c	end do
	
	close (2)
c	pause

100	format (101e20.12)
	end


	FUNCTION f(j,tv,i,a,b,c,ni,nd)
	real*8 tv(nd)
	if (j.eq.1) then
		f=-(tv((i-1)*ni+2)+tv((i-1)*ni+3))
c	    dx/dt=-(y+z)
	else if (j.eq.2) then
		f=tv((i-1)*ni+1)+a*tv((i-1)*ni+2)
c	    dy/dt=x+ay
	else if (j.eq.3) then
		f=b+tv((i-1)*ni+3)*(tv((i-1)*ni+1)-c)
c	    dz/dt=b+z(x-c)	
	end if
	end

	SUBROUTINE bound(tv,n)
	real*8 tv(-1:n+2)
	tv(-1)=tv(n-1)
	tv(0)=tv(n)
	tv(n+1)=tv(1)
	tv(n+2)=tv(2)
	return
	end

	SUBROUTINE colu(u,tv,nd)
	real*8 u(nd),tv(nd)
	do j=1,nd
	tv(j)=u(j)
	end do
	return
	end


	SUBROUTINE rk(tv,tv1,nd,n,ni,h,a,b,c)
	real*8 tv(nd),tv1(nd)
	do i=1,n
	  do j=1,ni
	    tv1((i-1)*ni+j)=h*f(j,tv,i,a,b,c,ni,nd)
	  end do
	end do
c	write (*,*) tv1(i)
c	pause
	return
	end






      SUBROUTINE ddqrdcmp(a,n,np,c,d,sing)
      INTEGER n,np
      REAL*8 a(np,np),c(n),d(n)
      LOGICAL sing
      INTEGER i,j,k
      REAL*8 scale,sigma,sum,tau
      sing=.false.
      do 17 k=1,n-1
        scale=0.
        do 11 i=k,n
          scale=max(scale,abs(a(i,k)))
11      continue
        if(scale.eq.0.)then
          sing=.true.
          c(k)=0.
          d(k)=0.
        else
          do 12 i=k,n
            a(i,k)=a(i,k)/scale
12        continue
          sum=0.
          do 13 i=k,n
            sum=sum+a(i,k)**2
13        continue
          sigma=sign(sqrt(sum),a(k,k))
          a(k,k)=a(k,k)+sigma
          c(k)=sigma*a(k,k)
          d(k)=-scale*sigma
          do 16 j=k+1,n
            sum=0.
            do 14 i=k,n
              sum=sum+a(i,k)*a(i,j)
14          continue
            tau=sum/c(k)
            do 15 i=k,n
              a(i,j)=a(i,j)-tau*a(i,k)
15          continue
16        continue
        endif
17    continue
      d(n)=a(n,n)
      if(d(n).eq.0.)sing=.true.
      return
      END




      SUBROUTINE shell(n,a)
      INTEGER n
      REAL*8 a(n)
      INTEGER i,j,inc
      REAL*8 v
      inc=1
1     inc=3*inc+1
      if(inc.le.n)goto 1
2     continue
        inc=inc/3
        do 11 i=inc+1,n
          v=a(i)
          j=i
3         if(a(j-inc).gt.v)then
            a(j)=a(j-inc)
            j=j-inc
            if(j.le.inc)goto 4
          goto 3
          endif
4         a(j)=v
11      continue
      if(inc.gt.1)goto 2
      return
      END