PROGRAM NewOp_ngl
! "ngl" ecuatii (ord. 2)
! Dubla Precizie
! Cu sau fara mediere pt. U3 - pondere "niu"
! Coeficientii limita: p,teta, p1,teta1, p2,teta2
! Eroarea + Eroarea medie patratica - in Deplasare si Viteza.
!
! Ecuatii:
! - Neliniara, de forma:
! MU'' + G(U,U') +F(U) =P(t)
! Functiile G, F, P; jacobienii JF, JG; derivata DP: in FGP_J
!
! - Liniara:
! MU'' +CU' +KU =P(t)
! G = C; F =K;
!
! Febr. 2006
!
! Fisierul de intrare - continut:
!===============================
! 1) Titlul problemei
! 2) ngl
! 3) e	! parametru
! 4) t0, h, tEnd	! 
! 5) u0, u10	! 
! 6) kodm			! cod mediere: NU: =0; DA: #0
!	6-1) kodm, niu	! Daca kodm #0; niu : factor de mediere.
! 7) k_rez			! cod rezultate: 0/#0: pentru scriere fisier deplasari & viteze.
!	7-1) tStart_H, tEnd_H, ipas_H ! pentru k_rez #0
! 8) true_solution	! logic
! ========================================

use dfwin
use work

use DFLIB

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

character(max_path) fname1
character(max_path), allocatable :: fname(:)
character(1) igl_char

real(4) s

!real(8) niu

common /coef/	beta_h3, gama_h2, delta_h, &
				teta_beta_h3, teta1_gama_h2, teta2_delta_h

data teta , beta  /.25_8, 6._8/       ! beta = 1/beta
data teta1, gama  / 3._8,  .5_8/        ! teta1= 1/teta1
data teta2, delta /.5_8,  1._8/

logical PrintIt
logical true_solution

real(8) norma_2

real(8), allocatable :: b(:)

real(8) C10,C20,C1,C2,C11,C21,da,ust
common/coef/C10,C20,C1,C2,C11,C21,da,ust

integer width


! coeficientii limita:
beta = 1d0/beta
teta1= 1d0/teta1
! delta = 1d0


call GetFile(fname1, status)

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

! Read data; title.
read(1,'(a)') title

read(1,*) ngl


! Alocari:
allocate( u0(ngl), u10(ngl), u20(ngl), u30(ngl) )

allocate( uEnd(ngl), u1End(ngl), u2End(ngl), u3End(ngl) )

allocate( m(ngl,ngl) )
allocate( f(ngl), g(ngl) )
allocate( jf(ngl, ngl), jg(ngl, ngl), jg1(ngl, ngl) )

allocate( p(ngl), dp(ngl) )

allocate( b(ngl) )
allocate( pp(ngl) )


read(1,*) e
read(1,*) t0, h, tEnd	!
read(1,*) u0, u10

read(1,*) eps, lnit		!, PrintIt

read(1,*)  kodm
	if(kodm /=0) read(1,*) niu

read(1,*) k_rez	! 0, #0 - for writing time-history: X, Xprime
	if(k_rez /=0) then
	read(1,*) tStart_H, tEnd_H, iPas_H
	endif

read(1,*) true_solution

!*
if(ngl ==1) then
C10 =0.; C20 =0.; C1 =0.; C2 =0.; C11 =0.; C21 =0.; da =0.; ust =0.;
call coef_lin()
endif
!*

! Coeficienti:
h2 = h*h; h3 = h2*h
beta_teta =beta*teta;	gama_teta1 =gama*teta1;	delta_teta2 =delta*teta2
beta_h3 =beta*h3; gama_h2= gama*h2; delta_h =delta*h
teta_beta_h3 = teta*beta_h3; teta1_gama_h2 =teta1*gama_h2; teta2_delta_h =teta2*delta_h


! Inital conditions:
! u0(ngl), u1(ngl)
!
! Calculeaza U20, U30:
! MU'' +G(U,U') +F(U) =P(t):  
!
! MU''' +JG(U,U')U' +JG1(U,U')U'' +JF(U)U' =DP(t)

call Gen_M()
                                              
call FGP_J(t0, u0, u10)	! Cu t0, u0, u10, etc.

! M*U20 = P - G - F
b = P - G - F
kod =0
call ELIM (ngl, M, b, prag, kod, kodPrint)	!
U20 =b

!*
write(1,*); write(1,*) 'U20:', u20;

! M*U30 = dP0 - jG*U10 - jG1*U20 - jF*U10
DO l=1,NGL
b1 =0d0
do j=1,ngl
b1 = b1 + jG(l,j)*U10(j) + jG1(l,j)*U20(j) + jF(l,j)*U10(j)
enddo

b(l) =dP(l) -b1
ENDDO

kod =0
call Gen_M()

call ELIM (ngl, M, b, prag, kod, kodPrint)	!
U30 =b

!*
write(1,*) 'U30:', u30

No_Intervals =(tEnd-t0)/h +0.5

allocate (u(0:No_Intervals, ngl), uPrime(0:No_Intervals, ngl))

! Initialize:
u(0,:) =u0; uPrime(0,:) =u10
u(1:No_Intervals,:) =0.0; uPrime(1:No_Intervals,:) =0.0


! Open for time-history:
IF(k_rez /=0) THEN
allocate(fname(1:ngl))

iPoint =index(fname1,'.',back=.TRUE.)

do igl =1, ngl		! ngl <=9
ic =inttochar(igl, igl_char, 10)

fname(igl) =fname1(1:iPoint-1)//'-'//igl_char//'.dat'

iunit =igl +1
open(iunit, file=fname(igl), status='unknown')

enddo


ENDIF


print *, Title
print *, 'ngl:', ngl
print *, 't0, Step:', t0,h
print *, 'Initial cond.(system):',u0, u10
print *, 'No.of Intervals:', No_Intervals

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


PrintIt =.false.

! Start Time:
s = secnds(0.0)
! Time-history integration:
call gen_M()	! pentru Newop -> sub_ab

NewOp_DO: DO i =1, No_Intervals
t =t0+ i*h

call NewOp(t, PrintIt, kod_elim)

if(kod_elim <0) exit NewOp_DO

u(i,:) =uEnd; uPrime(i,:) =u1End
u0 =uEnd; u10 =u1End; u20 =u2End; u30 =u3End

ENDDO NewOp_DO

s = secnds(s)

if(kod_elim <0) then

write(1,'(/, a)') ' Pivot < prag. Solutia nu se poate obtine.'
goto 500
endif

iw =width(No_Intervals)
write(1, '(/, a, I<iw>)') ' No.of Intervals: ', No_Intervals

if(kodm.ne.0) write(1, '(/, a, F6.2)') '  *** Mediere: niu =', niu

print '(/, a, f10.2, a)',' Integration Time =', s,' sec.'

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

write(1, '(/, a, g12.5)') ' tEnd:', t;

write(1,'(/,a, <ngl>(1pG21.11))') ' End: U : ', u(No_intervals, :)
write(1,'(a, <ngl>(1pG21.11))')   '      U'': ', uPrime(No_intervals, :)



IF(k_rez /=0) THEN
iStart =(tStart_H -t0)/h
iEnd =(tEnd_H -t0)/h

print '(/,a)',' Writing history ...'

IGL_DO: DO igl =1,ngl
iunit =igl+1

His_Do: DO i =iStart, iEnd, iPas_H

write(iunit, *) i*h, u(i, igl)

enddo His_Do

close(iunit)

ENDDO IGL_DO

!close(2); close(3); close(4); ;close(5)

ENDIF


deallocate( u0, u10, u20, u30 )
!deallocate( m )
deallocate( f, g )
deallocate( jf, jg, jg1 )
deallocate( p, dp )
deallocate( b )

! Now, calculate the error at tEnd and Maximumm error:

True_If: IF(true_solution) THEN
allocate( uExact(ngl), uError(ngl), uPrimeExact(ngl), uPrimeError(ngl) )

allocate(om(ngl), aa(ngl), fi(ngl,ngl))

! Error:
call exact(tEnd)		! uExact, uPrimeExact

!*
if(ngl <=6) then
write(1,'(/, a, <ngl>(1pG21.11))') ' U_exact :', uExact(:)
write(1,'(a, <ngl>(1pG21.11))') ' U''_exact:', uPrimeExact(:)
endif
!*

uError =uExact -uEnd
uPrimeError =uPrimeExact -u1End

write(1, '(/, a)') ' Error at tEnd:'
write(1,'(a, <ngl>(2x,1pD10.3))') ' U :', uError
write(1,'(a, <ngl>(2x,1pD10.3))') ' U'':', uPrimeError

! Maximum Error:
uError_max =0D0
uPrimeError_max =0D0

I_DO: DO i =1, No_Intervals
t =t0 +i*h
call exact(t)		! uExact, uPrimeExact

! Error:
uError =uExact -u(i,:)
uPrimeError =uPrimeExact -uPrime(i,:)

uNorm =norma_2(uError)

if(uNorm > uError_max) then
uError_max =uNorm
t_max =t
endif

uPrimeNorm =norma_2(uPrimeError)

if(uPrimeNorm > uPrimeError_max) then
uPrimeError_max =uPrimeNorm
t1_max =t
endif

END DO I_DO

write(1, '(/, a)') ' Maximum Error Norm in U and U'''
write(1,'(a, f10.6, a, 2x,1pD10.3)') ' U:', t_max,  ':', uError_max
write(1,'(a, f10.6, a, 2x,1pD10.3)') ' U''',t1_max, ':', uPrimeError_max

ENDIF True_If

500 continue
close (1)
deallocate( m )
deallocate(u, uPrime)
deallocate( uEnd, u1End, u2End, u3End )
deallocate( pp )

if(true_solution) then
deallocate( uExact, uError, uPrimeExact, uPrimeError )

deallocate(om, aa)
endif

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

end program NewOp_ngl
