close all;
clear all;
clc;

% PARAMETRY MEØÍCÍHO PØÍSTROJE
Dm=1.22*10^-1; % N*s/m, tlumení díky obtekaní ve snímaèi
ks=370.87; % N/m, tuhost pružiny ve snímaèi
kf=1.0286*10^6; % N/m, tuhost rámu
m=523.5*10^-6; % kg, hmotnost indentoru

Ei=1140*10^9; % Pa, modul pružnosti indentoru
vi=0.07; % Poisonovo èíslo indentoru

%PARAMETRY Z EXPERIMENTU
v=0.46; % Poisonovo cislo vzorku

%Naètení dat z excelu
filename='data_VZ1.xlsx'; sheet = 1; 

%Aplituda posunutí
u00Range='Q4:Q58';
u00=xlsread(filename,sheet,u00Range)*10^-9; %m
u0=u00(end:-1:1); %(srovnání podle frekvence od nejmenší po nejvìtší)

%Silova amplituda
F00Range='U4:U58';
F00=xlsread(filename,sheet,F00Range)*10^-6; %N
F0=F00(end:-1:1);

%Kontaktní plocha
Ac0Range='BI4:BI58';
Ac0=xlsread(filename,sheet,Ac0Range)*10^-18; % m^2
Ac=Ac0(end:-1:1);

%Tuhost indentoru
for i=1:length(Ac);
    ki(i)=((2*Ei)/(1-vi^2))*sqrt(Ac(i)/pi); % N/m
end

%Frekvence
f0Range = 'O4:O58';
f0=xlsread(filename,sheet,f0Range); 
f=f0(end:-1:1); %Hz

omega=2*pi*f; % rad/s

%Fázové úhly
fi0Range='S4:S58';
fi0=abs(xlsread(filename,sheet,fi0Range)); % °
fi=fi0(end:-1:1);

%dynamické tuhosti
for i=1:length(F0);
    DS(i)=F0(i)/u0(i); %N/m
end
%% 
%MINIMÁLNÍ, MAXIMÁLNÍ A STØEDNÍ HODNOTY
omega_min=omega(1); % rad/s
omega_max=omega(length(omega)); % rad/s
omega_str=omega(length(omega)/2+0.5); % rad/s

fi_min=fi(1); % °
fi_max=fi(length(fi)); % °
fi_str=fi(length(fi)/2+0.5); % °

DS_min=DS(1); % N/m
DS_max=DS(length(DS)); % N/m

%URÈENÍ P A Q
%f(omega_min)
fce_min=1/(DS_min*(1/(omega_min*cscd(fi_min)))-Dm);

%f(omega_max)
fce_max=1/(DS_max*(1/(omega_max*cscd(fi_max)))-Dm);

O=[1 omega_min^2;
    1 omega_max^2];
b=[fce_min;fce_max];
x=O^-1*b;
P=x(1,1);
Q=x(2,1);

%URÈENÍ ALPHA
syms alpha

Dc=((1+alpha)^2)/P;
k1=(1+alpha)/sqrt(P*Q);
k2=alpha*k1;

A=((k1+ks)/(k1+k2))*Dc+Dm;
B=Dc/(k1+k2);
G=((k1*k2)/(k1+k2))+ks;
H=(Dc*Dm)/(k1+k2);

g=vpa((omega_str*(A-B*G)+omega_str^3*B*H)/(G+omega_str^2*A*B-omega_str^2*(H+m)-omega_str^4*m*B^2)-tand(fi_str)); %

alpha=double(solve(g));

%URÈENÍ Dc, k1, k2, kc
Dc=((1+alpha)^2)/P % N*s/m, tlumení ve vzduchu
k1=(1+alpha)/sqrt(P*Q) % N/m, kombinovaná tuhost rámu, indentoru a kontaktní tuhosti
k2=alpha*k1 % N/m, tuhost viscoelastického elementu

for i=1:length(ki);
    kc(i)=1/(1/k1-1/kf-1/ki(i)); % N/m, kontaktni tuhost  
end

kc_prum=sum(kc)/length(kc);
kc=kc_prum

%STORAGE A LOSS MODUL
for j=1:length(omega);
    k_storage(j)=(kc*k2*(kc+k2)+omega(j)^2*Dc^2*kc)/((kc+k2)^2+omega(j)^2*Dc^2);
    E_storage(j)=(k_storage(j)/2)*sqrt(pi/Ac(j))*(1-v^2);
    
    k_loss(j)=(omega(j)*Dc*kc^2)/((kc+k2)^2+omega(j)^2*Dc^2);
    E_loss(j)=(k_loss(j)/2)*sqrt(pi/Ac(j))*(1-v^2);
end
%% 
%GRAFY
%Storage modul
figure
plot(f,(E_storage/10^6),'+ m')
xlabel('Frekvence [Hz]')
ylabel('Storage modul [MPa]')

%Loss modul
figure
plot(f,(E_loss/10^6),'+ m')
xlabel('Frekvence [Hz]')
ylabel('Loss modul [MPa]')

