clc, clear all,close all, warning off
%% Simulation 
t=1 ;%time of simulation [s] 
Rz=250;% external resistive load 
F=106.6; % Síla na amp
F_in=2*1364.3;% Síla na stack
delta_x=4e-4;
Q=30; % Quality factor   
% Amplifier
k_amp=F_in/delta_x;
ome_amp=364 ;% vlastni frekvence amplifieru V Hz ?!?!?
m_amp=(F/delta_x)/(ome_amp^2); % efektivni hmotnost Amplifieru
b_amp=2*0.02*m_amp*ome_amp;% tlumení amplifieru...odhad
%b_amp=1000;% tlumení amplifieru přepis aby to dobře vypadalo

%%Terfenol-D
ro_terf=9250 ;% Hustota terfenolu-D [Kg/m^3]
E=50e+10 ;% Yang?v modul  nebo tam mám dát modul v tahu což je n?jakých 30e+6 ??????????????????????????? 
d33 =2e-9; % magnetostrikční konstanta Terfenolu

% Rozměry terfenolu-D
d_terf=0.008 ;% vnitrní prumer terfenolového válce
l_terf=0.025 ;% Délka terfenolového válce
A=(pi*d_terf^2)/4 ;% prurez válcem terfenolu-D(stejný jako piezo)

%% Parametry cívky 
%počet závitů se přepočítá aby byl rovnoměrně rozložen po celé délce 
N_c=4000; % počet závitů (zhruba)
d_d=1e-4 ;% prumer drátu cívky 
ro_Cu=0.0169e-6; % rezistivita mědi [ohm*m]

l_c=l_terf; % délka cívky == delka terfenolu
d_c=d_terf; % vnitrní prumer cívky == prumer terfenolu
N_c_row=floor(l_c/d_d); % počet závitů v jednom řádku
N_c=N_c_row*floor(N_c/N_c_row);% přepočítaný počet závitů
D_c=((N_c/N_c_row)*2*d_d)+d_c; % vnější prumer cívky
d_c_stredni=(d_c+D_c)/2; % střední průměr závitu
A_c=((pi*D_c^2)/4)-((pi*d_c^2)/4) ;% prurez cívky
l_d=d_c_stredni*pi*N_c;% celková délka drátu[m]
Rc=4*(ro_Cu*l_d)/(pi*d_d^2); % Elektrický odpor cívky [Ohm]
% Indukčnost
perm=6.29e-6; % permeabilita terfenolu..realně není konstanta ale mění se podle B/H
perm_0=pi*4e-7 ;% permeabilita
L=(perm*(N_c^2)*A_c)/l_c; % indukcnost 

%% Mechanické vlastnosti Terfenolu-D
m_t=ro_terf*A*l_terf; % hmotnost terfenolu-D
k=E*A/l_terf; % tuhost terfenolu-D
ome=sqrt(k/m_t);
b_p=1/(2*Q);% pom?rný útlum terfenolu-D
b=2*b_p*m_t*ome; % tlumení terfenolu-D
%G=(N_c*d33*E*A)/l_c % Koueficient převodu !!???jak to má být spravně????
m=m_t+m_amp; %celková efektivní hmotnost 
%% Polynom závislosti megnetického toku na mach. napětí
fi=[211 152 40.5 34.5 8.85 2.3]*2e-8; % Hodnoty magnetického toku [190 140 38 32 8 2.1] z FEMM
napeti=[55 48 34 20 3 0]*1e+6; % Mechanické napětí v MPa
p=polyfit(napeti,fi,3);
%% Slábnutí mag pole se vzdáleností
r=d_c_stredni/2;

%% Vydrží terfenol ?
napeti_max=700e+6;% Pevnost v tlaku [pa] 
%[pa] % http://www.matweb.com/search/datasheet.aspx?matguid=aa68cad05c7c4d39932834d68e669a5d&ckck=1
format long
napeti=F_in/A;
if napeti<napeti_max
    fprintf('Terfenol vydrží')
    %sim const_model3.slx
    sim magnetorstikcni_generator_model.slx
else
    fprintf('Dojde k mechanickému poškození Terfenolu')
end
