	PROGRAM NEWMARK_1
! 1 singura ecuatie
! Newmark in dubla precizie
! Coeficientii: b, g
!
! Ecuatie:
! - Neliniara, de forma:
! mu'' + ga(u,u') +fk(u) =p(t)
!
! - Liniara:
! mu'' +cu' +ku =p(t)
! ga = c; fk =k;
!
! Febr. 2006
!
! Fisierul de intrare - continut:
!===============================
! 1) Titlul problemei
! 2) Lin: Caracter <=20. Tip ecuatie: "L,l" =Liniara; Alt car.: Neliniar.
! 3) m, zeta, k		! k=0 => k=4*pi^2
! 4) P0, omp		! raspuns fortat: omp<0: => omp =0.8*omega; P0<0: => P0 =omega**2
! 5) dt, tt, npas	! npas =pas de tiparire
! 6) t0, u0, u10	! u10 <0 => u10 =2*pi
! 7) krez			! 0/#0: nu/da - se scrie si fisierul "nume.dep" (date: t;u(t))
!	7-1) krez, npas_dep		! pentru krez #0; npas_dep =pas de tiparire pt. "nume.dep"
! 8) 1/beta, gama	! coeficientii: 0 => beta =1/6; 0 => gama =1/2
! 9) 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 b0, g0/6d0, 0.5d0/
data IND/' Er(U) **'/

character(20) lin

integer width
integer status

!logical kind_1, write_dep


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(1:1) =='L' .or. lin(1:1) =='l'

read(1, *) m, zeta, k
if(k ==0.) then
omega = pi2			! 2*pi
om2 = omega*omega	! 4*(pi)^2
k =m*om2
!*
!k =pi4

else 
om2 =k/m
omega =dsqrt(om2)
endif

c =2*m*omega*zeta	!	zeta=c/(2*m*omega)
oma =omega*dsqrt(dabs(1-zeta*zeta))

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
if(u10 <0.) u10 =pi2

read (1,*) krez
write_dep =krez /=0

if(write_dep) then
backspace 1
read(1,*) krez, npas_dep

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

! beta, gamma:
read(1, *) b, g
if(b ==0.) b =b0
if(g ==0.) g =g0

b =1._8/b

read(1,*) true_solution


! Coef. pt. Solutia exacta - liniar:
!if(liniar) then

	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
!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)
30 format(/, ' tetaP = ',G14.8,'; P0 = ',G14.8)

if(omp ==0. .and. P0 ==0.) then
write(1,'(/, a)') ' Raspuns liber'
else
write(1,30) omp,P0
endif

iw =width(npas)
!40	format(/, ' dt =', f10.4,';  tt =', f10.2, :, ';  npas = ', I<iw>)
40	format(/, ' dt =', 1pG10.2,';  tt =', G11.5, :, ';  npas = ', I<iw>)

write(1,40) dt, tt, npas

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

50 format(' Coeficienti:'/, ' beta = 1/', G9.3, 4x, 'gama =', G11.4, /)
write(1,50) 1_8/b, g

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(write_dep) then
write(3,70) t0, u0
endif

dt2=dt*dt

ave_erru =0.D0
ave_erru1 =0.D0

errU_max =0.D0
errV_max =0.D0

numar_pasi =nint( (tt-t0)/dt)
	
tol =1D-10; lnit =20;
call Newmark(TOL, lnit, ave_erru, ave_erru1, errU_max, errV_max, t_errU_max, t_errV_max)



if(true_solution) then

write(1, '(/, a, t11, 2(1pG21.11))') 'Exact:', ufct(tt), vfct(tt)

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


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


write(1, '(/, a, 1pG11.3, a, g12.6)') 'Eroarea maxima in U (in modul) =', errU_max, ', la t = ', t_errU_max
write(1, '(a, 1pG11.3, a, g12.6)') 'Eroarea maxima in V (in modul) =', errV_max, ', la t = ', t_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(write_dep) then
l3 =len_trim(fname3)
write(*, '(1x, a)') fname3(1:l3)
close(3)
endif

END



