! Exact solution:
subroutine exact(t)
use work

implicit real(8) (a-h, o-z)
real(8) y(ngl), dy(ngl)
real(8) m_fi(ngl), fi_y(ngl), fi_dy(ngl)		! _ = *; fi(ngl, ngl)
real(8) mm, d, y_st
real(8) pp0(ngl)

character(80) titlu



! Eigenvalues
! 1   13.89554166359515    
! 2   4.727337193155136    
! 3   32.13050484969981    
! 4   21.91639594497070    
! 5   28.15698096897788    

! Eigenvectors (of R)

! Defineste Constante:
Const_If: IF(t ==tEnd) THEN

!om(1) =  13.89554166359515
!om(2) =  4.727337193155136    
!om(3) =  32.13050484969981    
!om(4) = 21.91639594497070
!om(5) = 28.15698096897788

! Omega:
open (2, file ='L6-OMEGA.dat', status ='old', form ='formatted')
do it=1,3
read(2,'(a)') titlu
enddo

do igl =1, ngl

read(2,*) i, om(i)
enddo

close(2)

! Vect. proprii:
open (2, file ='L6-VEC_X.dat', status ='old', form ='formatted')
do it=1,3
read(2,'(a)') titlu
enddo

do igl =1, ngl
read(2,*) j		! vector #j

	do i =1, ngl
	read(2,*) fi(i,j)
	enddo
enddo

close(2)


!pp(1) =250D3
pp0 =matmul(transpose(fi), pp)

aa =0.; m_fi =0.;

do igl =1,ngl
d =1./( 1-(omp/om(igl))**2 )

m_fi =matmul(M, fi(:, igl))
mm =dot_product(fi(:,igl), m_fi)


y_st =pp0(igl)/( mm*om(igl)**2)

aa(igl) =y_st*d
enddo


!*
write(11,*) 'Omega:', om

write(11, '(/,a)') ' FI:'
do i =1,ngl
write(11,*); write(11,*) fi(i,:);
enddo

write(11,'(/,a)') 'PP0:'
write(11,*) pp0

write(11,'(/,a)') ' A:'
write(11,*) aa
!*

ENDIF Const_If


! Solutia:
do igl =1,ngl
y(igl)  =aa(igl)*dsin(omp*t) -omp*aa(igl)/om(igl)*dsin(om(igl)*t)
dy(igl) =omp*aa(igl)*( dcos(omp*t) -dcos(om(igl)*t) )
enddo

uExact = matmul(fi, y)

uPrimeExact = matmul(fi, dy)


end subroutine exact
