! Program pentru BDF - cu IMSL: DIVPAG
Program BDF_IMSL

use work
use DFLIB
!USE Numerical_Libraries

use ifwin

include 'link_fnl_shared.h'


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

character(max_path) lpstrTitle
integer status

character($maxpath) fname1, fname2, fname3
real*4 s

external derivs, Jfcn
integer, parameter :: nparam = 50
real(8) param(nparam)

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

lpstrTitle ='Selecteaza Fisierul de intrare'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
allocate(x0(n))

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

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(xEnd(n))

! Inital conditions:
!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(0,:) =x0
x(1:No_Intervals,:) =0.0
xEnd=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

! Initialize

CALL dSET (NPARAM, 0.0D0, PARAM, 1)
ido =1

if(TOL <=0) TOL =1D-10
param(1) =h
param(4) =1000000
param(10) =1	

!param(12) =1	! ADAMS-Multon
param(12) =2	! BDF (Gear)

!param(13) =0	! fixed-point iteration
param(13) =1	! Provided Jacobian
!param(13) =2	! Approximated Jacobian (divided difference)


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

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

CALL dIVPAG (IDO, N, derivs, Jfcn, AM, t, TEND, TOL, PARAM, xEnd)

x(i+1,:) =xEnd
ENDDO


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

s = secnds(s)
print '(/, a,f10.2,a)',' Integration Time=', s,' sec.'
write(1,'(/,a, I10)')  'No.of Function Evaluations:', int4(param(35))
write(1,'(a, I10, /)') 'No.of Steps taken:', int4(param(34))
write(1,'(a, 1pE10.2)')   'Current trial step: ', param(31)
write(1,'(a, 1pE10.2)')   'Computed minimum step: ', param(32)
write(1,'(a, 1pE10.2,/)') 'Computed maximum step: ', param(33)

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

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

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

close(2)
ENDIF

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

close (1)
deallocate(x)
!deallocate(fa)

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

600 continue
END Program BDF_IMSL


