PROGRAM NewOP_1
! 1 singura ecuatie
! NewOp in Dubla Precizie
! Improved Newmark cu sau fara mediere pt. U3 - pondere "niu"
! Coeficientii limita: p,teta, p1,teta1, p2,teta2
! Eroarea + Eroarea medie patratica - in Deplasare si Viteza.
!
! 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) kodm			! cod mediere: NU: =0; DA: #0
!	5-1) kodm, niu	! Daca kodm #0; niu : factor de mediere.
! 6) dt, tt, npas	! npas =pas de tiparire
! 7) t0, u0, u10	! u10 <0 => u10 =2*pi
! 8) krez			! 0/#0: pentru scriere fisier deplasari
! 9) true_solution	! logic
! ========================================

use dfwin
use work

implicit real(8) (a-h,o-z)
real(8) m1

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

logical mediere, 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

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

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

read(1, *) P0, omp

if(P0 <0)  then
P0 =om2
endif

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

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

teta =teta0
b = 1./b0
write(*,*) b, teta
bt =b*teta
	
teta1 =1./teta10
g = g0
write(*, *) g, teta1
gt =g*teta1

teta2 =teta20
d = d0
write(*,*) d, teta2
dtt =d*teta2

read(1,*) kodm
mediere =kodm /=0

if(mediere) then
backspace 1
read (1,*) kodm, niu
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

read(1,*) true_solution


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
u30 = (dp(t0)-c*u20 -k*u10)/m

else
! mu'' + ga(u,u') +fk(u) =p(t):
u20 = (p(t0) -ga(u0,u10) -fk(u0))/m
u30 = (dp(t0)-dga(u0,u10)*u10 -dga1(u0,u10)*u20 - dfk(u0)*u10)/m

endif

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


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

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

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

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,30) omp, p0
endif

40	format(/, ' dt = ', f10.4,'; tt =', f10.2, '; npas = ', i4)
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

!
write(1, '(a, 1pg15.7, a, 1pg15.7, /)') ' u20 =', u20, 'u30 =', u30
!

50 format(' Coeficienti:'/ &
'  beta =',F10.6,';  teta =',F10.6/ &
'  gama =',F10.6,'; teta1 =',F10.6/ &
' delta =',F10.6,'; teta2 =',F10.6//)
write(1,50) b,teta, g,teta1, d,teta2

51 format(' *** Mediere:   niu =',f5.3 / )
if(mediere) write(1,51) niu

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

dt2=dt*dt
dt3=dt2*dt

ave_erru =0.D0
ave_erru1 =0.D0

errU_max =0.
errV_max =0.

numar_pasi =nint( (tt -t0)/dt)
	
u  =u0
u1 =u10
u2 =u20
u3 =u30

lnit =20
TOL =1D-12

MAIN_DO: DO ipas= 1, numar_pasi
	t1= t0 +ipas*dt

	u0  =u
	u10 =u1
	u20 =u2
	u30 =u3

	u  =u0  +u10*dt +u20*dt2/2. +u30*b*dt3
	u1 =u10 +u20*dt +u30*g*dt2
	u2 =u20 +u30*d*dt
	u3 =u30

IF(liniar) THEN
	pp  =p(t1) -m*u2 -c*u1 -k*u
	m1  =m*dtt*dt +c*gt*dt2 +k*bt*dt3
	du3 = pp/m1

ELSE	! solve by Newton's method
	! iteration:
	w =0_8

	newton_do: do iter =1, lnit

! mu'' + ga(u,u') +fk(u) =p(t):
	x   = u + bt*dt3*w
	x1 = u1 + gt*dt2*w
	x2 = u2 + dtt*dt*w        

	f  =m*x2 +fk(x) +ga(x,x1) -p(t1)

	df =m*dtt*dt +bt*dt3*dfk(x) +bt*dt3*dga(x,x1) +gt*dt2*dga1(x,x1)

	w =w -f/df
		
	if( dabs(f/df) <= TOL) then
	exit newton_do
	endif
	
	enddo newton_do

77	du3 =w
!if(ipas ==1) print*, 'iter_1 =', iter

ENDIF


	u  = u +bt*dt3*du3
	u1 =u1 +gt*dt2*du3
    u2 =u2 +dtt*dt*du3

	u3 =u3 +du3


!Mediere
Med_If: IF(mediere) THEN

Lin_IF: if(liniar) then
u31 =(dp(t1) -c*u2 -k*u1)/m

else
u31 = (dp(t1)-dga(u,u1)*u1 -dga1(u,u1)*u2 - dfk(u)*u1)/m

endif Lin_If

u3  =niu*u3 + (1-niu)*u31


ENDIF Med_If


if(true_solution) then
ERU = ufct(t1)-u
ave_erru =ave_erru +eru*eru
errU_max =dmax1(errU_max, dabs(eru))

ERV = vfct(t1)-u1
ave_erru1 =ave_erru1 +erv*erv
errV_max =dmax1(errV_max, dabs(erv))
endif

IF(npas ==1 .or. mod(ipas,npas) ==0) THEN
if(true_solution) then
write(2,70) t1, u, u1, eru, erv

else
write(2,70) t1, u, u1
endif
ENDIF

if(write_dep) then
	if(npas ==1 .or. mod(ipas, npas_dep) ==0) then
	write(3,70) t1, u
	endif
endif


ENDDO MAIN_DO


! Write results at tEnd:
write(1,'(/, a)') ' Rezultate la tEnd:'
write(1,60)

if(true_solution) then
write(1,70) t1, u, u1, eru, erv
!*
write(1,'(a, 2(1pG21.11))') '     Exact:', ufct(t1), vfct(t1)
!*

else
write(1,70) t1, u, u1
endif


if(true_solution) then
write(1, '(/, a, 1pg13.3)') 'Eroarea medie patratica in U =', ave_erru
write(1, '(a, 1pg13.3)') 'Eroarea medie patratica in V =', ave_erru1


write(1, '(/, a, 1pg13.3)') 'Eroarea maxima in U (in modul) =', errU_max
write(1, '(a, 1pg13.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(write_dep) then
l3 =len_trim(fname3)
write(*, '(1x, a)') fname3(1:l3)
close(3)
endif

END

