! Metoda BDF - ordin "k", pas constant.
! Sistem de ordinul 1:
! x' =f(t,x); x(t0) =x0;
!
! Fisier de intrare (date):
! 1) Titlu:     <= 80 chars.
! 2) n, k:         nr. de ecuatii; ordinul metodei BDF
! 3) 3 parametri (ai sistemului): lambda, A, omega (notatia pt.: VDP)
! 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 (X, Xprime, Faz)
!   7-1) Numai pentru cod #0: tStart_H, tEnd_H, 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).
!
Program BDF

use work
use DFLIB

implicit real(8) (a-h,o-z)
character(80) title
!real(8), parameter :: PI=3.1415926535897932	!384626433832795

character($maxpath) fname1, fname2, fname3, lpstrTitle
integer status

real(4) s

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


! Ordin k:
c =0.
call gen_coef()
! k =5:
! Coeficienti:
!c(0) =137./60.
!c(1) =-5
!c(2) =5
!c(3) =-10./3.
!c(4) =5./4
!c(5) =-1./5

lpstrTitle ='Select Input File'C
call GetFile(fname1, lpstrTitle, status)
if(status <=0) goto 600

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

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

read(1,*) n, k		! k =ordinul metodei BDF
allocate(x0(n))		!, f(n))

read (1,*) lambda, A, omega
read(1,*) t0, h

read(1,*) x0	! cond.init. pt. sistem = calculate
read(1,*) tEnd	! No_Intervals
read(1,*) cod	! 0, #0 - for writing X, Xprime, Faz
	IF(cod /=0) THEN
	read(1,*) tStart_H, tEnd_H, iPas
	ENDIF

read(1,*) true_solution


! Alocari:
allocate(xa(-k+1:1, n))

! Inital conditions - test:
!x0(1) =1.; x0(2) =0.

No_Intervals =(tEnd -t0)/h +0.5

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

allocate (x(0:No_Intervals, n))

! Initialize:
x(1:No_Intervals,:) =0.0
f=0.

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

! Open for time-history:
IF(cod /=0) THEN
iPoint =index(fname1,'.',back=.TRUE.)
fname2 =fname1(1:iPoint-1)//'-X.dat'
open(2, file=fname2, status='unknown')

if(cod >=2) then
fname3 =fname1(1:iPoint-1)//'-Xprime.dat'
open(3, file =fname3, status='unknown')
endif

ENDIF

! Compute Start values:
call Start_Values()

! Start Time:
s = secnds(0.0)
! Time-history integration:

t =t0
x(0,:) =x0
xa(0,:) =x0

DO i=0, No_Intervals-1
tEnd =t0+ (i+1)*h

call sub_BDF(tEnd, kod)

if(kod <0) goto 500

x(i+1,:) =xa(1,:)

do l=-k+1,0
xa(l,:) =xa(l+1,:)
enddo


ENDDO


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

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

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

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


if(true_solution) then
allocate(x_True(n))
call True(tEnd, 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
iStart =(tStart_H -t0)/h
iEnd =(tEnd_H -t0)/h +0.5
print*,' Writing history ...'
do i=iStart, iEnd, iPas
write(2,*) t0+i*h, x(i,1)	!X
if(cod >=2) write(3,*) t0+i*h, x(i,2)	!X'
enddo

close(2)
ENDIF

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

500 continue

close (1)
deallocate(x)
deallocate(xa)

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

600 continue

END Program BDF

