! Metoda Runge-Kutta de ordinul 4.
! Sistem de ordinul 1
! X' =F(X)
! X ={x1, ..., xn}; F ={f1(X), ..., fn(X)}
! Membrii 2 "fi(X)", i=1,n se definesc in subrutina "derivs":
! X se stocheaza in tabloul "xEnd".
! X0 = conditiile initiale --> in tabloul "x0"
!
! Fisier de intrare (si iesire):
! 1) Titlu:     <= 80 chars.
! 2) n:         nr. de ecuatii
! 3) 4 parametri (ai sistemului): Param%k, Param%k1, Param%k3, Param%B 
! 4) t0, h:     mom. initial; pasul de integrare
! 5) x0(1:n): conditiile initiale
! 6) tEnd:      timpul final de raspuns
! 7) cod:       #0/0 = scrie/nu scrie - istoria miscarii
!   7-1) Numai pentru cod #0: tStart, tEnd, ipas --> pt. istorie
! 8) true_solution: T/F = Exista/Nu solutia exacta
!   Obs.: daca se introduce T (true) atunci: solutia exacta se defineste in subrutina "True(t, x_true)".
!         Pentru F (false) subrutina "True" poate fi oricare (ex.: vida).
! =======================================================================
! Rez. se scriu in continuarea datelor de intrare.

PROGRAM RK4
! Main for Runge-Kutta 4
use work
use DFLIB

implicit real*8 (a-h,o-z)
character(80) title

!real(8), allocatable :: x0(:), xEnd(:)
character($maxpath) fname1, fname2,fname3
real(4) s

logical true_solution
real(8), allocatable :: x_True(:)
integer status


call GetFile(fname1, status)
if(status <=0) goto 600

open(1,file =fname1,status ='old', form ='formatted')

! Read: title; data.
read(1,'(a)') title

read(1,*) n
allocate(x0(n))
x0 =0D0

read(1,*) Param%k, Param%k1, Param%k3, Param%B

read(1,*) t0, h

read(1,*) x0	! cond.init. pt. sistem

read (1,*) tEnd

read(1,*) cod	! 0, #0 - for writing X, Xprime, Faz
	IF(cod /=0) THEN
	read(1,*) tStart, tEnd_2, iPas
	ENDIF

read(1,*) true_solution

if(true_solution) then
allocate(x00(n))
x00 =x0
endif

No_Intervals =(tEnd -t0)/h +0.5

print '(/, a)', Title
print *, 'n:', n
print *, 't0, tEnd, Step:', t0, tEnd, h
print *, 'Initial cond.(system):',x0
print *, 'No.of Intervals:', No_Intervals

allocate (x(0:No_Intervals, n))
allocate (xEnd(n), f(n))
allocate(w1(n),w2(n),w3(n),w4(n))
! Initialize:

x(0,:) =x0
x(1:No_Intervals,:) =0.0

xEnd=0.; ; f=0.;
w1=0.; w2=0.; w3=0.; w4=0.;

print '(/, a)', ' * Working ...'
! Open for time-history:
if(cod /=0) then
iPoint =index(fname1,'.',back=.TRUE.)
fname2 =fname1(1:iPoint-1)//'-X.dat'
fname3 =fname1(1:iPoint-1)//'-XPrime.dat'

open(2, file=fname2, status='unknown')
open(3, file=fname3, status='unknown')
endif

! Start Time:
s = secnds(0.0)
! Time-history integration:
	xEnd =x0
DO i=0, No_Intervals-1
t =t0+ i*h
call RKSYST(t)	!, x0, xEnd)

	x(i+1,:) =xEnd
	x0 =xEnd
ENDDO

s = secnds(s)

print '(a,f10.2,a)',' Integration Time=', s,' sec.'

write(1,'(/, a, I12)') ' No.of Intervals:', No_Intervals

write(*,'(/,a, t8, f10.3)') ' tEnd =',t+h
write(*,'(/,a, (t8, 2(1pE19.10),/))') ' X, X'':', x(No_intervals,1:n)

write(1,'(/,a, t8, f10.3)') ' tEnd =',t+h
write(1,'(/,a, (t8, 2(1pE19.10),/))') ' X, X'':', x(No_intervals,1:n)


if(true_solution) then
allocate(x_True(n)); x_True(:) =0D0;
call True(tEnd, x_True)


write(1, '(a, t8, 2(1pE19.10))') ' True:', x_True

write(*, '(a)') ' Errors at tEnd:'
write(*, '(a, t8, 2(1pG16.4))') ' X, X'':', x_True -x(No_intervals,:)

write(1, '(/, a)') ' Errors at tEnd:'
write(1, '(a, t8, 2(1pG16.4))') ' X, X'':', x_True -x(No_intervals,:)

endif


IF(cod /=0) THEN
print*,' Writing history ...'

iStart = (tStart -t0)/h
iEnd = (tEnd_2 -t0)/h
do i=iStart, iEnd, iPas
write(2,*) t0+i*h, x(i,1)
write(3,*) t0+i*h, x(i,2)
enddo

close(2)
close(3)
ENDIF

write(1,*)
write(1,'(a,f10.2,a)') ' Integration Time=',s,' sec.'
close(1)

deallocate(x)
deallocate(x0, xEnd, f)
deallocate(w1,w2,w3,w4)

print '(a)', ' * Finished.'

600 continue

end program RK4
