! Calculeaza:
! Functii: f(u); g(u,u'); p(t);
! Jacobieni: jf = df/du; jg = dg/du; jg1 =dg/du1;
! Derivata: dp;
!
! Cadru etajat - 6 grade de lib. Excitatie armonica.
! ngl =6

subroutine fgp_j(t, uu, v)
use work

implicit real(8) (a-h, o-z)
real(8) uu(ngl), v(ngl)

real(8) KK(ngl,ngl)		!, PP(ngl)
real(8) ki(ngl), r0
real(8) NO

zeta =e

! Define KK:
KK =0.

r0 =1.2D5	! kN/m
NO =10

!
do i=1, ngl
ki(i) =r0 +(8-i)*NO*1D3
enddo

! Gen. K = tridiagonala;
! 1
kk(1, 1:2) =(/ki(1), -ki(1)/)

! ngl
kk(ngl, ngl-1:ngl) =(/-ki(ngl-1), ki(ngl)+ki(ngl-1)/)

! 2:(ngl-1)
do igl =2, ngl-1

kk(igl, igl-1:igl+1) =(/-ki(igl-1), ki(igl)+ki(igl-1), -ki(igl)/)
enddo


! f(u) =K*U
f =matmul(KK,UU)


! Jacobianul lui f:
JF =KK

! g(u,u'):
g =0.

! Jacobianul Dg/Du:
jg =0.

! Jacobianul Dg/Du':
jg1 =0.


! p(t):
pp =0D0;
pp(1) =250_8;

omp =4.8679456016374711		!omp =1.1*om(1)

p =pp*dsin(omp*t)

! P derivat:
dp =omp*pp*cos(omp*t)

END subroutine FGP_J
