! Exact solution:
! x =
! y =
! x'=
! y'=
! 
!
subroutine exact(t)		!, uExact, uPrimeExact)
use work

implicit real(8) (a-h, o-z)
!real(8) xExact(ngl), xPrimeExact(ngl)


!Define exact solution:
uExact(1) = ufct(t)


uPrimeExact(1)= vfct(t)


end subroutine exact


!----------------------------------
real(8) function ufct(t)
use work

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


	if(zeta.ne.0) goto 200

!raspuns neam. liber sau fortat
ufct0 = C20*dcos(omega*t) +C10*dsin(omega*t)
	if(omp.eq.0.) then
	ufct=ufct0
	return
	else
ufct = ufct0 + ust*DA*(dsin(omp*t) -omp/omega*dsin(omega*t))
	endif
	return

!raspuns liber am. cu zeta<=> 1
200	if(zeta ==1) then
ufct=dexp(-omega*t)*(C1*t +C2)
return

else if(zeta <1) then
ufct=dexp(-zeta*omega*t)*( C1*dsin(oma*t) +C2*dcos(oma*t))

else
ufct=dexp(-zeta*omega*t)*( C11*dexp(-oma*t) +C21*dexp(oma*t))
endif

return
end

!----------------------------
real(8) function vfct(t)
use work

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


if(zeta.ne.0) goto 200

!raspuns neam. liber sau fortat
! ufct0= C20*dcos(omega*t) +C10*dsin(omega*t)
vfct0 = -C20*omega*dsin(omega*t) +C10*omega*dcos(omega*t)

if(omp.eq.0.) then
vfct=vfct0

else
!ufct= ufct0 + ust*DA*(dsin(omp*t) -omp/omega*dsin(omega*t))
vfct = vfct0 + ust*DA*omp*(dcos(omp*t) -dcos(omega*t))
endif

RETURN

!raspuns liber am. cu zeta<=> 1
200	if(zeta ==1) then
!fct=dexp(-omega*t)*(C1*t +C2)
vfct=dexp(-omega*t)*C1 -omega*dexp(-omega*t)*(C1*t +C2)

else if(zeta <1) then
ufct=dexp(-zeta*omega*t)*( C1*dsin(oma*t) +C2*dcos(oma*t))
!
vfct=-zeta*omega*ufct +dexp(-zeta*omega*t)*( C1*oma*dcos(oma*t) -C2*oma*dsin(oma*t))

else
ufct=dexp(-zeta*omega*t)*( C11*dexp(-oma*t) +C21*dexp(oma*t))
!
vfct=-zeta*omega*ufct +dexp(-zeta*omega*t)*( -oma*C11*dexp(-oma*t) +oma*C21*dexp(oma*t))
endif

return
end


!-----------------------------------------
subroutine coef_lin()
use work

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


zeta =e

omega =pi2
om2 =pi2*pi2
k =om2

p0 =k
c =2*zeta*omega
omp =0.8*omega

oma =omega*dsqrt(dabs(1-zeta*zeta))

if(omp /=0) then
da =1./(1-omp*omp/om2)
ust =P0/k
endif


c20 =u0(1)
c10 =u10(1)/omega
if(c ==0.) RETURN
c2 =u0(1)

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

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

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

END
