function bublina_tlak_pole()
dbstop if error
clear all; close all; clc;
addpath(genpath('./'))
filename='5um_01f_240kpa.mat'; %kam se uloží výsledky
uloz=1; % uložit výsledky (0/1)
xRP=1; % RP = Rayleigh-Plesset
xMBM=1; % MBM = mikrobublina s membránou

xO=1; % zobrazit jednotlivá øešení
%% ZADÁNÍ PARAMETRÙ %%
    %%POÈÁTEÈNÍ POLOMÌR BUBLINY
    R_0=3e-6; %[m] %R_0=R_20=vnìjší polomìr pro bublinu s membránou
    
    %%BUZENÍ
    fb=1; %budící frekvence (Pokud fb<=10: f=fb*f0, jinak: f=fb) [Hz];
    A=1200e3; %amplituda buzení [Pa]  
    
    %%ÈAS VÝPOÈTU
    t_0=0;
    t_max=1e-5; %[s]

    %%TEPLOTA A TLAK
    T_0=309.65; %poèáteèní teplota [K] telesna teplota
    p_inf0=117321; %tlak kapaliny [Pa] 80mmHg-111989 120mmHg-117321

    %%Další parametry: MBM
    tlm=5e-9; %tlouš�ka membrány [m]  OK
    rho_S=1210; %hustota membrány [kg/m^3]    PEG 2000
    G_S=5e6/(2*(1+0.5)); %modul pružnosti ve smyku membrány [Pa] OK
    G_L=0; %modul pružnosti ve smyku okolí [Pa]  OK
    mu_S=0.00005*rho_S; %dyn. viskozita membrány [Pa.s] OK
    sigma2=0.056; %povrchové napìtí na R_2 [N/m] OK
    sigma1=0.04; %povrchové napìtí na R_1 [N/m]  OK
    R20=R_0; %vnìjší polomìr [m]
    R10=R20-tlm; %vnitøní polomìr [m]
    V_S=R20^3-R10^3;
    %----------------------------------------%
    p_v=IAPWS_IF97('psat_T',T_0)*1e6; %tlak sytých par [Pa]
    sigma=0,05; %povrchové napìtí (0.073)[N/m] povrchové napìtí krve   OK
    p_g0=p_inf0+2*sigma/R_0-p_v; %poèáteèní parciální tlak plynu uvnitø bubliny [Pa]  OK
    mu=0.003; %dynamická viskozita krve [Pa.s] 0,003-0,0036
    kappa=1.4; % Poissonova konstanta [-]
    rho=1060; %hustota krve [kg/m^3]
    c_inf=1581; %rychlost zvuku v kapalinì [m.s^-1]
    B=304.9e6; n=7.15; %konstanty Taitovy stavové rovnice [Pa,-]
    b=0.1735; %van der Waalsova konstanta hexanu [l/mol]
    V_m=22.4; %molární objem [l/mol]
    p_g0_MBM=(p_inf0+2*sigma1/R10+2*sigma2/R20)*(1/(1-b/V_m))^(kappa);
    
    %%BUZENÍ    
    if fb<=10
        f=fb*nat_freq();
    else
        f=fb;
    end
    function p=p_inf(t) 
        p=p_inf0+A*sin(2*pi*f*t); %budicí funkce [Pa]
    end
    function p=dp_inf(t) 
        p=A*cos(2*pi*f*t)*2*pi*f; %derivace budicí funkce [Pa]
    end
%% Vlastní frekvence
    function f_0=nat_freq() %vlastní frekvence bubliny
        f_0=1/(2*pi*R_0)*(1/rho*(3*kappa*p_g0-2*sigma/R_0))^0.5;
    end
    nf=nat_freq();
    
    function f_0=nat_freq_MBM() %vlastní frekvence MBM
    pg0mbm=(p_inf0+2*sigma1/R10+2*sigma2/R20)*(1/(1-b/V_m))^(kappa);
    f_0=1/(2*pi*R20)*(1/rho*(3*kappa*pg0mbm+12*G_S*tlm/R20))^0.5;
    end
    nf_MBM=nat_freq();

%% Mechanický index    
    M=(A*10^-6)/sqrt(nat_freq_MBM*10^-6)
%% VÝPIS
fprintf(['PARAMETRY:\n',...
        '-----------\n',...
        'Poèáteèní teplota:\t%0.2f K \n',...
        'Poèáteèní tlak krve:\t%0.3f Pa \n',...
        'Poèáteèní tlak plynu uvnitø bubliny:\t%0.3f Pa \n\n',...
        'Tlak sytých par:\t%0.3f Pa \n',...
        'Povrchové napìtí:\t%0.3f mN/m \n',...
        'Dynamická viskozita krve:\t%0.3e Pa.s \n',...
        'Hustota krve:\t%0.3f kg/m^3 \n\n',...
        'Modul pružnosti membrány:\t%0.3e N/m\n',...
        'Viskozita membrány:\t%0.3e N/m\n\n'],...
        T_0,p_inf0,p_g0,p_v,sigma*1000,mu,rho,G_S,mu_S)
    
    fprintf('Vlastní frekvence bubliny: %0.6e Hz\n',nat_freq())
    fprintf('Vlastní frekvence MBM: %0.6e Hz\n',nat_freq_MBM())
    fprintf('Budicí frekvence: %0.6e Hz\n\n',f)
    fprintf('Mechanický index: %0.6e \n\n',M)
 
 
%% VOLÁNÍ FUNKCE PRO VÝPOÈET %%
% POÈÁTEÈNÍ PODMÍNKY %
init_R=R_0; % poèáteèní polomìr
init_dRdt=0; % poèáteèní rychlost

reltol=1e-9; % relativní tolerance
abstol=1e-12; % absolutní tolerance
maxdt=1e-6; % maximální èasový krok
options = odeset('RelTol',reltol,'AbsTol',abstol,'Stats','on','MaxStep',maxdt);

if xRP==1 
display('Výpoèet RP:');
tic
[t_RP,Rb_RP]=ode113( @rhs_RP, [t_0 t_max], [init_R init_dRdt], options);
R_RP=Rb_RP(:,1);
dR_RP=Rb_RP(:,2);
toc
display(' ');
end

if xMBM==1
display('Výpoèet MBM:');
tic
[t_MBM,Rb_MBM]=ode113( @rhs_MBM, [t_0 t_max], [init_R init_dRdt], options);
R_MBM=Rb_MBM(:,1);
dR_MBM=Rb_MBM(:,2);
toc
display(' ');
end

fprintf('EQ \t max(R/R0)\t min(R/R0)\t max(dR)\t max(T) \t max(p)\n')
%% VÝPOÈTY RP %%
if xRP==1
for i=1:(length(t_RP)-1)
    dt_RP(i)=t_RP(i+1)-t_RP(i); % zjištìní èasového kroku
end
dt_RP(i+1)=dt_RP(i);

for i=1:length(t_RP)
p_prubeh_RP(i)=p_inf(t_RP(i)); % prùbìh budícího tlaku
end

for i=1:length(R_RP)
p_i_RP(i)=p_g0*(R_0/R_RP(i))^(3*kappa); %tlak uvnitø bubliny
end

for i=1:length(R_RP)
T_i_RP(i)=(R_0/R_RP(i))^(3*kappa-3)*T_0; %teplota uvnitø bubliny [K]
end
fprintf('RP \t %0.3e \t %0.3e \t %0.3e \t %0.3e \t %0.3e\n',max(R_RP)/R_0,min(R_RP)/R_0,max(abs(dR_RP))/R_0,max(T_i_RP),max(p_i_RP))
end
%% VÝPOÈTY MBM %%
if xMBM==1
for i=1:(length(t_MBM)-1)
    dt_MBM(i)=t_MBM(i+1)-t_MBM(i); % zjištìní èasového kroku
end
dt_MBM(i+1)=dt_MBM(i);

for i=1:length(t_MBM)
p_prubeh_MBM(i)=p_inf(t_MBM(i)); % prùbìh budícího tlaku
end

for i=1:length(R_MBM)
p_i_MBM(i)=p_g0_MBM*(R_0/R_MBM(i))^(3*kappa); %tlak uvnitø MBM
end
R1_MBM=(R_MBM.^3-V_S).^(1/3);
for i=1:length(R_MBM)
T_i_MBM(i)=T_0/((R1_MBM(i)/R10)^3-b/V_m)^(kappa-1); %teplota uvnitø MBM [K]
end
fprintf('MBM \t %0.3e \t %0.3e \t %0.3e \t %0.3e \t %0.3e\n',max(R_MBM)/R_0,min(R_MBM)/R_0,max(dR_MBM)/R_0,max(T_i_MBM),max(p_i_MBM))
end

%% VÝSTUP %%
lnwdth=0.8; %tlouš�ka èar
co = {[0 0.4470 0.7410];[0.8500 0.3250 0.0980];...
    [0.9290 0.6940 0.1250];[0.4940 0.1840 0.5560];...
    [0.4660 0.6740 0.1880];[0.3010 0.7450 0.9330];...
    [0.6350 0.0780 0.1840]};

if xRP==1 && xO==1
figure(1)
subplot(3,1,1)
plot(t_RP,R_RP/R_0,'Color',co{1});
hold on

title({'RP',sprintf('R_0=%0.3e m   f=%0.3e Hz   p_A=%0.3e Pa   f/f_0=%0.4f',R_0,f,A,f/nat_freq())})
xlabel('t [s]');ylabel('R/R_0 [-]')
grid on

subplot(3,1,2)
plot(t_RP,p_prubeh_RP,'Color',co{1})
xlabel('t [s]');ylabel('p_i_n_f [Pa]')
grid on

subplot(3,1,3)
plot(t_RP,abs(dR_RP),'Color',co{1})
axis([0 inf -100 1.5*c_inf])
xlabel('t [s]');ylabel('|dR/dt| [m.s^-^1]')
grid on

end

if xMBM==1 && xO==1
figure(2)
subplot(3,1,1)
plot(t_MBM,R_MBM/R_0,'Color',co{1});
title({'MBM',sprintf('R_0=%0.3e m   f=%0.3e Hz   p_A=%0.3e Pa   f/f_0=%0.4f',R_0,f,A,f/nat_freq_MBM())})
xlabel('t [s]');ylabel('R/R_0 [-]')
grid on

subplot(3,1,2)
plot(t_MBM,p_prubeh_MBM,'Color',co{1})
xlabel('t [s]');ylabel('p_i_n_f [Pa]')
grid on

subplot(3,1,3)
plot(t_MBM,abs(dR_MBM),'Color',co{1})
axis([0 inf -100 1.5*c_inf])
xlabel('t [s]');ylabel('|dR/dt| [m.s^-^1]')
grid on
end

%%Porovnání MODELÙ
figure(3)
hold on
pli=1;
PL=struct;

L{pli}='Rayleighùv èas';

if xRP==1
    pli=pli+1;
    PL.pli=plot(t_RP,R_RP/R_0,'LineWidth',lnwdth,'Color',co{pli-1});
    L{pli}='RP';
end
if xMBM==1
    pli=pli+1;
    PL.pli=plot(t_MBM,R_MBM/R_0,'LineWidth',lnwdth,'Color',co{pli-1});
    L{pli}='MBM';
end

xlabel('t [s]')
ylabel('R/R_0 [-]')
title(sprintf('R_0=%3.0f\\mum,  f/f_0=%3.2e,  A=%3.2e kPa',R_0*1e6,f/nat_freq(),A/1e3))
set(findall(gcf,'-property','FontSize'),'FontSize',12)
%% FUNKCE %%
%%RAYLEIGH-PLESSET%%
function dRdt=rhs_RP(t,R) 
dRdt_1 = R(2);
dRdt_2 = 1/(rho*R(1))*(p_v-p_inf(t)+p_g0*(R_0/R(1))^(3*kappa)-2*sigma/R(1)-4*mu*R(2)/R(1))-3/2*(R(2))^2/R(1);
dRdt=[dRdt_1; dRdt_2];
end

%%MIKROBUBLINA S MEMBRÁNOU%%
function dRdt=rhs_MBM(t,R)
R2=R(1);
dR2=R(2);
R1=(R2^3-V_S)^(1/3);

pgt=(p_inf0+2*sigma1/R10+2*sigma2/R20)*(1/((R1/R10)^3-b/V_m))^(kappa);
    
dRdt_1 = R(2);
dRdt_2 = -(dR2^2*((3*rho*(dR2/(3*c_inf) - 1))/2 + rho_S*(dR2/c_inf + 1)*(R2^4/(2*R1^4) - (2*R2)/R1 + 3/2)) - (dR2/c_inf + 1)*(p_inf0 + p_inf(t) - pgt + (2*sigma1)/R1 + (2*sigma2)/R2 - (4*G_L*(R20^3/R2^3 - 1))/3 + (4*dR2*mu)/R2 + (4*G_S*V_S*(R20^3/R2^3 - 1))/(3*(V_S - R2^3)) - (4*V_S*dR2*mu_S)/(R2*(V_S - R2^3))) + (3*R2^3*dR2*kappa*pgt)/(c_inf*R10^3*(b/V_m - R1^3/R10^3)))/(R2*(rho*(dR2/c_inf - 1) - rho_S*(R2/R1 - 1)*(dR2/c_inf + 1)));

dRdt=[dRdt_1; dRdt_2];
end

if uloz==1
    save(filename)
end
end