function [Omega0ndof,Theta_NDOF_vect,e31_NDOF_vect,cPVDFmat,sPVDFmat,Theta_expression_NDOF]= couplingNDOF(L,Le,Lz,Lp,Offset,elem,be,b,hs,hp,Es,Ep,mi_p,rho_s,rho_p,mM,Rz,psi,accel,Cp,f_mer,V_mer)
N=fix(L/elem);
nnel=2;                 %pocet dof na element
uzoldof=2;              % pocet dof na uzol
nnode=(nnel-1)*N+1;     % celkovy pocet nodov
Ndof=nnode*uzoldof;     % celkovy pocet supnov volnosti 

rho_z = mM/(Lz*b*hs);

%prvky
elemE = [round(Offset/elem) round((Le+Offset)/elem)];%cislo prvku z L, kde zacina a konci elektroda
elemP = [round(Offset/elem) round((Lp+Offset)/elem)]; %cislo prvku z L, kde zacina a konci piezo
elemZ = [round((L-Lz)/elem) N]; %cislo prvku z L, kde zacina a konci zavazie

sPVDFmat = [    1/Ep -mi_p/Ep  -mi_p/Ep 0 0 0;
                -mi_p/Ep 1/Ep -mi_p/Ep 0 0 0;
                -mi_p/Ep -mi_p/Ep 1/Ep 0 0 0;
                0 0 0 (2*(1+mi_p))/Ep 0 0;
                0 0 0 0 (2*(1+mi_p))/Ep 0;
                0 0 0 0 0 (2*(1+mi_p))/Ep];
cPVDFmat=inv(sPVDFmat);

%________________________________________________________________________________________________________
%% NEUTRALNA OS
SS = (hs+hp)/2; %geom. stred
zTs = -(SS-(hs/2)); %vzd. tazisko substratu -geom. stred
zTp = hs+(hp/2)-SS; %vzd. tazisko pieza -geom. stred
zn = ((Es*zTs*hs)+(Ep*zTp*hp))/((Es*hs)+(Ep*hp)); %neutralna os

% PRIEREZY
as=b*hs; % substrat
ap=b*hp; % PVDF

%________________________________________________________________________________________________________
%% osove kvadraticke momenty 
Js=((b*hs^3)/12); % substrat

%________________________________________________________________________________________________________
%% HMOTNOST
%HMOTNOST elementu
me=[    156 22*elem 54 -13*elem;...
        22*elem 4*elem^2 13*elem -3*elem^2;...
        54 13*elem 156 -22*elem;...
        -13*elem -3*elem^2 -22*elem 4*elem^2];
ms =((rho_s*as*elem)/420)*me; %substrat
mps =( (((rho_p*ap)+(rho_s*as) )*elem)/420)*me;%PVDF  
mz =(((rho_z+rho_s)*as*elem)/420)*me; %zavazie

% GLOBALNA HMOTNOST
M_G=zeros(Ndof,Ndof);

for k = 1:N
    if (k <= elemP(1) )   ||   ((k > elemP(2)) && (k <= elemZ(1)) )
      M_G((k-1)*2+1:(k-1)*2+4,(k-1)*2+1:(k-1)*2+4)=M_G((k-1)*2+1:(k-1)*2+4,(k-1)*2+1:(k-1)*2+4)+ms;
    elseif ( (k > elemP(1)) && (k <= elemP(2) )  )
        M_G((k-1)*2+1:(k-1)*2+4,(k-1)*2+1:(k-1)*2+4)=M_G((k-1)*2+1:(k-1)*2+4,(k-1)*2+1:(k-1)*2+4)+mps;
    elseif((k > elemZ(1)) && (k <= elemZ(2)) )
        M_G((k-1)*2+1:(k-1)*2+4,(k-1)*2+1:(k-1)*2+4)=M_G((k-1)*2+1:(k-1)*2+4,(k-1)*2+1:(k-1)*2+4)+mz;
    end
end
M=M_G; % Globalna hmotnost 

%________________________________________________________________________________________________________
%% TUHOST
%laminatova teoria - konstanty
konst_A = (Es*hs)+(Ep*hp);
konst_B = 0.5*hs*hp*(Ep-Es);
konst_D = (1/12)*(    (  Es*((hs^3)+(3*hs*hp^2))  ) + (  Ep*((hp^3)+(3*hp*hs^2))  )      );
EspJsp = ((konst_D-(konst_B^2/konst_A))*b);

% TUHOST ELEMENTU
ke=[    12 6*elem -12 6*elem;...
        6*elem 4*elem^2 -6*elem 2*elem^2;...
        -12 -6*elem 12 -6*elem;...
         6*elem 2*elem^2 -6*elem 4*elem^2];  
ks=((Es*Js)/elem^3)*ke;% substrat 
kps=(EspJsp/elem^3)*ke; %PVDF

% GLOBALNA TUHOST
K_G=zeros(Ndof,Ndof);
for k = 1:N
     if ( (k > elemP(1)) && (k <= elemP(2) )  )
        K_G((k-1)*2+1:(k-1)*2+4,(k-1)*2+1:(k-1)*2+4) = K_G((k-1)*2+1:(k-1)*2+4,(k-1)*2+1:(k-1)*2+4)+kps;
     else
        K_G((k-1)*2+1:(k-1)*2+4,(k-1)*2+1:(k-1)*2+4) = K_G((k-1)*2+1:(k-1)*2+4,(k-1)*2+1:(k-1)*2+4)+ks;
     end
end
K=K_G; %globalna tuhost

%%
OP=3;
M=M(OP:end,OP:end); %okr. podmienka
K=K(OP:end,OP:end); %okr. podmienka

%% VLASTNA FREKVENCIA
Omega0 = min(abs(sqrt(eig(K,M))));
Omega0ndof = Omega0;
f_0 = Omega0/(2*pi);
fprintf('f0_NDOF %.4f Hz \n',f_0(1))

%% VLASTNY TVAR
[Eigenvectors,Eigenvalues] = eig(K,M);
[~,ind]=sort(diag(Eigenvalues),'ascend');
Modeshapes_sorted = Eigenvectors(:,ind);
FI = Modeshapes_sorted(:,1);
FI = FI/(sqrt(FI'*M*FI));

%________________________________________________________________________________________________________
%% COUPLING
phi = FI(1:2:2*N); %neparne-posuvy
phi_d = FI(2:2:2*N); %parne - natocenia

Theta_expression_NDOF= be* (zTp-zn)* (phi_d(elemE(2))- phi_d(elemE(1))) /phi(elemZ(2));
omega_b=2*pi.*f_mer;

e31_ndof_vect=zeros(numel(V_mer),1);
Theta_ndof_vect=zeros(numel(V_mer),1);
%kvadraticka rovnica
for j=1:numel(V_mer)
        menovatel = ((-omega_b(j)^2)+(2*psi*i*omega_b(j)*Omega0)+Omega0^2);
        Koef_A = (V_mer(j)*i*omega_b(j));
        Koef_B = +(i*omega_b(j)*accel);
        Koef_C = ((V_mer(j)/Rz)+(V_mer(j)*Cp*i*omega_b(j)))*menovatel;
        discr=sqrt(Koef_B^2 - (4*Koef_A*Koef_C));
        roots = zeros(2,1);
        roots(1)= ( -Koef_B + discr ) / (2*Koef_A);
        roots(2)= ( -Koef_B - discr ) / (2*Koef_A);
        [~,th]=min(abs(roots));
        Theta_ndof = roots(th);
        e31_ndof_vect(j,1) = -Theta_ndof/Theta_expression_NDOF;
        Theta_ndof_vect(j,1) = Theta_ndof;
end
[~,idx_mean] = max(V_mer);
e31_NDOF_vect=-e31_ndof_vect(idx_mean-2:idx_mean+2);
Theta_NDOF_vect=-Theta_ndof_vect(idx_mean-2:idx_mean+2);

end