 !    Reshavanje j-ne kretanja za 1D LHO - Runge-Kutta metoda 
!     Program
  Program LHO_rk 
!  zadavanje pochetnih uslova
real*8, parameter:: k=1.0d0
   real*8, parameter:: x0=1.d0, v0=0.d0, tk=20.d0, dt=0.001, tp=0 !pochetni polozhaj, brzina
   ! konachno vreme, korak u vremenu, pochetno vreme
! deklaracija tipa svih promenljivih 
 	  real*8      x(2), y(2), xout(2), energija, yy  ! redom : trenutni polozhaj, brzina, ubrzanje, energija, f-ja koja ce se izrachunati
	  real*8 t ! vreme
! Otvaranje fajla (datoteke) u koji ce se smestiti rezultati izrachunavanja
    open(30,file='lhork4.dat',status='unknown')
! prorachunavanje polozhaja, brzine i energije
  t=tp
  x(1)=x0
  x(2)=v0
  yy=1.
! rachunati sve dok vreme ne bude jednako konachnom vremenu
  do while(t.le.tk)
! energija lho jedinicne mase
  energija=1.d0/2.d0*x(2)*x(2)+1.d0/2.d0*k*x(1)*x(1)
! egzaktna trajektorija lho sa parametrima izabranim u knjizi i programu
    yy=dcos(t)
!  ispisivanje podataka u datoteku, 100 je broj naredbe koja pokazuje kako se ispisuju podaci
  write(30,100) t,x(1),x(2),energija,yy
  call rk4(x,t,dt)
  enddo
100 format(10(e13.7,2x))
  stop
end	

subroutine derivs(t,x,y)
real*8, parameter:: k=1.0d0
real*8 x(2),y(2),t
y(1)=x(2)
y(2)=-k*x(1)
return
end


! Runge-Kutta 4tog reda
SUBROUTINE rk4(y,x,h)
INTEGER, parameter:: n=2,NMAX=2
REAL*8 h,x,dydx(n),y(n),yout(n)
!EXTERNAL derivs
!Given values for the variables y(1:n) and their derivatives dydx(1:n) known at x, use
!the fourth-order Runge-Kutta method to advance the solution over an interval h and return
!the incremented variables as yout(1:n), which need not be a distinct array from y. The
!user supplies the subroutine derivs(x,y,dydx), which returns derivatives dydx at x.
INTEGER i
REAL*8 h6,hh,xh,dym(NMAX),dyt(NMAX),yt(NMAX)
hh=h*0.5
h6=h/6.
call derivs(x,y,dydx)
do  i=1,n !First step.
yt(i)=y(i)+hh*dydx(i)
enddo 
xh=x+hh
call derivs(xh,yt,dyt) !Second step.
do  i=1,n
yt(i)=y(i)+hh*dyt(i)
enddo 
call derivs(xh,yt,dym) !Third step.
do  i=1,n
yt(i)=y(i)+h*dym(i)
dym(i)=dyt(i)+dym(i)
enddo 
call derivs(x+h,yt,dyt) !Fourth step.
do  i=1,n !Accumulate increments with proper weights.
yout(i)=y(i)+h6*(dydx(i)+dyt(i)+2.*dym(i))
y(i)=yout(i)
enddo 
x=x+h
return
END
