Commit 48ca32e2 authored by Lena Privat Laptop's avatar Lena Privat Laptop
Browse files

Add evaluation routines for exoplanet core paper (Noack and Lasbleis)

parent 0679eee8
% need to start statistics package first:
% pkg load optim
function core_regression
% 1 2 3 4 5 6 7 8 9 10 11 12 13 14
% M C_FePl(wt%) C_FeM (mole%) Psurf Rp Rc Rrc gs rho_c rho_m Cp_c Cp_m alpha_c alpha_m
% 15 16 17 18 19 20 21 22 23 24 25 26 27 28
% k_m shear_m gamma_m Ts Tcmb p_c DTc T_center p_center gamma_c r_ICB p_ICB RR_ICB nr???
% [fname, fpath, fltidx] = uigetfile("data_IS.res");
data1 = load('-ascii',"C:/OwnCloud_FUB/Literature_Projects/Eigene_Projekte/Project_Core_Exoplanets/ProfilesTemp_Jan2019/Ini_With_DTcmb/data_IS.res");
data2 = load('-ascii',"C:/OwnCloud_FUB/Literature_Projects/Eigene_Projekte/Project_Core_Exoplanets/ProfilesTemp_Jan2019/Ini_With_DTcmb_moreFeCases/data_IS.res");
data = [data1;data2];
%%% CMF
% get mass fraction of iron in entire planet -> need to subtract iron in mantle
X_Fe = data(:,2)/100; % between 0.15 and 0.8
M_Fe = 55.845; %g/mol
M_Mg = 24.305;
M_Si = 28.0855;
M_O = 15.999;
NrFeM = data(:,3); % iron number in mantle, ratio based on M2SiO4+Fe2SiO4 minerals
%mass fraction iron in mantle
X_FeM = 2*NrFeM/100*M_Fe ./ ( 2*( (1-NrFeM/100)*M_Mg + NrFeM/100*M_Fe ) + M_Si + 4*M_O ); % between 0 and 0.1457
% X_CMF*M_p = C_Fe*M_p-X_FeM*M_m; M_m=(1-X_CMF)*M_p, hence:
X_CMF = (X_Fe-X_FeM)./(1-X_FeM); % between 0.005 and 0.8
%%% Planet radius
%% 1. version -> which exponent for Mass?
% start with loop over Mp
Mp = data(:,1);
Mass = [0.8,0.9,1.0,1.1,1.2,1.3,1.4,1.5,1.6,1.8,2.0];
b_1 = zeros(size(Mass));
b_2 = zeros(size(Mass));
y_all = data(:,5)/1000; % Rp in km
for i=1:11
x = X_Fe(Mp==Mass(i));, y = y_all(Mp==Mass(i));
modelfun = @(b, x) ((b(1) - b(2)*x*100));
beta0 = [7000;18];
[beta, R, J, CovB, MSE] = nlinfit(x, y, modelfun, beta0);
% fprintf('Rp: %d %f %f \n',Mass(i),beta)
b_1(i)=beta(1); b_2(i)=beta(2);
end
modelfun = @(b, x) (b(1)*x.^b(2));
beta0 = [7000;0.3];
[beta] = nlinfit(Mass, b_1, modelfun, beta0);
% fprintf('b1: %f %f \n',beta)
modelfun = @(b, x) (b(1)*x.^b(2));
beta0 = [18;0.3];
[beta] = nlinfit(Mass, b_2, modelfun, beta0);
% fprintf('b2: %f %f \n',beta)
% Rp = 7030*Mp.^0.28 - 18.4*X_Fe*100.*Mp.^0.295;
% dy = abs(Rp-y_all)./y_all; % error w.r.t. original data y
% fprintf('Rp: av/max error: %2.1f%% %2.1f%%\n',mean(dy)*100,max(dy)*100)
% Rp = Rp*1000;
%% 2. version -> fixed exponent bexp
bexp = 0.28;
y_all = data(:,5)/1000; % Rp in km
Rp_Mp = y_all ./ data(:,1).^bexp; % divide by Mp^bexp
modelfun = @(b, x) ((b(1) - b(2)*x*100));
beta0 = [6000;20];
[beta] = nlinfit(X_Fe, Rp_Mp, modelfun, beta0);
fprintf('Rp: %f %f \n',beta)
beta0 = [7030;18.4]; % rounded values
Rp = modelfun(beta0, X_Fe);
Rp = Rp.*Mp.^bexp;
dy = abs(Rp-y_all)./y_all; % error w.r.t. original data y
fprintf('Rp: av/max error: %2.1f%% %2.1f%% for %f/%f/%f\n',mean(dy)*100,max(dy)*100,beta0,bexp)
Rp = Rp*1000;
%%% rho_c_av
modelfun = @(b, x) (b(1) * x.^b(2));
y = data(:,9);
beta0 = [12300;0.2];
[beta, R, J, CovB, MSE] = nlinfit(Mp, y, modelfun, beta0);
fprintf('rho_c,av: %f %f\n',beta)
rho_c_av = modelfun(beta0, Mp);
dy = abs(y-rho_c_av)./y; % error w.r.t. original data y
fprintf('rho_c,av: av/max error: %2.1f%% %2.1f%% for %f/%f\n',mean(dy)*100,max(dy)*100,beta0)
% figure
% subplot(1,4,1), plot(Mp(Mp==1),y(Mp==1),'*',Mp(Mp==1),rho_c_av(Mp==1),'ro')
% subplot(1,4,2), plot(1-X_CMF(Mp==1),y(Mp==1),'*')
% subplot(1,4,3), x=data(:,5);, plot(x(Mp==1),y(Mp==1),'*')
% subplot(1,4,4), x=data(:,19);, plot(x(Mp==1),y(Mp==1),'*')
% figure
% for k=1:5
% for j=1:5
% i = (k-1)*5+j;
% subplot(5,5,i)
% x=data(:,i+1);
% plot(x(Mp==1),y(Mp==1),'*')
% end
% end
% new version, include Rp and Mp dependence:
% start with loop over Mp
Mass = [0.8,0.9,1.0,1.1,1.2,1.3,1.4,1.5,1.6,1.8,2.0];
b_1 = zeros(size(Mass));
b_2 = zeros(size(Mass));
y_all = data(:,9); % rho_c
Rp_data = Rp/6371000; % planet radius in Earth radii
for i=1:11
x = Rp_data(Mp==Mass(i));, y = y_all(Mp==Mass(i));
modelfun = @(b, x) ((b(1) - b(2)*x));
beta0 = [15000;2000];
beta = nlinfit(x, y, modelfun, beta0);
% fprintf('rho_c(M): %d %f %f \n',Mass(i),beta)
b_1(i)=beta(1); b_2(i)=beta(2);
end
modelfun = @(b, x) (b(1)*x.^b(2));, beta0 = [1;1];, [beta] = nlinfit(Mass, b_1, modelfun, beta0);
fprintf('b1: %f %f \n',beta)
modelfun = @(b, x) (b(1)*x.^b(2));, beta0 = [1;1];, [beta] = nlinfit(Mass, b_2, modelfun, beta0);
fprintf('b2: %f %f \n',beta)
rho_c_av = 14400*Mp.^0.18 - 2170*Rp_data.*Mp.^-0.2;
dy = abs(rho_c_av-y_all)./y_all; % error w.r.t. original data y
fprintf('rho_c_av: av/max error: %2.1f%% %2.1f%%\n',mean(dy)*100,max(dy)*100)
%% 3rd version
rho_div_M = data(:,9)./Mp.^0.2;
modelfun = @(b, x) (b(1)+x.*b(2));, beta0 = [1;1];, [beta] = nlinfit(Rp_data, rho_div_M, modelfun, beta0);
fprintf('b1: %f %f \n',beta)
rho_c_av = (13100 - 780*Rp_data).*Mp.^0.2;
dy = abs(rho_c_av-y_all)./y_all; % error w.r.t. original data y
fprintf('rho_c_av: av/max error: %2.1f%% %2.1f%%\n',mean(dy)*100,max(dy)*100)
return
%%% Core radius
Rc = (X_CMF.*Mp*5.982*10^24./(4/3*pi*rho_c_av)).^(1/3);
y = data(:,6);
dy = abs(y-Rc)./y;
fprintf('Rc: av/max error: %2.1f%% %2.1f%%\n',mean(dy)*100,max(dy)*100)
% average mantle density
rho_m_av = (1.0-X_CMF).*Mp*5.982*10^24 ./ (4/3 * pi * (Rp.^3-Rc.^3));
y = data(:,10);
dy = abs(rho_m_av-y)./y; % error w.r.t. original data y
fprintf('rho_m,av: av/max error: %2.1f%% %2.1f%%\n',mean(dy)*100,max(dy)*100)
% gravitational acceleration
Gconst = 6.67384*10^-11;
g0 = Gconst*Mp*5.982*10^24./Rp.^2;
dy = abs(g0-data(:,8))./data(:,8); % error w.r.t. original data y
fprintf('g0: av/max error: %2.1f%% %2.1f%%\n',mean(dy)*100,max(dy)*100)
gc = Gconst*X_CMF.*Mp*5.982*10^24./Rc.^2;
gc_data = Gconst*X_CMF.*Mp*5.982*10^24./data(:,6).^2;
dy = abs(gc-gc_data)./gc_data; % error w.r.t. original data y
fprintf('gc: av/max error: %2.1f%% %2.1f%%\n',mean(dy)*100,max(dy)*100)
gm = 0.5*(g0+gc);
% CMB pressure
P_CMB = gm.*rho_m_av.*(Rp-Rc)*10^-9; % pressure in GPa
dy = abs(P_CMB-data(:,20))./data(:,20); % error w.r.t. original data y
fprintf('Pc: av/max error: %2.1f%% %2.1f%%\n',mean(dy)*100,max(dy)*100)
% temperature at CMB
T_CMB = 5400*(P_CMB/140).^0.48 ./ (1-log(1-NrFeM/100)); % log = logarithm naturalis (ln)
dy = abs(T_CMB-data(:,19))./data(:,19); % error w.r.t. original data y
fprintf('Tc: av/max error: %2.1f%% %2.1f%%\n',mean(dy)*100,max(dy)*100)
%%% Thermal expansion coefficient
% start with loop over Mp
Mass = [0.8,0.9,1.0,1.1,1.2,1.3,1.4,1.5,1.6,1.8,2.0];
b_1 = zeros(size(Mass));
b_2 = zeros(size(Mass));
y_all = data(:,13); % alpha_c
for i=1:11
x = rho_c_av(Mp==Mass(i));, y = y_all(Mp==Mass(i));
modelfun = @(b, x) ((b(1) - b(2)*x/1000));
beta0 = [10;0.01];
beta = nlinfit(x, y, modelfun, beta0);
% fprintf('alpha: %d %f %f \n',Mass(i),beta)
b_1(i)=beta(1); b_2(i)=beta(2);
end
modelfun = @(b, x) (b(1)*x.^b(2));
beta0 = [1;-0.5];
[beta] = nlinfit(Mass, b_1, modelfun, beta0);
% fprintf('b1: %f %f \n',beta)
modelfun = @(b, x) (b(1)*x.^b(2));
beta0 = [1;-1];
[beta] = nlinfit(Mass, b_2, modelfun, beta0);
% fprintf('b2: %f %f \n',beta)
alpha = 5.2*Mp.^0.15 - 0.28/1000*rho_c_av.*Mp.^0.19;
dy = abs(alpha-y_all)./y_all; % error w.r.t. original data y
% fprintf('alpha: av/max error: %2.1f%% %2.1f%%\n',mean(dy)*100,max(dy)*100)
alpha = 8.48*Mp.^-0.69 - 0.000555*rho_c_av.*Mp.^-0.95;
dy = abs(alpha-y_all)./y_all; % error w.r.t. original data y
% fprintf('alpha: av/max error: %2.1f%% %2.1f%%\n',mean(dy)*100,max(dy)*100)
x = rho_c_av;, y = y_all;
modelfun = @(b, x) ((b(1) - b(2)*x));
beta0 = [10;0.01];
beta = nlinfit(x, y, modelfun, beta0);
% fprintf('alpha: %f %f \n',beta)
beta0 = [5;0.00027];
alpha_c = modelfun(beta0, rho_c_av);
dy = abs(alpha_c-y_all)./y_all; % error w.r.t. original data y
% fprintf('alpha_c: av/max error: %2.1f%% %2.1f%%\n',mean(dy)*100,max(dy)*100)
% figure
% plot(rho_c_av,data(:,13),'*')
end
\ No newline at end of file
% need to start statistics package first:
% pkg load optim
function core_regression_v2
% 1 2 3 4 5 6 7 8 9 10 11 12 13 14
% M C_FePl(wt%) C_FeM (mole%) Psurf Rp Rc Rrc gs rho_c rho_m Cp_c Cp_m alpha_c alpha_m
% 15 16 17 18 19 20 21 22 23 24 25 26 27 28
% k_m shear_m gamma_m Ts Tcmb p_c DTc T_center p_center gamma_c r_ICB p_ICB RR_ICB nr???
% [fname, fpath, fltidx] = uigetfile("data_IS.res");
data1 = load('-ascii',"C:/OwnCloud_FUB/Literature_Projects/Eigene_Projekte/Project_Core_Exoplanets/ProfilesTemp_Jan2019/Ini_With_DTcmb/data_IS.res");
data2 = load('-ascii',"C:/OwnCloud_FUB/Literature_Projects/Eigene_Projekte/Project_Core_Exoplanets/ProfilesTemp_Jan2019/Ini_With_DTcmb_moreFeCases/data_IS.res");
data = [data1;data2];
%%% CMF
% get mass fraction of iron in entire planet -> need to subtract iron in mantle
X_Fe = data(:,2)/100; % between 0.15 and 0.8
M_Fe = 55.845; %g/mol
M_Mg = 24.305;
M_Si = 28.0855;
M_O = 15.999;
NrFeM = data(:,3); % iron number in mantle, ratio based on M2SiO4+Fe2SiO4 minerals
%mass fraction iron in mantle
X_FeM = 2*NrFeM/100*M_Fe ./ ( 2*( (1-NrFeM/100)*M_Mg + NrFeM/100*M_Fe ) + M_Si + 4*M_O ); % between 0 and 0.1457
% X_CMF*M_p = C_Fe*M_p-X_FeM*M_m; M_m=(1-X_CMF)*M_p, hence:
X_CMF = (X_Fe-X_FeM)./(1-X_FeM); % between 0.005 and 0.8
Mp = data(:,1);
% for loops over mass
Mass = [0.8,0.9,1.0,1.1,1.2,1.3,1.4,1.5,1.6,1.8,2.0];
%%% Planet radius
%Rp = a * M^c
% a = a_1 + X_Fe*a_2
% c = c_1 + X_Fe*c_2
% For each mass, get Rp = a_1(M) + X_Fe*a_2(M)
Fe_arr = 0.15:0.05:0.8;
b_1 = zeros(size(Fe_arr));
b_2 = zeros(size(Fe_arr));
x_all = Mp;
y_all = data(:,5)/1000; %Rp in km
tol = 0.01;
for i=1:length(Fe_arr)
x = x_all(X_Fe>Fe_arr(i)-tol & X_Fe<Fe_arr(i)+tol);
y = y_all(X_Fe>Fe_arr(i)-tol & X_Fe<Fe_arr(i)+tol);
modelfun = @(b, x) ((b(1)*x.^b(2)));
beta = nlinfit(x, y, modelfun, [1;1]);
% fprintf('Rp: %d %f %f \n',Fe_arr(i),beta)
b_1(i)=beta(1); b_2(i)=beta(2);
end
% prefactor a:
modelfun = @(b, x) ((b(1)+x.*b(2)));
a_Rp = nlinfit(Fe_arr, b_1, modelfun, [1;1]);
c_Rp = nlinfit(Fe_arr, b_2, modelfun, [1;1]);
Rp = (a_Rp(1)+a_Rp(2).*X_Fe).*Mp.^(c_Rp(1)+c_Rp(2).*X_Fe);
dy = abs(Rp-y_all)./y_all; % error w.r.t. original data y
fprintf('Rp: av/max error: %2.1f%% %2.1f%% for %f/%f/%f/%f\n',mean(dy)*100,max(dy)*100,a_Rp,c_Rp)
Rp = (a_Rp(1)+a_Rp(2).*X_Fe).*Mp.^(0.28);
dy = abs(Rp-y_all)./y_all; % error w.r.t. original data y
fprintf('Rp: av/max error: %2.1f%% %2.1f%% for %f/%f/%f\n',mean(dy)*100,max(dy)*100,a_Rp,0.28)
%%% Core radius
%Rc = a * M^c
% a = a_1 + X_CMF*a_2
% c = c_1 + X_CMF*c_2
b_1 = zeros(size(Mass));
b_2 = zeros(size(Mass));
x_all = X_CMF;
y_all = data(:,6)/1000; %Rc in km
tol = 0.01;
for i=1:length(Mass)
x = x_all(Mp>Mass(i)-tol & Mp<Mass(i)+tol);
y = y_all(Mp>Mass(i)-tol & Mp<Mass(i)+tol);
modelfun = @(b, x) ((b(1)*x.^b(2)));
beta = nlinfit(x, y, modelfun, [6000;0.3]);
Rc(Mp>Mass(i)-tol & Mp<Mass(i)+tol) = beta(1).*x'.^beta(2);
dy = abs(Rc(Mp>Mass(i)-tol & Mp<Mass(i)+tol)-y')./y'; % error w.r.t. original data y
%fprintf('Rc: %d %f %f %f -> %f %f\n',i,Mass(i),beta,mean(dy)*100,max(dy)*100)
b_1(i)=beta(1); b_2(i)=beta(2);
end
% include mass:
modelfun = @(b, x) ((b(1)*x.^b(2)));
b_M1 = nlinfit(Mass, b_1, modelfun, [1;1]);
a_Rc = b_M1(1);
b_Rc = 0.33;
c_Rc = b_M1(2);
Rc = a_Rc*X_CMF.^b_Rc.*Mp.^c_Rc;
dy = abs(Rc-y_all)./y_all; % error w.r.t. original data y
fprintf('Rc: av/max error: %2.1f%% %2.1f%% for %f/%f/%f\n',mean(dy)*100,max(dy)*100,a_Rc,b_Rc,c_Rc)
% rho_c
y_all = data(:,9);
rho_c_av = X_CMF.*Mp*5.982*10^24./(4/3*pi*(Rc.*1000).^3);
dy = abs(rho_c_av-y_all)./y_all; % error w.r.t. original data y
fprintf('rho_c,av: av/max error: %2.1f%% %2.1f%% \n',mean(dy)*100,max(dy)*100)
% average mantle density
rho_m_av = (1.0-X_CMF).*Mp*5.982*10^24 ./ (4/3 * pi * (Rp.^3-Rc.^3).*1000.^3);
y_all = data(:,10);
dy = abs(rho_m_av-y_all)./y_all; % error w.r.t. original data y
fprintf('rho_m,av: av/max error: %2.1f%% %2.1f%%\n',mean(dy)*100,max(dy)*100)
% gravitational acceleration
Gconst = 6.67384*10^-11;
g0 = Gconst*Mp*5.982*10^24./(1000*Rp).^2;
dy = abs(g0-data(:,8))./data(:,8); % error w.r.t. original data y
fprintf('g0: av/max error: %2.1f%% %2.1f%%\n',mean(dy)*100,max(dy)*100)
gc = Gconst*X_CMF.*Mp*5.982*10^24./(1000*Rc).^2;
gc_data = Gconst*X_CMF.*Mp*5.982*10^24./data(:,6).^2;
dy = abs(gc-gc_data)./gc_data; % error w.r.t. original data y
fprintf('gc: av/max error: %2.1f%% %2.1f%%\n',mean(dy)*100,max(dy)*100)
gm = 0.5*(g0+gc);
% CMB pressure
P_CMB = gm.*rho_m_av.*(Rp-Rc)*1000*10^-9; % pressure in GPa
dy = abs(P_CMB-data(:,20))./data(:,20); % error w.r.t. original data y
fprintf('Pc: av/max error: %2.1f%% %2.1f%%\n',mean(dy)*100,max(dy)*100)
% temperature at CMB
T_CMB = 5400*(P_CMB/140).^0.48 ./ (1-log(1-NrFeM/100)); % log = logarithm naturalis (ln)
dy = abs(T_CMB-data(:,19))./data(:,19); % error w.r.t. original data y
fprintf('Tc: av/max error: %2.1f%% %2.1f%%\n',mean(dy)*100,max(dy)*100)
%%% Thermal expansion coefficient
% start with loop over Mp
b_1 = zeros(size(Mass));
b_2 = zeros(size(Mass));
y_all = data(:,13); % alpha_c
for i=1:11
x = rho_c_av(Mp==Mass(i));, y = y_all(Mp==Mass(i));
modelfun = @(b, x) ((b(1) - b(2)*x/1000));
beta0 = [10;0.01];
beta = nlinfit(x, y, modelfun, beta0);
%fprintf('alpha: %d %f %f \n',Mass(i),beta)
b_1(i)=beta(1); b_2(i)=beta(2);
end
modelfun = @(b, x) (b(1)*x.^b(2));
beta0 = [1;-0.5];
[beta] = nlinfit(Mass, b_1, modelfun, beta0);
fprintf('b1: %f %f \n',beta)
modelfun = @(b, x) (b(1)*x.^b(2));
beta0 = [1;-1];
[beta] = nlinfit(Mass, b_2, modelfun, beta0);
fprintf('b2: %f %f \n',beta)
alpha = 5.2*Mp.^0.15 - 0.28/1000*rho_c_av.*Mp.^0.19;
dy = abs(alpha-y_all)./y_all; % error w.r.t. original data y
fprintf('alpha: av/max error: %2.1f%% %2.1f%%\n',mean(dy)*100,max(dy)*100)
alpha = 8.48*Mp.^-0.69 - 0.000555*rho_c_av.*Mp.^-0.95;
dy = abs(alpha-y_all)./y_all; % error w.r.t. original data y
fprintf('alpha: av/max error: %2.1f%% %2.1f%%\n',mean(dy)*100,max(dy)*100)
alpha = 9.45*Mp.^-0.93 - 0.63/1000*rho_c_av.*Mp.^-1.23;
dy = abs(alpha-y_all)./y_all; % error w.r.t. original data y
fprintf('alpha_c: av/max error: %2.1f%% %2.1f%%\n',mean(dy)*100,max(dy)*100)
x = rho_c_av;, y = y_all;
modelfun = @(b, x) ((b(1) - b(2)*x));
beta0 = [10;0.01];
beta = nlinfit(x, y, modelfun, beta0);
% fprintf('alpha: %f %f \n',beta)
beta0 = [5;0.00027];
alpha_c = modelfun(beta0, rho_c_av);
dy = abs(alpha_c-y_all)./y_all; % error w.r.t. original data y
fprintf('alpha_c: av/max error: %2.1f%% %2.1f%%\n',mean(dy)*100,max(dy)*100)
% figure
%% plot(rho_c_av,data(:,13),'*') %rho_c_av
%% plot(rho_c_av,data(:,19),'*') %T_CMB
% plot(rho_c_av,X_CMF,'*') %X_CMF
% new try via X_CMF
for i=1:11
x = X_CMF(Mp==Mass(i));, y = y_all(Mp==Mass(i));
modelfun = @(b, x) ((b(1) * exp(x).^b(2)));
beta0 = [10;0.01];
beta = nlinfit(x, y, modelfun, beta0);
error = (modelfun(beta, x)-y)./y*100;
% fprintf('alpha: %d %f %f Err: %f%% %f%% \n',Mass(i),beta,mean(error),max(error))
b_1(i)=beta(1); b_2(i)=beta(2);
end
figure
subplot(2,2,1)
plot(Mass,b_1)
subplot(2,2,2)
plot(Mass,b_2)
modelfun = @(b, x) (b(1)*x.^b(2));
beta0 = [1;-0.5];
[beta] = nlinfit(Mass, b_1, modelfun, beta0);
fprintf('b1: %f %f \n',beta)
modelfun = @(b, x) (b(1)*x.^b(2));
beta0 = [1;-1];
[beta] = nlinfit(Mass, b_2, modelfun, beta0);
fprintf('b2: %f %f \n',beta)
subplot(2,2,3)
plot(Mass,modelfun([1.78,-0.514], Mass))
subplot(2,2,4)
plot(Mass,modelfun([-0.176,-0.909], Mass))
alpha = 1.78*Mp.^-0.514 .* exp(X_CMF) .^(-0.176*Mp.^-0.909);
dy = abs(alpha-y_all)./y_all; % error w.r.t. original data y
fprintf('alpha: av/max error: %2.1f%% %2.1f%%\n',mean(dy)*100,max(dy)*100)
end
\ No newline at end of file
This diff is collapsed.
This diff is collapsed.
This source diff could not be displayed because it is too large. You can view the blob instead.