PROGRAM Central_Dif
!  Dubla Precizie
!  Metoda dif. centrale
!  Eroarea in Deplasare si Viteza !
! u =x(0); u1=x'(0); u2 =x"(0)
! Xant----X----X1:
! Xant =x(t-dt); X =x(t); X1=x(t+dt)

! Ian 2005
!
! Fisierul de intrare - continut:
!===============================
! 1) Titlul problemei
! 2) Tip ecuatie (caracter): "L,l" =Liniara; Alt car.: Neliniar
! 3) m
! 4) zeta, k		! k=0 => k=4*pi^2
! 5) omp, P0		! raspuns fortat: omp<0: => omp =0.8*omega; P0<0: => P0 =omega**2
! 7) dt, tt, npas	! npas =pas de tiparire
! 8) t0, u0, u10	! u10 <0 => u10 =2*pi
! 9) cod_rez		! #0: se scrie si fisierul "nume.dep"
! 10) True_solution ! logic T/F
! ========================================

use dfwin
use work

implicit real(8) (a-h,o-z)
common C10,C20,C1,C2,C11,C21,da,ust

character(80) titlu
character(260) fname1, fname2, fname3, IND*10

data IND/' Er(U) **'/
character(1) lin

real(8) m1, k1, k1ant


call GetFile(fname1, status)

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

iPoint =index(fname1,'.',back=.TRUE.)
fname2 =fname1(1:iPoint-1)//'.rez'
open(2, file =fname2, status='unknown')

read(1, '(a80)') titlu

read(1,*) lin
liniar =lin =='L' .or. lin =='l'

read(1, *) m

read(1, *) zeta,k
if(k ==0.) k=pi4


!IF(liniar) THEN
om2 =k/m
omega =dsqrt(om2)
c =2*m*omega*zeta	!	zeta=c/(2*m*omega)
oma =omega*dsqrt(dabs(1-zeta*zeta))
!ENDIF

read(1, *) P0, omp

if(omp <0.) omp =.8*omega

! se omite daca P0 real este <0.
if(P0 <0)  then
P0 =om2
endif
!

IF(omp /=0) THEN
da =1./(1-omp**2/om2)
ust =P0/k
ENDIF

read(1,*) dt,tt, npas

read(1,*) t0, u0, u10

read (1,*) krez
if(krez /=0) then
read(1,*) npas_dep

iPoint =index(fname1,'.',back=.TRUE.)
fname3 =fname1(1:iPoint-1)//'.dep'
open(3, file =fname3)
endif

read(1,*) true_solution

if(u10 <0.) u10 =pi2

	c20 =u0
	c10 =u10/omega
	if(c ==0.) goto 66
	c2 =u0

	if(zeta <1.) then
	c1 =(u10 +zeta*omega*u0)/oma

	else if(zeta ==1.) then
	c1 =(u10 +zeta*omega*u0)

	else
	u1c =(u10 +zeta*omega*u0)/oma
	c11 =(u0 -u1c)/2.
	c21 =(u0 +u1c)/2.
	endif

66 continue

if(liniar) then
u20 =(p(t0)-c*u10 -k*u0)/m

else	! mu'' + ga(u)u' +fk(u) =p(t)
u20 =(p(t0) -ga(u0)*u10 -fk(u0))/m
endif


write(2, '(/, 1x, a)') trim(titlu)

ltitlu =len_trim(titlu)
10 format(1x, <ltitlu>(1h=), /)
write(2,10)

write(1,*)
write(1, '(a, f10.5)') ' m =', m

if(liniar) then
write(1, '(2(a, 1pg15.7))') ' c =', c, '; k =', k
write(1, '(a, f12.7)') ' omega =', omega
endif

write(1, '(a, f12.7)') ' zeta =', zeta


30 format(/, ' tetaP = ',f12.7,'; P0 = ',f16.7)
if(P0 ==0.) then
write(1,'(/, a)') ' Raspuns liber'

else
write(1,'(/, a)') ' Raspuns fortat'
write(1,30) omp,P0
endif

40	format(/, ' dt = ', f10.4,'; tt =', f10.2, :, '; npas = ', i)
write(1,40) dt, tt, npas

41 format(/, ' t0 =',F10.2, 3x,' u0 =',F10.5, 3x,' u10 =',F10.5, /)
write(1,41) t0, u0, u10

60	format &
('      T         U                   V                            ERR(U)               ERR(V)', /)
write(2,60)

70	format(' ',f10.3, 2(1pG21.11), 2(1pG21.3))
write(2,70) t0, u0, u10, 0., 0.

if(krez /=0) then
write(3,70) t0, u0
endif

! Coeficienti:
dt2=dt*dt
a0 =1./dt2
a1 =1./(2*dt)
a2 =2*a0

! M^:
m1 = a0*m +a1*c
k1 = k-a2*m
k1ant =a0*m -a1*c

! X(t-dt):
Xant =u0 -dt*u10 +u20*dt2/2	! Prima iteratie: Xant =X => Xant
X =u0

! Calculeaza raspunsul:
ave_erru =0.D0
ave_erru1 =0.D0

errU_max =0.
errV_max =0.

numar_pasi =nint(tt/dt)

ipas=0
100	ipas=ipas+1
if(ipas >numar_pasi) goto 500

t1= ipas*dt
!	if(t1.gt.tt+dt/2) goto 500

! Membrul II = R^:
!r1 = p(t1) -(k-a2*m)*X -(a0*m -a1*c)*Xant
r1 = p(t1) -k1*X -k1ant*Xant

! x(t+dt):
X1 =r1/m1
! x'(t):
Xprime =(X1-Xant)*a1

! Erori:
errU = ufct(t1)-X
errV = vfct(t1)-Xprime

ave_erru =ave_erru +errU*errU
errU_max =dmax1(errU_max, dabs(errU))

ERV = vfct(t1)-Xprime
ave_erru1 =ave_erru1 +errV*errV
errV_max =dmax1(errV_max, dabs(errV))

IF(npas ==1 .or. mod(ipas, npas) ==0) THEN
write(2,70) t1, X, Xprime, errU, errV
ENDIF

if(krez /=0) then
	if(npas ==1 .or. mod(ipas, npas_dep) ==0) then
	write(3,70) t1, X
	endif
endif

Xant =X
X = X1
GOTO 100


500	continue

! Write results at tEnd:
write(1,'(/, a)') ' Rezultate la tEnd:'
write(1,60)
if(true_solution) then
write(1,70) t1, Xant, Xprime, errU, errV
else
write(1,70) t1, Xant, Xprime
endif

if(true_solution) then

ave_erru =dsqrt(ave_erru/numar_pasi)
ave_erru1 =dsqrt(ave_erru1/numar_pasi)

write(1,*)
write(1, '(a, E11.3)') 'Eroarea medie patratica in U =', ave_erru
write(1, '(a, E11.3)') 'Eroarea medie patratica in V =', ave_erru1

write(1,*)
write(1, '(a, E11.3)') 'Eroarea maxima in U (in modul) =', errU_max
write(1, '(a, E11.3)') 'Eroarea maxima in V (in modul) =', errV_max

endif

l1 =lstrlen(fname1); l2 =len_trim(fname2)

write(*, '(a)') ' Fisiere de iesire:'
write(*, '(1x, a)') fname1(1:l1)
write(*, '(1x, a)') fname2(1:l2)

close(1)
close(2)

if(krez /=0) then
l3 =len_trim(fname3)
write(*, '(1x, a)') fname3(1:l3)
close(3)
endif

END

