function [Omega0_SDOF,Theta_SDOF_vec,e31_SDOF_vec,cPVDFmat,sPVDFmat,Theta_expression_SDOF]= couplingSDOF(L,Le,Lp,Offset,elem,be,b,hs,hp,Es,Ep,mi_p,rho_s,rho_p,mM,Rz,psi,accel,Cp,f_mer,V_mer,usek)
Ls=L-usek;                                              % rozdelenie na nonsika hmotu
mM = mM + ((usek)*b*hs*rho_s);                          % prepocitana hmotnost zavazia s castou nosnika
elemE = [round(Offset/elem) round((Le+Offset)/elem)];   % cislo prvku z L, kde zacina a konci elektroda
sPVDFmat = [    1/Ep -mi_p/Ep  -mi_p/Ep 0 0 0;          % matica poddajnosti
                -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);                                 % matica tuhosti
%________________________________________________________________________________________________________
%% 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
Jp=((b*hp^3)/12);

%________________________________________________________________________________________________________
%% EKVIVALENTNA HMOTNOST
ms =rho_s*as*Ls;    %substrat
mp =rho_p*ap*Lp;    %PVDF
M = ((33/140)*(ms+mp))+mM; %ekvivalentna hmotnost

%% EKVIVALENTNA TUHOST
ks=3*Es*Js/(Ls^3);  %tuhost substrat
kp=3*Ep*Jp/(Lp^3);  %tuhost PVDF
K = ks+kp;          %ekvivalentna tuhost

%% VLASTNA FREKVENCIA
Omega0 = sqrt(K/M);
f_0 = sqrt(K/M)/(2*pi);
fprintf('f0_SDOF %.4f Hz \n',f_0)
Omega0_SDOF = Omega0;


%%  NORMOVANY VLASTNY TVAR
fun = @(LAM) ((cos(LAM)*cosh(LAM))+1);
x0 = 0.1;
LAM=abs(fzero(fun,x0));
syms FIshape(x)
FIshape(x) = ((cosh(LAM*x) - cos(LAM*x))-(((1/sqrt(M))*(sinh(LAM) - sin(LAM))/(cosh(LAM) + cos(LAM)))*(sinh(LAM*x) - sin(LAM*x))));

oldnorm=double(int(power(FIshape,2),x,0,L));
FIn = FIshape/sqrt(oldnorm)*sqrt(1/M);  %normovany vlastny tvar priehybu
FIn_d = diff(FIn,x);                    %normovany vlastny tvar natocenia

%% ELEKTROMECHANICKY COUPLING
Theta_expression_SDOF = be*(zTp-zn)* ...
    (double(FIn_d(elemE(2)*elem))-double(FIn_d(elemE(1)*elem))) ...
    /(double(FIn(L))); 

omega_b=2*pi.*f_mer;
e31_sdof_vect=zeros(numel(V_mer),1);
Theta_sdof_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_sdof = roots(th);
        Theta_sdof_vect(j,1) = Theta_sdof;
        e31_sdof_vect(j,1) = -Theta_sdof/Theta_expression_SDOF;
end

[~,idx_mean] = max(V_mer);
e31_SDOF_vec=-e31_sdof_vect(idx_mean-2:idx_mean+2);
Theta_SDOF_vec=-Theta_sdof_vect(idx_mean-2:idx_mean+2);

end