% This is the wrapper for fitting exercise 3, the simplified steady-erosion
% one. 
%
% Greg Balco - Berkeley Geochronology Center - 2018
% Not licensed for distribution

clear all; close all;

% Get data
d = create_BCO_data(20); % Don't use upper 20 g/cm2

% Atmospheric pressure
p.h = antatm(2183);

% Surface production rate not needed for this one

% P10_St = 33.2; % This is inferred from v3 calculators with CRONUS "Primary" data set
% P21_St = P10_St.*4; % Ratio is basically the average of lots of studies...

% p.P0 = P21_St;

% Other constants for muon  production

p.fC = 0.296; % Heisinger 2002b Table 2 for Si in SiO2
p.fD = 0.6559; % Same
p.Natoms = 1.003e22; % For Si in SiO2

% Other constants for nucleogenic production
p.l232 = -log(0.5)./1.405e10;
p.l235 = -log(0.5)./7.038e8;
p.l238 = -log(0.5)./4.468e9;

% Other constants
p.rho = 2.25;
p.Lsp = 140;
p.texp = 14.5e6;

% range of muon x-sections

fstars = linspace(0,0.005,10);
sigma0s = 1e-27.*linspace(0,0.015,11);

% Initial guess
ee0 = 0.1; % m/Myr
tc0 = 150; % Ma;
N21sp0 = 5; % atoms/g

% Constraints
x0 = [ee0 tc0 N21sp0];
xmin = [0 0 0];
xmax = [0.2 250 Inf];

opts = optimset('fmincon');
opts = optimset(opts,'display','none');

% Uncertainty expansions
p.dnfactor = 0.0; % Define uncert on nucleogenic Ne-21
p.dmfactor = 1; % Increase measurement uncertainty to test...

% Optimization loop

for row = 1:length(sigma0s);
    for col = 1:length(fstars);

        % We need to precompute muon stuff, otherwise it takes too long. 
        p.predz = logspace(0,5,100);
        mc.Natoms = p.Natoms;
        mc.sigma0 = sigma0s(row); 
        mc.k_neg = fstars(col).*p.fC.*p.fD;
        for a = 1:length(p.predz);
            p.Pmu(a) = P_mu_total_alpha1(p.predz(a),p.h,mc,'no');
        end;

        [optx,fval,exitflag] = fmincon(@(x) BCO_objective_steadyerosion_simple(x,d,p,0),x0,[],[],[],[],xmin,xmax,[],opts);

        eflags(row,col) = exitflag;
        fvals(row,col) = fval;
        
        opt.ee(row,col) = optx(1);
        opt.tc(row,col) = optx(2);
        opt.N21sp0(row,col) = optx(3).*1e8;
        disp([sprintf('%0.1f',fval) ' E = ' sprintf('%0.2f',optx(1)) ' m/Myr; tc = ' sprintf('%0.1f',optx(2)) ' Ma']);
        disp([int2str(row) ' ' int2str(col)]);
    end
end;

% save results

eval(['save simple_results_dnfactor_' int2str(round(p.dnfactor.*100))]);

%%

clear all; load simple_results_dnfactor_0; % Load something


dof = length(d.N21m) - 3;
fitstat = fvals./dof; % if you want to plot reduced chi-2
%fitstat = fvals; % regular chi-2

cs = [4.8 5 6 7 8 9 10 11 12]; % appropriate for zero P21 uncert
%cs = [1.5 1.6 1.7 1.8 1.9 2 4 6 8 10]; % appropriate for 20%

sigma0s_mb = sigma0s./1e-27;

figure; 
[c,h] = contour(fstars,sigma0s_mb,fitstat,cs,'linewidth',1.5);
clabel(c,h);
xlabel('f^{*}_{21}');
ylabel('\sigma_{0,21} (mb)');

% Plot closure age
hold on;

cs2 = [0:20:260];
%figure; 
[c1,h1] = contour(fstars,sigma0s_mb,opt.tc,cs2,'linestyle','--');
clabel(c1,h1,'color','red');
xlabel('f^{*}_{21}');
ylabel('\sigma_{0,21} (mb)');


% Plot various cross-section estimates
plot(0.0029,0.79./190,'ko','markersize',10,'markerfacecolor',[0.5 0.5 0.5]);
plot(0,0.0112,'ks','markersize',10,'markerfacecolor','k');
plot(0,0.0033,'ks','markersize',10,'markerfacecolor','k');

set(h,'color','k');
set(h1,'color','r');
set(gca,'fontsize',14);

% Compute and plot total production due to muons. This is not in the figure
% in the paper. 

P_mu_tot = zeros(size(opt.tc));
for row = 1:length(sigma0s)
    for col = 1:length(fstars)
        mc.Natoms = p.Natoms;
        mc.sigma0 = sigma0s(row); 
        mc.k_neg = fstars(col).*p.fC.*p.fD;
        P_mu_tot(row,col) = P_mu_total_alpha1(0,1013.25,mc);
    end
end

cs3 = [0:0.1:1];
[c3,h3] = contour(fstars,sigma0s_mb,P_mu_tot,cs3);
clabel(c3,h3,'color','green');
set(h3,'color','green');



%% Erosion rate space check

cs2 = [0 0.001 0.005 0.01 0.05 0.1 0.19 0.29 0.99];

figure; 
[c2,h2] = contour(fstars,sigma0s_mb,opt.ee,cs2,'linewidth',1.5);
clabel(c2,h2);
xlabel('f^{*}_{21}');
ylabel('\sigma_{0,21} (mb)');

hold on;
plot(0.0029,0.79./190,'ro','markerfacecolor','r');
plot(0,0.0112,'bo','markerfacecolor','b');
plot(0,0.0033,'bo','markerfacecolor','b'); % min constraint


