clear all; close all; clc; g=9.81; eps_0=8.854e-12;
set(groot,'defaulttextinterpreter','latex');  
set(groot, 'defaultAxesTickLabelInterpreter','latex');  
set(groot, 'defaultLegendInterpreter','latex');
set(0,'defaultAxesFontSize',20)
set(0, 'DefaultLineLineWidth', 4);
set(0, 'DefaultAxesLineWidth', 1);
set(0,'DefaultLineMarkerSize',5);
modra = [0.3010, 0.7450, 0.9330];
zelena = [0.4660, 0.6740, 0.1880];
cervena = [0.8350, 0.1, 0.0];

%_______________________________________________________________________________________________________________________________________________________________
%% DATA Z MERANIA - NAPATIE,FREKVENCIA,RYCHLOST,TLMENIE
%% meranie 4 =50mg
accel=50e-3*g;  % Amplituda zrychlenia zakladne [m/s^2]
load('Vf_meranie4vz2.mat');

%% ulozenie merania do premennych
f_mer = vars(:,1);      % Merana frekvencia [Hz]
V_mer = vars(:,2);      % Merane napatie [V]
Udot_mer = vars(:,3);   % Merana rychlost [m/s]
load('KSI_vz2.mat')
fprintf('Pomerny utlm %.4f \n',xi);
%_______________________________________________________________________________________________________________________________________________________________
%% VSTUPNE PARAMETRE 
L=71.9e-3;      % celkova dlzka nosnika [m]
Le=48.2e-3;     % dlzka elektrody [m]
Lp=63e-3;       % dlzka pieza [m]
Lz=4e-3;        % dlzka konca, kde je zavazie [m]
be=8e-3;        % sirka elektrody [m]
elem = 0.1e-3;  % velkost elementu, plati L/elem==integer [m]
Offset = 1e-3;  % vzdialenost od voknutia [m]
mM = 2.6e-3;    % hmotnost zavazia %+( paska vazi cca 0.2g !)[m]
Rz=1e7;         % Odporova zataz
%SDOF
usek=2.92e-3;   % dlzka so zavazim

%_______________________________________________________________________________________________________________________________________________________________
%% DALSIE PARAMETRE v nasom pripade vacsinou rovnake parametre
b=10e-3;        % sirka nosnika
hs=.3e-3;       % hrubka substratu
hp=.13e-3;      % hrubka PVDF
Es=2.1e11;      % Youngov modul substrat
Ep=2e9;         % Youngov modul PVDF
mi_p=0.34;      % Poissonovo cislo PVDF
rho_s=7850;     % hustota substratu [kg/m^3]
rho_p=1760;     % hustota PVDF [kg/m^3]
eps33 = 12;     % relativna permitivita

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%_____VYPOCTY____%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%_______________________________________________________________________________________________________________________________________________________________
%% HODNOTY ROVNAKE PRE OBA MODELY
Cp=(Le*be*eps33*eps_0)/hp ; % kapacita
eps33t=eps33*eps_0;         % permitivita pri konstantnom napati

%_______________________________________________________________________________________________________________________________________________________________
%% MATLAB NDOF VYPOCET COUPLINGU Z MERANIA VF
%VSTUPY  = merane napatie,frekvencia
%VYSTUPY = Omega0, Theta, e31
[Omega0_NDOF,Theta_NDOF_vec,e31_NDOF_vec,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,xi,accel,Cp,f_mer,V_mer);

% PIEZOELEKTRICKE KONSTANTY
Theta_NDOF = mean(Theta_NDOF_vec);
e31_NDOF = mean(e31_NDOF_vec);
d31_NDOF = (sPVDFmat(1,1)+sPVDFmat(1,2))*e31_NDOF;
k31_NDOF = abs(sqrt((d31_NDOF^2*(cPVDFmat(1,1)+cPVDFmat(2,1)))/eps33t));

% REPRODUKCIA NAPATIA PODLA THETA Z MERANIA
d_Om=0.01;
V_NDOF=zeros(1,numel(Omega0_NDOF*0.9:d_Om:Omega0_NDOF*1.1)); U_NDOF=zeros(1,numel(Omega0_NDOF*0.9:d_Om:Omega0_NDOF*1.1));
j=0;
for omega_b= Omega0_NDOF*0.9:d_Om:Omega0_NDOF*1.1  %budiaca frekvencia
    j=j+1;
    clen_NDOF = ((Cp*i*omega_b)+(1/Rz));
%   U_NDOF(j) = (-accel*clen)/((((-omega_b^2)+(2*psi*i*omega_b*Omega0_NDOF)+Omega0_NDOF^2)*clen)+((Theta_NDOF^2)*i*omega_b));
    V_NDOF(j)= (Theta_NDOF*i*omega_b*accel) / (  ((-omega_b^2+(2*xi*i*omega_b*Omega0_NDOF)+Omega0_NDOF^2)*clen_NDOF) +((Theta_NDOF^2)*i*omega_b) );
end


%_______________________________________________________________________________________________________________________________________________________________
%% MATLAB SDOF VYPOCET COUPLINGU Z MERANIA VF
%VSTUPY  = merane napatie,frekvencia, variabilna dlzka nosnika
%VYSTUPY = Omega0, Theta, e31
[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,xi,accel,Cp,f_mer,V_mer,usek);

% PIEZOELEKTRICKE KONSTANTY
Theta_SDOF = mean(Theta_SDOF_vec);
e31_SDOF = mean(e31_SDOF_vec);
d31_SDOF = (sPVDFmat(1,1)+sPVDFmat(1,2))*e31_SDOF;
k31_SDOF = abs(sqrt((d31_SDOF^2*(cPVDFmat(1,1)+cPVDFmat(2,1)))/eps33t));

d_Om=0.01;
V_SDOF=zeros(1,numel(Omega0_SDOF*0.9:d_Om:Omega0_SDOF*1.1 ));
j=0;
for omega_b= Omega0_SDOF*0.9:d_Om:Omega0_SDOF*1.1 
    j=j+1;
    clen_SDOF = ((Cp*i*omega_b)+(1/Rz));
%     U_SDOF(j)=(-accel*konst)/((((-omega_b^2)+(2*psi*i*omega_b*Omega0_SDOF)+Omega0_SDOF^2)*clen_SDOF)+((Theta_SDOF^2)*i*omega_b));%
    V_SDOF(j)= (Theta_SDOF*i*omega_b*accel) / (  ((-omega_b^2+(2*xi*i*omega_b*Omega0_SDOF)+Omega0_SDOF^2)*clen_SDOF) +((Theta_SDOF^2)*i*omega_b) );
end

%_______________________________________________________________________________________________________________________________________________________________
%% PLOT POROVNANIE MERANEHO A REPRODUKOVANEHO NAPATIA Z VYPOCITANEJ THETA VF
figure('units','normalized','outerposition',[0 0 1 1])
    [~,idx_mean] = max(V_mer);
    plot((0.9*Omega0_NDOF:d_Om:1.1*Omega0_NDOF) /(2*pi),(abs(V_NDOF)),'-b','Color',modra','MarkerFaceColor',modra); %Matlab nDOF
    hold on
    plot((0.9*Omega0_SDOF:d_Om:1.1*Omega0_SDOF) /(2*pi),(abs(V_SDOF)),'--','Color',cervena','MarkerFaceColor',cervena); %Matlab SDOF
    plot(f_mer,V_mer,'-ob','Color',zelena','MarkerFaceColor',zelena) %meranie
    grid on
    set(gca,'box','off')
    xlabel('f [Hz]');
    ylabel('V [V]');
    title('Vzorka 2')
    legend('Matlab NDOF','Matlab SDOF', 'Meranie','EdgeColor','none');
    set(gca,'XMinorTick','on')
fprintf('____PARAMETRE Z NAPATIA A FREKVENCIE = KVADRATICKA ROVNICA_____\n');
fprintf('Theta_SDOF     %.3fE-7 []\n'   ,abs(Theta_SDOF)*1e6)
fprintf('Theta_NDOF     %.3fE-7 []\n'   ,abs(Theta_NDOF)*1e6)    
fprintf('e31_SDOF       %.4f [C/m]\n'   ,-abs(e31_SDOF))
fprintf('e31_NDOF       %.4f [C/m]\n'   ,-abs(e31_NDOF))
fprintf('d31_SDOF       %.4f [pC/N]\n'  ,-abs(d31_SDOF)*1e12)
fprintf('d31_NDOF       %.4f [pC/N]\n'  ,-abs(d31_NDOF)*1e12)
fprintf('k31_SDOF       %.4f [%%] \n'   ,k31_SDOF*100)
fprintf('k31_NDOF       %.4f [%%] \n'   ,k31_NDOF*100)
fprintf('V_SDOF_max     %.2f [V] \n'    ,max(abs(V_SDOF)))
fprintf('V_NDOF_max     %.2f [V]\n'     ,max(abs(V_NDOF)))
fprintf('V_mer_max     %.2f [V]\n'      ,max(V_mer))
fprintf('f_0_mer        %.2f [Hz]\n'      ,f_mer(idx_mean))
 
% saveas(gca,'C:\Users\simon\OneDrive - Vysoké učení technické v Brně\VUT_Dokumenty\DP\LaTeX\DP_latex\grafy\VfzE31poc_vz2.png')
%_______________________________________________________________________________________________________________________________________________________________
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%_______________________________________________________________________________________________________________________________________________________________
%% VYPOCET Z NAMERANEJ RYCHLOSI
%VSTUPY  = merane napatie-Vmer,frekvencia-f_mer,rychlost-Udot
%VYSTUPY = Theta_mer
omega_b=2*pi*f_mer;
for j = 1:numel(V_mer)
clen = ((Cp*i*omega_b(j))+(1/Rz));
Theta_mer(j) = (V_mer(j)*clen)/Udot_mer(j);
end
Theta_MER=mean(Theta_mer);
%PLOT THETA Z ELEKTRICKEJ ROVNICE
figure('units','normalized','outerposition',[0 0 1 1])
    grid on
    scatter(Udot_mer,abs(Theta_mer),200,'filled','MarkerFaceColor',modra)
    hold on
    xlabel('$\dot U_{mer}$ [m/s]')
    ylabel('$\vartheta$ [C]')
    grid on
    ylim([0 2e-6])
    title('Vzorka 2')
    legend('\rmA_{0}= 50 mg','Location','Northwest','Interpreter','tex','color','white','EdgeColor', get(legend, 'Color' ),'EdgeColor','none')
% saveas(gca,'C:\Users\simon\OneDrive - Vysoké učení technické v Brně\VUT_Dokumenty\DP\LaTeX\DP_latex\grafy\VUdot_vz2.png')

%_______________________________________________________________________________________________________________________________________________________________
%% THETA Z ELEKTRICKEJ ROVNICE DO MATLAB NDOF 
%PIEZO KONSTANTY
%VSTUP: zavisi na modalnom tvare z modelu - Theta_expression_NDOF = be* (zTp-zn)* (phi_d(elemE(2))- phi_d(elemE(1))) /phi(elemZ(2))
e31_NDOF_zUdot = -Theta_MER./Theta_expression_NDOF;
d31_NDOF_zUdot = (sPVDFmat(1,1)+sPVDFmat(1,2))*mean(e31_NDOF_zUdot);
k31_NDOF_zUdot = abs(sqrt((d31_NDOF_zUdot^2*(cPVDFmat(1,1)+cPVDFmat(2,1)))/eps33t));

% REPRODUKCIA NAPATIA Z MERANEJ THETY Z RYCHLOSTI
d_Om=0.01;
V_NDOF_zUdot=zeros(1,numel(Omega0_NDOF*0.9:d_Om:Omega0_NDOF*1.1)); U_NDOF_zUdot=zeros(1,numel(Omega0_NDOF*0.9:d_Om:Omega0_NDOF*1.1));
j=0;
for omega_b= Omega0_NDOF*0.9:d_Om:Omega0_NDOF*1.1 
    j=j+1;
    clen = ((Cp*i*omega_b)+(1/Rz));
%     U_NDOF_zUdot(j) = (-accel*clen)/((((-omega_b^2)+(2*psi*i*omega_b*Omega0_NDOF)+Omega0_NDOF^2)*clen)+((Theta_mer^2)*i*omega_b));
    V_NDOF_zUdot(j)= (mean(Theta_mer)*i*omega_b*accel) / (  ((-omega_b^2+(2*xi*i*omega_b*Omega0_NDOF)+Omega0_NDOF^2)*clen) +((mean(Theta_mer)^2)*i*omega_b) );
end

%_______________________________________________________________________________________________________________________________________________________________
%% THETA Z ELEKTRICKEJ ROVNICE DO MATLAB SDOF 
%PIEZO KONSTANTY
%VSTUP: zavisi na modalnom tvare z modelu - Theta_expression_SDOF = be* (zTp-zn)* (phi_d(elemE(2))- phi_d(elemE(1))) /phi(elemZ(2))
e31_SDOF_zUdot = -Theta_MER./Theta_expression_SDOF;
d31_SDOF_zUdot = (sPVDFmat(1,1)+sPVDFmat(1,2))*e31_SDOF_zUdot;
k31_SDOF_zUdot = abs(sqrt((d31_SDOF_zUdot^2*(cPVDFmat(1,1)+cPVDFmat(2,1)))/eps33t));

% REPRODUKCIA NAPATIA Z MERANEJ THETY Z RYCHLOSTI
d_Om=0.01;
V_SDOF_zUdot=zeros(1,numel(Omega0_SDOF*0.9:d_Om:Omega0_SDOF*1.1)); U_SDOF_zUdot=zeros(1,numel(Omega0_SDOF*0.9:d_Om:Omega0_SDOF*1.1));
j=0;
for omega_b= Omega0_SDOF*0.9:d_Om:Omega0_SDOF*1.1 
    j=j+1;
    clen = ((Cp*i*omega_b)+(1/Rz));
%     U_NDOF_zUdot(j) = (-accel*clen)/((((-omega_b^2)+(2*psi*i*omega_b*Omega0_NDOF)+Omega0_NDOF^2)*clen)+((Theta_mer^2)*i*omega_b));
    V_SDOF_zUdot(j)= (mean(Theta_mer)*i*omega_b*accel) / (  ((-omega_b^2+(2*xi*i*omega_b*Omega0_SDOF)+Omega0_SDOF^2)*clen) +((mean(Theta_mer)^2)*i*omega_b) );
end

% PLOT REPRODUKOVANE NAPATIE
figure('units','normalized','outerposition',[0 0 1 1])
    plot((0.9*Omega0_NDOF:d_Om:1.1*Omega0_NDOF) /(2*pi),(abs(V_NDOF_zUdot)),'-b','Color',modra','MarkerFaceColor',modra); %Matlab nDOF
    hold on
    plot((0.9*Omega0_SDOF:d_Om:1.1*Omega0_SDOF) /(2*pi),(abs(V_SDOF_zUdot)),'--','Color',cervena','MarkerFaceColor',cervena); %Matlab SDOF
    plot(f_mer,V_mer,'-ob','Color',zelena','MarkerFaceColor',zelena) %meranie
    grid on
    set(gca,'box','off')
    xlabel('f [Hz]');
    ylabel('V [V]');
    title('Vzorka 2')
    legend('Matlab NDOF','Matlab SDOF', 'Meranie','EdgeColor','none');
    set(gca, 'XMinorTick','on')
% saveas(gca,'C:\Users\simon\OneDrive - Vysoké učení technické v Brně\VUT_Dokumenty\DP\LaTeX\DP_latex\grafy\VfzE31Udot_vz2.png')
fprintf('____PARAMETRE Z NAPATIA, RYCHLOSTI A FREKVENCIE = ELEKTRICKA ROVNICA_____\n');
fprintf('Theta z Udot   %.2fE-7 []\n'               ,abs(Theta_MER)*1e6)
fprintf('e31_SDOF       %.4f [C/m]\n'               ,-abs(e31_SDOF_zUdot))
fprintf('e31_NDOF       %.4f [C/m]\n'               ,-abs(e31_NDOF_zUdot))
fprintf('d3_SDOF        %.4f  [pC/N]\n'             ,-abs(d31_SDOF_zUdot)*1e12)
fprintf('d31_NDOF       %.4f [pC/N]\n'              ,-abs(d31_NDOF_zUdot)*1e12)
fprintf('k31_SDOF       %.4f [%%] \n'               ,k31_SDOF_zUdot*100)
fprintf('k31_NDOF       %.4f [%%]\n'                ,k31_NDOF_zUdot*100)
fprintf('V_SDOF_max     %.2f [V] \n'                ,max(abs(V_SDOF_zUdot)))
fprintf('V_NDOF_max     %.2f [V] \n'                ,max(abs(V_NDOF_zUdot))) 
fprintf('V_mer_max      %.2f [V]\n'                 ,max(V_mer)) 