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

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)



zeta =e

! ngl =5
n0 =6

! Set ki:
r0 =1.2D5		! kN/m

do igl =1,ngl

ki(igl) =r0 +(8D0 -igl)*1D3*n0		! kN/m
enddo

! Define KK:
KK =0.

!kk(1, 1:2) =(/ki(1), -ki(1)/)

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

!kk(3, 2:4) =(/-ki(2), ki(3)+ki(2), -ki(3)/)

!kk(4, 3:5) =(/-ki(3), ki(4)+ki(3), -ki(4)/)

!kk(5, 4:5) =(/-ki(4), ki(5)+ki(4)/)

do igl =1, ngl

if(igl ==1) then
kk(1, 1:2) =(/ki(1), -ki(1)/)

elseif(igl ==ngl) then

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

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

enddo


!KK =KK*1D3		! N

! 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(1) =250_8; pp(2:5) =0.;	! N (25 kN); ngl =5

omp =14_8		!omp =1.1*om(1)

p =pp*dsin(omp*t)


END subroutine FGP_J
