clear all; close all; clc;
set(groot,'defaulttextinterpreter','latex');  
set(groot, 'defaultAxesTickLabelInterpreter','latex');  
set(groot, 'defaultLegendInterpreter','latex');
elem = 5e-4;
g=9.81;
%% vstupne parametre
% DLZKY
L = 56e-3;      % cely substrat [m]
Lp= 43e-3;      % PVDF [m]
Le = 43e-3;     % elektroda [m]
Lz = 6e-3;      % zavazie [m]

usek=3.4e-3;
% usek=4.5e-3;
Ls = L-usek;    % substrat bez zavazia [m] 
OffsetE = 1e-3; % vzdialenost elektrody od votknutia[m]
OffsetP=OffsetE;% vzdialenost PVDF od votknutia[m]

%PRIRADENIE CISLA ELEMENTU DLZKE
N=L/elem;       % pocet elementov
elemE = [round(OffsetE/elem) round((Le+OffsetE)/elem)]; % cislo prvku z L, kde zacina a konci elektroda
elemP = [round(OffsetP/elem) round((Lp+OffsetP)/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

%SIRKY
b = 10e-3;      % nosnik[m]
be = 7.5e-3;    % elektroda[m]

%HRUBKY
hs=.3e-3;       % substrat [m]
hp=.13e-3;      % PVDF[m]

%Youngove moduly
Es=2.1e11;      % substrat[Pa]
Ep=2e9;         % PVDF[Pa]

%PERMITIVITY
eps_0 = 8.854187e-12;   % vakuum [F/m]
eps33_t=12;             % konstantne napatie [-]
eps33_s=11.75272;       % konstantne pretvorenie [-]

%HUSTOTY
rho_s=7850;     % substrat[kg/m^3]
rho_p=1760;     % PVDF [kg/m^3]

%ZAVAZIE 

mM =5.3e-3 + ((usek)*b*hs*rho_s);
%PIEZOELEKTRICKE KONSTANTY
d31 = -1.9e-11;     % [C/N]
e31 = -5.6903e-2;   % [C/m]

%ODPOR ZATAZE
Rz = 1e7;           % [Ohm]

%TLMENIE
xi = 0.01;        % [-]

%ZRYCHLENIE
accel = 50e-3*g; 

%SIMULACIA CASOVA ODOZVA
t_end = 0.2;        % cas simulacie
U_0 = -2e-3;        % pociatocna vychylka


%% Prierezy, osove kvadraticke momenty
%PRIEREZY
as=b*hs;            % substrat
ap=b*hp;            % PVDF 

%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

%OSOVE KVADRATICKE MOMENTY 
Js=(b*hs^3)/12;     % substrat  
Jp=(b*hp^3)/12;     % PVDF

%OHYBOVA TUHOST LAMINAT
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); %ohybova tuhost

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

%% 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('Vlastna frekvencia sustredne parametre sqrt(K/M) %.4f Hz \n',f_0)

%% TLMENIE
C = 2*xi*M*sqrt(K/M);

%% KAPACITA
Cp = (Le*be*eps33_s*eps_0)/hp;

%%  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,Ls));
FIn = FIshape/sqrt(oldnorm)*sqrt(1/M);  % normovany vlastny tvar priehybu
FIn_d = diff(FIn,x);                    % normovany vlastny tvar natocenia

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

%__________________________________________________________________________
%% AMPLITUDOVO FREKVENCNA CHARAKTERISTIKA
freq_start=20; 
freq_end=21.5; 
Om_start=2*pi*freq_start; 
Om_end = freq_end*2*pi;
d_Om=0.01;

U_sdof=zeros(1,numel(Om_start:d_Om:Om_end));
V_sdof=zeros(1,numel(Om_start:d_Om:Om_end));
j=0;
for omega=Om_start:d_Om:Om_end 
j=j+1;
    konst = ((Cp*i*omega)+(1/Rz));
    U_sdof(j) = (-accel*konst)/((((-omega^2)+(2*xi*i*omega*Omega0)+Omega0^2)*konst)+((Theta^2)*i*omega));
    V_sdof(j)= (Theta*i*omega*accel) / (  ((-omega^2+(2*xi*i*omega*Omega0)+Omega0^2)*konst) +((Theta^2)*i*omega) );
   
end
Omega0_sdof = Omega0;
save('SDOFMatlabPlotAF','U_sdof','V_sdof','Omega0_sdof');

%__________________________________________________________________________
%% VOLNE KMITANIE V CASOVEJ OBLASTI
t_end = 2;
U0=-2e-3;
out = sim('SDOF_model_TR.slx',t_end);
t_sdof = out.ScopeData.time;
U_sdof = out.ScopeData.signals(1).values;
V_sdof = out.ScopeData.signals(2).values;
save ('SDOFMatlabPlotTR','t_sdof','U_sdof','V_sdof')








