clear;clc;close all;
% 设置各区域的同位素年龄数组
tdata_raw        = [4.35; 4.10; 3.90; 3.91; 3.91; 3.70; 3.72; 3.53; 3.28; 3.18; 3.40; 3.30; 0.85; 0.109; 0.05; 0.026; 0.375;2.03];  % 原始同位素年龄
tdata_Hie        = [      4.10; 3.90; 3.91; 3.91; 3.70; 3.72; 3.53; 3.28; 3.18; 3.40; 3.30; 0.85; 0.109; 0.05; 0.026; 0.375;2.03]; % Hiesinger et al. (2012)
% tdata_Hie        = [ 3.70; 3.72; 3.53; 3.28; 3.18; 3.40; 3.30; 0.85; 0.109; 0.05; 0.026; 0.375;2.03]; % Hiesinger et al. (2012)


cdata_raw        = [0.36; 0.12; 0.034; 0.035; 0.037; 0.01; 0.009; 0.0064; 0.0032; 0.0036; 0.0033; 0.003; 0.000668;  0.00009;   0.000044;  0.000021; 0.00036; 0.00174]; % 原始N(1)值  Jia:0.00174;   Wu:0.00124;   Qian:0.00128
cdata_Hie        = [      0.12; 0.034; 0.035; 0.037; 0.01; 0.009; 0.0064; 0.0032; 0.0036; 0.0033; 0.003; 0.000668;  0.00009;   0.000044;  0.000021; 0.00036; 0.00174]; % Hiesinger et al. (2012)
% cdata_Hie        = [0.01; 0.009; 0.0064; 0.0032; 0.0036; 0.0033; 0.003; 0.000668;  0.00009;   0.000044;  0.000021; 0.00036; 0.00174]; % Hiesinger et al. (2012)

% -------------------------------------------------------------------------
% 用Neukum的值作为初值，用最小二乘法求解指数+线性函数；x0表示最小二乘法的初值，x表示得到的系数，f表示利用最小二乘法反求得到的N(1)值
x0 = [5.44*10^(-14); 6.93; 8.38*10^(-4)];
% x0 = [5.355688889153016e-14; 6.929988330588613e+00; 8.900405084373684e-04];
% 求出未知数的值
options = optimoptions('lsqcurvefit','MaxIter',1400,'TolFun', 1e-36, 'TolPCG',0.1,'TolX',1e-36,'Diagnostics','on',...
    'Display','iter-detailed','FunValCheck','on','PlotFcn','optimplotresnorm','MaxFunEvals',5000);
% options = optimoptions('lsqcurvefit','Display','iter-detailed','PlotFcn','optimplotresnorm','TolFun', 1e-36,'TolX',1e-36,'MaxFunEvals',5000);
[x,resnorm,residual,exitflag,output,lambda,jacobian] = lsqcurvefit('curvefun1',x0,tdata_Hie,cdata_Hie,[],[],options);       % x1 表示添加新数据后的拟合

% -------------------------------------------------------------------------
% 将拟合函数结果和原始数据画出来
plotdata = min(tdata_raw):0.01:max(tdata_raw);  % 将数据加密
plotNew  = curvefun1(x,plotdata);          % 将拟合函数画出来
plotNeu  = curvefun1(x0,plotdata);         % 将Neukum函数画出来

% 这里要画出误差
variance_Age = zeros(1,size(tdata_raw,1));
variance_Age = variance_Age + 0.1;
variance_Age(8) = 0.05;
variance_Age(11) = 0.04;
variance_Age(13) = 0.2;
variance_Age(14) = 0.004;
variance_Age(15) = 0.0014;
variance_Age(16) = 0.0008;
variance_Age(17) = 0.075;
variance_Age(18) = 0.004;

variance_N1 = zeros(1,size(tdata_raw,1));
variance_N1(1) = 1.1*10^(-1); variance_N1(2) = 0.4*10^(-1); variance_N1(3) = 0.7*10^(-2);
variance_N1(4) = 0.5*10^(-2); variance_N1(5) = 0.7*10^(-2); variance_N1(6) = 0.3*10^(-2);
variance_N1(7) = 1.8*10^(-3); variance_N1(8) = 2.0*10^(-3); variance_N1(9) = 1.1*10^(-3);
variance_N1(10) = 1.1*10^(-3); variance_N1(11) = 1.0*10^(-3); variance_N1(12) = 0.6*10^(-3);
variance_N1(13) = 0.3*10^(-3); variance_N1(14) = 1.8*10^(-5); variance_N1(15) = 1.1*10^(-5);
variance_N1(16) = 0.5*10^(-5); variance_N1(17) = 1.1*10^(-4); variance_N1(18) = 3.2*10^(-4);






figure (2)
h2=semilogy(plotdata,plotNew,'r','lineWidth',2);  % 画出新函数拟合

set(gca,'XDir','reverse')              % 对X方向反转
xlabel('Age (Ga)','fontsize',16)
ylabel('N(1) (km^-^2)','fontsize',16)
g2=get(h2,'Parent');
set(g2,'LineWidth',1.5)
hold on

semilogy(plotdata,plotNeu,'b-.','lineWidth',2)       % 画出Neukum公式
% semilogy(tdata_raw,cdata_raw,'k+','lineWidth',1.5)                   % 画出原始数据
errorbar(tdata_raw,cdata_raw,variance_N1,variance_N1,variance_Age,variance_Age,'ko','lineWidth',1.5);

handle = legend('This study','Neukum (1983)','Radiometric ages');
set(handle,'FontSize',12)
hold off;


% 将新曲线与Neukum曲线之间的差画出来：目的是得到在相同的N(1)情况下，年龄结果的差异
% 横轴是N(1)，需要离散化
N1 = -5:0.01:0;
N1 = 10.^N1;
% 然后分别由新函数和Neukum函数反算年龄
num = size(N1,2);
plot_dif = zeros(1,num);
plot_Neu = zeros(1,num);


for i=1:num
    % 先用Neukum公式求得N(1)对应的年龄
    N_age = Bisection(x0,N1(i));
    % 再用新的公式求得N(1)对应的年龄
    F_age = Bisection(x,N1(i));
    % 然后将二者做差得到相同N1下的年龄差值
    plot_dif(i) = F_age - N_age;
    plot_Neu(i) = N_age;
end

% 将结果画出来
figure(3)
h3=semilogx(N1,plot_dif,'r','lineWidth',2);  % 画出线性函数拟合
set(gca,'XDir','reverse')              % 对X方向反转
xlabel('N(1) (km^-^2)','fontsize',16)
ylabel('Age difference (Ga)','fontsize',16)
g3=get(h3,'Parent');
set(g3,'LineWidth',1.5)

xlim=get(gca,'Xlim');
hold on
plot(xlim,[0,0],'k-.','LineWidth',2)

figure(4)
% 以Neukum的年龄为横坐标，以Jia模型年龄与Neukum模型年龄之差作为纵坐标画出来

h4=plot(plot_Neu,plot_dif,'r','lineWidth',2);  % 画出线性函数拟合
set(gca,'XDir','reverse')              % 对X方向反转
xlabel('Age (Ga)','fontsize',16)
ylabel('Age difference (Ga)','fontsize',16)
g4=get(h4,'Parent');
set(g4,'LineWidth',1.5)

% 画出在1 km处的撞击频率随时间变化
F_rate = curvefun3(x,plotdata);         % 得到新函数的撞击率；x-新函数拟合系数；plotdata：年龄步长；F_rate：新函数撞击率

H_rate = curvefun3(x0,plotdata);         % 不用新数据点得到的函数的撞击率；x2-Neukum函数系数；plotdata：年龄步长；H_rate：Hiesinger函数撞击率

xlim=get(gca,'Xlim');
hold on
plot(xlim,[0,0],'k-.','LineWidth',2)


figure(5)
h5=semilogy(plotdata,F_rate,'r','lineWidth',2);  % 画出新函数的撞击率
hold on
semilogy(plotdata,H_rate,'b-.','lineWidth',2)  % 画出新函数的撞击率
set(gca,'XDir','reverse')              % 对X方向反转
xlabel('Age (Ga)','fontsize',16)
ylabel('Cumulative cratering rate (km^-^2 Ga^-^1)','fontsize',16)
handle = legend('This study','Neukum (1983)');
set(handle,'FontSize',12)
g5=get(h5,'Parent');
set(g5,'LineWidth',1.5)
hold off

figure(6)

h6=plot(plotdata,F_rate - H_rate,'r','lineWidth',2);  % 画出Neukum公式
   
set(gca,'XDir','reverse')              % 对X方向反转
xlabel('Age (Ga)','fontsize',16)
ylabel('Cumulative cratering rate difference (km^-^2 Ga^-^1)','fontsize',12)
g6=get(h6,'Parent');
set(g6,'LineWidth',1.5)

xlim=get(gca,'Xlim');
hold on
plot(xlim,[0,0],'k-.','LineWidth',2)


hold off

