! Program pentru Predictor-Corector ADAMS - ordinul 4 si 5
! Problema celor doua corpuri

Program Adams45

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, fname4, fname5, dir
character($maxpath) fname6, fname7, fname8		! fazX, fazY, orbita

real(4) s

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

real(8) xExact(4), xError(4)
integer order

call GetFile(fname1, status)

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

! Read dummy; data; title.
read(1,'(a)') title
read(1,*) n		! nr. ec. de ord 1
read(1,*) e
read(1,*) t0, h, TOL, order		! ordin =4,5

if(order >5) then
order =5
endif

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

allocate(x0(n))
! Inital conditions:
x0(1) =1 -e
x0(2) =0.
x0(3) =0.
x0(4) =dsqrt( (1+e)/(1-e) )

No_Intervals =tEnd/h +0.5

print '(1x, a)', 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))
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.;


! 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'
fname4 =fname1(1:iPoint-1)//'-Y.dat'
fname5 =fname1(1:iPoint-1)//'-Yprime.dat'

open(2, file=fname2, status='unknown')
open(3, file=fname3, status='unknown')
open(4, file=fname4, status='unknown')
open(5, file=fname5, status='unknown')

	if(cod >=2) then
	fname6 =fname1(1:iPoint-1)//'-XFaz.dat'
	fname7 =fname1(1:iPoint-1)//'-YFaz.dat'
	fname8 =fname1(1:iPoint-1)//'-orb.dat'
	open(6, file=fname6, status='unknown')
	open(7, file=fname7, status='unknown')
	open(8, file=fname8, status='unknown')
	endif
ENDIF


if(order ==4) then
kend = 3			!allocate(fa(-3:1,n))
else
kend =4				!allocate(fa(-4:1,n))
endif

allocate(fa(-kend:1,n))
!if(order ==4) then
!allocate(fa(-3:1,n))
!else
!allocate(fa(-4:1,n))
!endif

! START VALUES - by R-K:
! x0 = exista; x1, x2, x3
xEnd =x0

!kend =3; if(order ==5) kend =4

print '(/, a, I3)', ' ADAMS Metod, order', order
print '(a)', ' * Working ...'

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

do i = 0, kend-1		! x1, x2, x3 (,x4)
t =t0+ i*h
call RKSYST(t)
	x(i+1,:) =xEnd
	x0 =xEnd
enddo

deallocate(w1,w2,w3,w4)		! End R-K
!nEval =12; if(order ==5) nEval =nEval +4
nEval =kend*4

! Pregateste primele "fa": (-4), -3, -2, -1, 0:
!iend1 = iend+1	!3; 4

do ia =0, kend	!3;4
xEnd =x(ia,:)
ta = t0 +ia*h
call derivs(ta)
la =ia -kend	!la =ia -iend
fa(la,:) =f(:)	! la =(-4,) -3, -2, -1, 0.
enddo

!nEval =nEval +4; if(order ==5) nEval =nEval+1
nEval =nEval +kend +1

DO i =kend, No_Intervals-1
ta =t0+ i*h

if(i ==kend) goto 150
! Pregateste fa: (-4);-3, -2, -1 = Translatie:
!if(order ==5) then
!fa(-4,:) =fa(-3,:)
!endif
!fa(-3,:) =fa(-2,:); fa(-2,:)=fa(-1,:); fa(-1,:)=fa(0,:);

do la = -kend, -1
fa(la,:) =fa(la+1,:)
enddo

150 continue
! Calculeaza noul fa(0):
xEnd =x(i,:)
call derivs(ta)
fa(0,:) =f(:)

nEval =nEval +1

x0 =x(i,:)
t =ta+h
if(order ==4) then
call ADAMS4(t, TOL, iter)		!TOL =0 => 1 iteratie
else
call ADAMS5(t, TOL, iter)
endif

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


s = secnds(s)

print '(/, a,f10.2,a)',' Integration Time=', s,' sec.'
print '(a, I)', ' No. of function evaluations: ', nEval

write(1, '(/, a, I)') ' No.of Intervals:', No_Intervals
write(1, '(/, a, I3)') ' ADAMS Method, order', order
write(1,'(a, I)') ' No. of function evaluations: ', nEval

write(1,'(/,a,2(1pE16.7))') ' End: X,  Y : ', x(No_intervals,1:2)
write(1,'(a,2(1pE16.7))')   '      X'', Y'': ', x(No_intervals,3:4)

IF(cod /=0) THEN
iStart =tStart_H/h
iEnd =tEnd_H/h
print*,' Writing history ...'
do i=iStart, iEnd, iPas
write(2,*) i*h, x(i,1)	!X
write(3,*) i*h, x(i,3)	!Xprime
write(4,*) i*h, x(i,2)	!Y
write(5,*) i*h, x(i,4)	!Yprime

	if(cod >=2) then
	write(6,*) x(i,1), x(i,3) !XFaz
	write(7,*) x(i,2), x(i,4) !YFaz
	write(8,*) x(i,1), x(i,2) !Orbita
	endif
enddo

close(2)
close(3)
close(4)
close(5)
	if(cod >=2) then
	close(6)
	close(7)
	close(8)
	endif
ENDIF

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

! Now, calculate the error at i*T, where T=Period:
True_If: IF(true_solution) THEN

Period =2*PI
nStepinPeriod= Period/h
NoofPeriod =No_Intervals/nStepinPeriod

write(1, '(/,a)') ' Error in xEnd, yEnd, x''End, y''End:'
write(1, '(a)')   ' Period	Time		Error'

I_DO: DO I =1, NoofPeriod
IStep =I*nStepinPeriod
tEnd =t0 +IStep*h
call exact(tEnd, xExact)

! Error:
xError =xExact -x(IStep,:)
write(1,'(I3,2x, f10.6, 2x, 4(2x,1pD10.3))') I, tEnd, xError
END DO I_DO

ENDIF True_If

close (1)
deallocate(x)
deallocate(x0, xEnd, f)
deallocate(fa)

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


END Program Adams45

