program lylorenz
	parameter(n=3,m=12,st=100)
	integer::i,j,k
    real y(m),z(n),s(n),yc(m),h,y1(m),a,b,r,f(m),k1,k2,k3
	y(1)=10.
	y(2)=1.
	y(3)=0.
	a=10.
    b=8./3.
    r=28.
	t=0.
	h=0.01

!!!!!initial conditions
do i=n+1,m
	y(i)=0.
end do
do  i=1,n
	y((n+1)*i)=1.
    s(i)=0
end do

open(1,file='lorenz1.dat')
open(2,file='ly lorenz.dat')
do 100 k=1,st   !!!!!!!!st iterations
 call rgkt(m,h,t,y,f,yc,y1)

    !!!!normarize vector 123
	z=0.   
    do i=1,n
	 do j=1,n
	z(i)=z(i)+y(n*j+i)**2
	enddo
	if(z(i)>0.)z(i)=sqrt(z(i))
	do j=1,n
    y(n*j+i)=y(n*j+i)/z(i)
	enddo
    end do
   !!!!generate gsr coefficient
   k1=0.
   k2=0.
   k3=0.
do i=1,n	
	k1=k1+y(3*i+1)*y(3*i+2)
	k2=k2+y(3*i+3)*y(3*i+2)
	k3=k3+y(3*i+1)*y(3*i+3)
end do
   !!!!conduct new vector2 and 3
do i=1,n
	y(3*i+2)=y(3*i+2)-k1*y(3*i+1)	
    y(3*i+3)=y(3*i+3)-k2*y(3*i+2)-k3*y(3*i+1)
end do
   
    !!!generate new vector2 and 3,normarize them
do i=2,n
    z(i)=0.
    do j=2,n
	z(i)=z(i)+y(n*j+i)**2
	enddo
    if(z(i)>0.)z(i)=sqrt(z(i))
	do j=2,n
    y(n*j+i)=y(n*j+i)/z(i)
	end do
end do

    !!!!!!!update lyapunov exponent
do i=1,n
   if(z(i)>0)s(i)=s(i)+log(z(i))
enddo

100 continue
do i=1,n
  s(i)=s(i)/(1.*st*h*1000.)
write(2,*)s(i)
enddo
end

!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine rgkt(m,h,t,y,f,yc,y1)

real y(m),f(m),y1(m),yc(m),a,b,r
integer::i,j
do j=1,1000
	call df(m,t,y,f)
	t=t+h/2.0
	do i=1,m
		yc(i)=y(i)+h*f(i)/2.0
		y1(i)=y(i)+h*f(i)/6.0		
		end do			
	call df(m,t,yc,f)
	do i=1,m
		yc(i)=y(i)+h*f(i)/2.0
		y1(i)=y1(i)+h*f(i)/3.0		
	end do	
	call df(m,t,yc,f)
	t=t+h/2.0
	do i=1,m
		yc(i)=y(i)+h*f(i)
		y1(i)=y1(i)+h*f(i)/3.0
	end do
	call df(m,t,yc,f)
	do i=1,m
		y(i)=y1(i)+h*f(i)/6.0
		end do
if(j>500)write(1,*)t,y(1),y(2),y(3)
end do
return
end
!!!!!!!!!!!!!!!!!!!!!!!!
subroutine df(m,t,y,f)
real y(m),a,b,r,f(m)
common a,b,r
    a=10.
    b=8./3.
    r=28.
	f(1)=a*(y(2)-y(1))
	f(2)=y(1)*(r-y(3))-y(2)
	f(3)=y(1)*y(2)-b*y(3)

    do  i=0,2
	f(4+i)=a*y(7+i)-y(4+i)
	f(7+i)=y(4+i)*(r-y(3))-y(7+i)-y(1)*y(10+i)
	f(10+i)=y(2)*y(4+i)-b*y(10+i)+y(1)*y(7+i)
    enddo
	return
end