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)
% Rovnice, které se budou počítat (0/1):
xRP=1; % RP = Rayleigh-Plesset
xHR=1; % HR = Herring
xGM=1; % GM = Gilmore
xKM=1; % KM = Keller-Miksis
xTS=1; % TS = Tomita-Shima
xMBM=1; % MBM = mikrobublina s membránou

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

    %%TEPLOTA A TLAK
    T_0=293.15; %počáteční teplota [K]
    p_inf0=101325; %tlak kapaliny [Pa]

    %%Další parametry: MBM
    tlm=1.95e-9; %tloušťka membrány [m]
    rho_S=1100; %hustota membrány [kg/m^3]
    G_S=2.8e9/(2*(1+0.33)); %modul pružnosti ve smyku membrány [Pa]
    G_L=0; %modul pružnosti ve smyku okolí [Pa]
    mu_S=1.77; %dyn. viskozita membrány [Pa.s]
    sigma2=0.056; %povrchové napětí na R_2 [N/m]
    sigma1=0.04; %povrchové napětí na R_1 [N/m]
    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=surf_tension(T_0); %povrchové napětí (0.073)[N/m]
    p_g0=p_inf0+2*sigma/R_0-p_v; %počáteční parciální tlak plynu uvnitř bubliny [Pa]
    mu=IAPWS_IF97('mu_pT',p_inf0/1e6,T_0); %dynamická viskozita vody [Pa.s]
    kappa=1.4; % Poissonova konstanta [-]
    rho=1/IAPWS_IF97('v_pT',p_inf0/1e6,T_0); %hustota vody [kg/m^3]
    c_inf=1481; %rychlost zvuku v kapalině [m.s^-1]
    B=304.9e6; n=7.15; %konstanty Taitovy stavové rovnice [Pa,-]
    b=0.1727; %van der Waalsova konstanta [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_MBM();
    
%% Rayleighův čas
    tau=0.915*R_0*(rho/(p_inf0-p_v))^(0.5); % [s]
    
%% VÝPIS
fprintf(['PARAMETRY:\n',...
        '-----------\n',...
        'Počáteční teplota:\t%0.2f K \n',...
        'Počáteční tlak vody:\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 vody:\t%0.3e Pa.s \n',...
        'Hustota vody:\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)
 
%% 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 xHR==1
display('Výpočet HR:');
tic
[t_HR,Rb_HR]=ode113( @rhs_HR, [t_0 t_max], [init_R init_dRdt], options);
R_HR=Rb_HR(:,1);
dR_HR=Rb_HR(:,2);
toc
display(' ');
end

if xGM==1 
display('Výpočet GM:');
tic
[t_GM,Rb_GM]=ode113( @rhs_GM, [t_0 t_max], [init_R init_dRdt], options);
R_GM=Rb_GM(:,1);
dR_GM=Rb_GM(:,2);
toc
display(' ');
end

if xKM==1
display('Výpočet KM:');
tic
[t_KM,Rb_KM]=ode113( @rhs_KM, [t_0 t_max], [init_R init_dRdt], options);
R_KM=Rb_KM(:,1);
dR_KM=Rb_KM(:,2);
toc
display(' ');
end

if xTS==1 
display('Výpočet TS:');
tic
[t_TS,Rb_TS]=ode113( @rhs_TS, [t_0 t_max], [init_R init_dRdt], options);
R_TS=Rb_TS(:,1);
dR_TS=Rb_TS(:,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 HR %%
if xHR==1
for i=1:(length(t_HR)-1)
    dt_HR(i)=t_HR(i+1)-t_HR(i); % zjištění časového kroku
end
dt_HR(i+1)=dt_HR(i);

for i=1:length(t_HR)
p_prubeh_HR(i)=p_inf(t_HR(i)); % průběh budícího tlaku
end

for i=1:length(R_HR)
p_i_HR(i)=p_g0*(R_0/R_HR(i))^(3*kappa); %tlak uvnitř bubliny
end

for i=1:length(R_HR)
T_i_HR(i)=(R_0/R_HR(i))^(3*kappa-3)*T_0; %teplota uvnitř bubliny [K]
end
fprintf('HR \t %0.3e \t %0.3e \t %0.3e \t %0.3e \t %0.3e\n',max(R_HR)/R_0,min(R_HR)/R_0,max(abs(dR_HR))/R_0,max(T_i_HR),max(p_i_HR))
end
%% VÝPOČTY GM %%
if xGM==1
for i=1:(length(t_GM)-1)
    dt_GM(i)=t_GM(i+1)-t_GM(i); % zjištění časového kroku
end
dt_GM(i+1)=dt_GM(i);

for i=1:length(t_GM)
p_prubeh_GM(i)=p_inf(t_GM(i)); % průběh budícího tlaku
end

for i=1:length(R_GM)
p_i_GM(i)=p_g0*(R_0/R_GM(i))^(3*kappa); %tlak uvnitř bubliny
end

for i=1:length(R_GM)
T_i_GM(i)=(R_0/R_GM(i))^(3*kappa-3)*T_0; %teplota uvnitř bubliny [K]
end
fprintf('GM \t %0.3e \t %0.3e \t %0.3e \t %0.3e \t %0.3e\n',max(R_GM)/R_0,min(R_GM)/R_0,max(abs(dR_GM))/R_0,max(T_i_GM),max(p_i_GM))
end
%% VÝPOČTY KM %%
if xKM==1
for i=1:(length(t_KM)-1)
    dt_KM(i)=t_KM(i+1)-t_KM(i); % zjištění časového kroku
end
dt_KM(i+1)=dt_KM(i);

for i=1:length(t_KM)
p_prubeh_KM(i)=p_inf(t_KM(i)); % průběh budícího tlaku
end

for i=1:length(R_KM)
p_i_KM(i)=p_g0*(R_0/R_KM(i))^(3*kappa); %tlak uvnitř bubliny
end

for i=1:length(R_KM)
T_i_KM(i)=(R_0/R_KM(i))^(3*kappa-3)*T_0; %teplota uvnitř bubliny [K]
end
fprintf('KM \t %0.3e \t %0.3e \t %0.3e \t %0.3e \t %0.3e\n',max(R_KM)/R_0,min(R_KM)/R_0,max(abs(dR_KM))/R_0,max(T_i_KM),max(p_i_KM))
end
%% VÝPOČTY TS %%
if xTS==1
for i=1:(length(t_TS)-1)
    dt_TS(i)=t_TS(i+1)-t_TS(i); % zjištění časového kroku
end
dt_TS(i+1)=dt_TS(i);

for i=1:length(t_TS)
p_prubeh_TS(i)=p_inf(t_TS(i)); % průběh budícího tlaku
end

for i=1:length(R_TS)
p_i_TS(i)=p_g0*(R_0/R_TS(i))^(3*kappa); %tlak uvnitř bubliny
end

for i=1:length(R_TS)
T_i_TS(i)=(R_0/R_TS(i))^(3*kappa-3)*T_0; %teplota uvnitř bubliny [K]
end
fprintf('TS \t %0.3e \t %0.3e \t %0.3e \t %0.3e \t %0.3e\n',max(R_TS)/R_0,min(R_TS)/R_0,max(abs(dR_TS))/R_0,max(T_i_TS),max(p_i_TS))
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
plot([tau tau], [0 max(R_RP)/R_0],'Color','k');
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 xHR==1 && xO==1
figure(2)
subplot(3,1,1)
plot(t_HR,R_HR/R_0,'Color',co{1});
hold on
plot([tau tau], [0 max(R_HR/R_0)],'Color','k');
title({'HR',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_HR,p_prubeh_HR,'Color',co{1})
xlabel('t [s]');ylabel('p_i_n_f [Pa]')
grid on

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

if xGM==1 && xO==1
figure(3)
subplot(3,1,1)
plot(t_GM,R_GM/R_0,'Color',co{1});
hold on
plot([tau tau], [0 max(R_GM/R_0)],'Color','k');
title({'GM',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_GM,p_prubeh_GM,'Color',co{1})
xlabel('t [s]');ylabel('p_i_n_f [Pa]')
grid on

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

if xKM==1 && xO==1
figure(4)
subplot(3,1,1)
plot(t_KM,R_KM/R_0,'Color',co{1});
hold on
plot([tau tau], [0 max(R_KM/R_0)],'Color','k');
title({'KM',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_KM,p_prubeh_KM,'Color',co{1})
xlabel('t [s]');ylabel('p_i_n_f [Pa]')
grid on

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

if xTS==1 && xO==1
figure(5)
subplot(3,1,1)
plot(t_TS,R_TS/R_0,'Color',co{1});
hold on
plot([tau tau], [0 max(R_TS/R_0)],'Color','k');
title({'TS',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_TS,p_prubeh_TS,'Color',co{1})
xlabel('t [s]');ylabel('p_i_n_f [Pa]')
grid on

subplot(3,1,3)
plot(t_TS,abs(dR_TS),'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(6)
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(7)
hold on
pli=1;
PL=struct;
PL.pli=plot([tau tau], [0 5],'Color','k');
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 xHR==1
    pli=pli+1;
    PL.pli=plot(t_HR,R_HR/R_0,'-','LineWidth',lnwdth,'Color',co{pli-1});
    L{pli}='Herring';
end
if xGM==1
    pli=pli+1;
    PL.pli=plot(t_GM,R_GM/R_0,'LineWidth',lnwdth,'Color',co{pli-1});
    L{pli}='Gilmore';
end
if xKM==1
    pli=pli+1;
    PL.pli=plot(t_KM,R_KM/R_0,'LineWidth',lnwdth,'Color',co{pli-1});
    L{pli}='Keller-Miksis';
end
if xTS==1
    pli=pli+1;
    PL.pli=plot(t_TS,R_TS/R_0,'LineWidth',lnwdth,'Color',co{pli-1});
    L{pli}='Tomita-Shima';
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
legend(L)
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
%%HERRING%%
function dRdt=rhs_HR(t,R)
P=p_v+p_g0*(R_0/R(1))^(3*kappa)-2*sigma/R(1)-4*mu*R(2)/R(1);
dPx=-3*kappa*p_g0*(R_0/R(1))^(2*kappa)*R_0*R(2)/R(1)^2+2*sigma*R(2)/R(1)^2+4*mu*R(2)^2/R(1)^2;
dRdt_1 = R(2);
dRdt_2=-((3*R(2)^2*((4*R(2))/(3*c_inf) - 1))/2 - (p_inf(t) - P + (R(1)*dPx*(R(2)/c_inf - 1))/c_inf)/rho)/(R(1)*((2*R(2))/c_inf - 1) + (4*mu*(R(2)/c_inf - 1))/(c_inf*rho));
dRdt=[dRdt_1; dRdt_2];
end
%%GILMORE%%
function dRdt=rhs_GM(t,R)    
P=p_v+p_g0*(R_0/R(1))^(3*kappa)-2*sigma/R(1)-4*mu*R(2)/R(1);
dPx=-3*kappa*p_g0*(R_0/R(1))^(2*kappa)*R_0*R(2)/R(1)^2+2*sigma*R(2)/R(1)^2+4*mu*R(2)^2/R(1)^2;
C=c_inf*((P+B)/(p_inf(t)+B))^((n-1)/(2*n));
H=1/rho*n/(n-1)*(p_inf(t)+B)*(((P+B)/(p_inf(t)+B))^((n-1)/n)-1);

dRdt_1 = R(2);
dRdt_2=-(H*(R(2)/C + 1) + (3*R(2)^2*(R(2)/(3*C) - 1))/2 + (R(1)*(R(2)/C - 1)*((dp_inf(t)*n)/(rho*(n - 1)) - (n*(dp_inf(t)*((B + P)/(B + p_inf(t)))^((n - 1)/n) - (p_inf(t)*(n - 1)*(dp_inf(t)*(B + P) - dPx*(B + p_inf(t))))/(n*(B + p_inf(t))^2*((B + P)/(B + p_inf(t)))^(1/n))))/(rho*(n - 1)) + (B*(dp_inf(t)*(B + P) - dPx*(B + p_inf(t))))/(rho*(B + p_inf(t))^2*((B + P)/(B + p_inf(t)))^(1/n))))/C)/(R(1)*(R(2)/C - 1) + (R(1)*((4*B*mu)/(R(1)*rho*(B + p_inf(t))*((B + P)/(B + p_inf(t)))^(1/n)) + (4*mu*p_inf(t))/(R(1)*rho*(B + p_inf(t))*((B + P)/(B + p_inf(t)))^(1/n)))*(R(2)/C - 1))/C);
dRdt=[dRdt_1; dRdt_2];
end
%%KELLER-MIKSIS%%
function dRdt=rhs_KM(t,R)    
P=p_v+p_g0*(R_0/R(1))^(3*kappa)-2*sigma/R(1)-4*mu*R(2)/R(1);
dPx=-3*kappa*p_g0*(R_0/R(1))^(2*kappa)*R_0*R(2)/R(1)^2+2*sigma*R(2)/R(1)^2+4*mu*R(2)^2/R(1)^2;
C=c_inf*((P+B)/(p_inf(t)+B))^((n-1)/(2*n));

dRdt_1 = R(2);
dRdt_2=-(R(2)^2*(R(2)/(2*C) - 3/2) + ((P - p_inf(t))*(R(2)/C + 1))/rho + (R(1)*dPx)/(C*rho))/(R(1)*(R(2)/C - 1) - (4*mu)/(C*rho));
dRdt=[dRdt_1; dRdt_2];
end
%%TOMITA-SHIMA%%
function dRdt=rhs_TS(t,R)    
P=p_v+p_g0*(R_0/R(1))^(3*kappa)-2*sigma/R(1)-4*mu*R(2)/R(1);
dPx=-3*kappa*p_g0*(R_0/R(1))^(2*kappa)*R_0*R(2)/R(1)^2+2*sigma*R(2)/R(1)^2+4*mu*R(2)^2/R(1)^2;
C=c_inf*((P+B)/(p_inf(t)+B))^((n-1)/(2*n));

dRdt_1 = R(2);
dRdt_2=-(R(2)^2*((21*R(2)^2)/(10*C^2) - (2*R(2))/C + 3/2) - ((P - p_inf(t))*(R(2)^2/(2*C^2) - ((3*P)/2 - (3*p_inf(t))/2)/(C^2*rho) + 1))/rho + (R(1)*dPx*((2*R(2))/C - 1))/(C*rho))/(R(1)*((23*R(2)^2)/(10*C^2) - (2*R(2))/C + 1) - (4*mu*((2*R(2))/C - 1))/(C*rho));
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