jon_justification.m

clear;close all;clc;
addpath('../IMRv2/src/forward_solver')

expDataPrep
radial = 2;
stress = 1; % 1: models 1,2,3; 2: models 4,6; 3: models 5,7
bubtherm = 1;
medtherm = 0;
masstrans = 1;
vapor = 1;
Nt = 100;
Mt = 100;

tvector = linspace(0,tvector(end),500);

[t, R, U, p,~,~,~,S] = f_imr_fd('radial', radial, 'stress', stress, ...
        'bubtherm', bubtherm, 'medtherm', medtherm, ...
        'masstrans', masstrans, 'vapor', vapor, 'tvector', tvector, ...
        'r0', R0, 'req', Req, 'mu', 1e-1, 'g', 1e4, 'lambda1', 0, ...
        'lambda2', 0, 'alphax', 0, 'collapse', 0, 'nt', Nt, 'mt', Mt);

strainRate = -2*U./R;
strain = log(R./(Req/R0));
stress = S;

%% 
figure;

% Left y-axis
yyaxis left
plot(S, -2*U./R, 'b', 'LineWidth', 1.5);
ylabel('strain rate');
xlabel('stress');
grid on;

% Right y-axis
yyaxis right
plot(S, log(R./(Req/R0)), 'r', 'LineWidth', 1.5);
ylabel('strain');

% Title and legend
legend({'strain rate', 'strain'}, 'Location', 'best');

%%
figure; hold on;
plot(t,stress./max(abs(stress)));
plot(t,strain./max(abs(strain)));
plot(t,strainRate./max(abs(strainRate)));
legend('stress','strain','strainRate')

%%
figure; hold on;
% plot(t,stress);
% plot(t,strain);
plot(strain,strainRate,'k','LineWidth',1)
% plot(t,strainRate);
% legend('stress','strain','strainRate')
%%
close all;
G = 1e4;
mu = 1e-1;

S_G = (-G/2).*(5-((Req/R0)./R).^4-4*((Req/R0)./R));
S_mu = -4*mu.*U./R;

threshold = 1e5;
Stotal = S_G + S_mu;
filter = nan(size(Stotal));
for i = 1:length(Stotal)
    if (abs(strainRate(i))/mean(tc)>threshold) 
        filter(i) = 1;
    end
    if (abs(strain(i)) > 0.1*2.5)
        filter(i) = 1;
    end
end
Stotal_filter = Stotal.*filter;
figure; hold on;
plot(t,Stotal,'k')
plot(t,Stotal_filter,'rs')
legend('total','filter')