function out = trial_forward_integration(d,dFlag)

if nargin < 2; dFlag = 0; end

% This function does a forward integration for Be-10 and C-14
% concentrations in a core, compares them to data, and returns a misfit
% value. 
%
% Notes:
% Things that have to be passed to this
% 
% input arg d is structure with:
% d.data10 - data structure for Be-10 measurements in core
% d.data14 - data structure for C-14 measurements in core
% d.site_elv_now - meters present bedrock surface (core top) elevation
% d.t1 - output of generate_te_history
% d.l14 - C-14 decay const
% d.l10 - Be-10 decay const
% d.Dmass, d.Lmass - parameters for firn density approximation
% d.rhoice - density of glacier ice 
% d.P - precalculated production rate variable structure
% d.N10inh - inherited Be-10 concentration; normally zero but available for
%   experimenting with 
%
% Optional arg: d.plot_z_cm, which triggers return of predicted depth
% profile at z-values given
%
% Returns either a misfit statistic (if dFlag = 0, for use as an objective
% function) or a misfit statistic plus a bunch of diagnostic and plotting
% data (if dFlag = 1). 
%
% Greg Balco
% Berkeley Geochronology Center
% February, 2022
% 
% Accompanies a paper entitled 'REVERSIBLE ICE SHEET THINNING IN THE
% AMUNDSEN SEA EMBAYMENT DURING THE LATE HOLOCENE,' by Balco, Brown,
% Nichols, Venturelli, and others. 
%
% This code is provided solely for review of the accompanying paper. It
% is not licensed for other use or distribution. 


%% Now start function

% This script starts with a theoretical ice thickness change history and then
% integrates through it to compute the nuclide concentration at the end of
% the history. 

ice_thickness_m_v = d.t1.surface_elvs - d.site_elv_now;
ice_thickness_cm_v = ice_thickness_m_v .* 100;

% This is the equation for mass depth as function of ice thickness over site
mass_depth_gcm2_v = ice_thickness_cm_v.*d.rhoice - d.Dmass.*(1 - exp(-ice_thickness_cm_v./d.Lmass)); 

% Suppress negative values
iszero = find(mass_depth_gcm2_v < 0);
mass_depth_gcm2_v(iszero) = zeros(size(iszero));

% Now we have times, surface elevations, and mass thicknesses. 

% Perform forward integration. 
% The formula for computing the nuclide concentration is:
%
% Integral from 0 to t_max of P(t).*exp(-lambda.*t)
% 
% And P is actually a function of z (mass depth), so really this is:
% 
% Integral from 0 to t_max of P(z(t)).*exp(-lambda.*t)
%
% The function being integrated is provided  by PofZofT.m.
 

% Obtain sample depths from core data

% This needs to account for missing data on both sides now. 

% Get all possible z values
core_z_cm = unique([unique(d.data14.z); unique(d.data10.z)]);
core_z_gcm2 = core_z_cm.*2.77;

% Do integration at all sample depths. 

% Obtain integration bounds
maxt = max(d.t1.times);

% This integrates once for each unique sample depth
for a = 1:length(core_z_cm)
    N14(a) = integral(@(t) PofZofT(t,d.t1.times,d.t1.surface_elvs,(mass_depth_gcm2_v+core_z_gcm2(a)),d.P,14).*exp(-d.l14.*t),0,maxt);
    N10(a) = d.N10inh + integral(@(t) PofZofT(t,d.t1.times,d.t1.surface_elvs,(mass_depth_gcm2_v+core_z_gcm2(a)),d.P,10).*exp(-d.l10.*t),0,maxt);
    if isfield(d,'data26')
        N26(a) = d.N26inh + integral(@(t) PofZofT(t,d.t1.times,d.t1.surface_elvs,(mass_depth_gcm2_v+core_z_gcm2(a)),d.P,26).*exp(-d.l26.*t),0,maxt);
    end
end

% Here we need a fit statistic to see how well the model fits the data.

% Normally fit metrics look something like this:
% if you have some measurements Nm,i and some model predictions Np,i, then you
% can measure the fit between the two with something like:
%
% M = sum((Nm,i - Np,i).^2)
% 
% That expression is unweighted. Error-weighting would be something like: 
% M = sum ( ((Nm,i-Np,i)./delNm,i).^2 )

% As we have symmetrical uncertainties for the Be-10 data, we can use this
% type of a statistic for Be-10. 

% Line up Be-10 data with predictions. 
depth10_index = zeros(size(d.data10.z));
for a = 1:length(d.data10.z)
    depth10index(a) = find(d.data10.z(a) == core_z_cm);
end

dN10 = d.data10.Natomsg_bounds(:,3) - d.data10.Natomsg_bounds(:,2);
% Note: this is actually redundant with d.data10.dN, but keep it this way
% to be more general

% Calculate misfit to Be-10 data
M10s = (d.data10.Natomsg_bounds(:,3) - N10(depth10index)')./dN10;

% As we have asymmetrical uncertainties for C-14 data, we need to use some kind of an
% asymmetrical error-weighting. To do this, assume that the stated
% distribution is just asymmetric normal with SD equal to the width of the
% 65% confidence interval on the appropriate side of the data. 

% "high" applies when prediction > observation, approx uncertainty is 84-50
dN14_high = d.data14.Natomsg_bounds(:,4) - d.data14.Natomsg_bounds(:,3);
% "low" applies when prediction < observation, approx uncertainty is 50-16
dN14_low = d.data14.Natomsg_bounds(:,3) - d.data14.Natomsg_bounds(:,2);

% Line up C-14 data with predictions. 
depth14_index = zeros(size(d.data14.z));
for a = 1:length(d.data14.z)
    depth14index(a) = find(d.data14.z(a) == core_z_cm);
end

% Figure out which predictions are high/low
% "high" applies when prediction > observation
is_high = (d.data14.Natomsg_bounds(:,3) < N14(depth14index)');
is_low = ~is_high;

% Calculate misfit to C-14 data
M14s = zeros(length(depth14index),1);
M14s(is_high) = (d.data14.Natomsg_bounds(is_high,3) - N14(depth14index(is_high))')./dN14_high(is_high);
M14s(is_low) = (d.data14.Natomsg_bounds(is_low,3) - N14(depth14index(is_low))')./dN14_high(is_low);

% If there are Al-26 data, do this for them as well

if isfield(d,'data26')
    % Line up Al-26 data with predictions. 
    depth26_index = zeros(size(d.data26.z));
    for a = 1:length(d.data26.z)
        depth26index(a) = find(d.data26.z(a) == core_z_cm);
    end

    dN26 = d.data26.Natomsg_bounds(:,3) - d.data26.Natomsg_bounds(:,2);
    % Note: this is actually redundant with d.data26.dN, but keep it this way
    % to be more general

    % Calculate misfit to Al-26 data
    M26s = (d.data26.Natomsg_bounds(:,3) - N26(depth26index)')./dN26;
else
    M26s = [];
end
    


% Make chi-squared-like misfit
% Note: this overweights C-14, just because there are more data, but then
% again the uncertainties for C-14 are larger, so it is basically a wash. 
%
% Note: Al-26 data not included here if present

M = sum(M10s.^2) + sum(M14s.^2);





if dFlag == 1
    % Case not used as simple objective function, spit out a lot of data
    out.M = M;
    out.N10 = N10; 
    out.N14 = N14;
    out.miss10 = M10s;
    out.miss14 = M14s;
    out.core_z_cm = core_z_cm;
    out.core_z_gcm2 = core_z_gcm2;
    % Also do integration at plotting depths, return depth profile. 
    % Note: also returns predicted Al-26 profile for comparison; Al-26 not
    % used above in fitting algorithm  
    if isfield(d,'plot_z_cm')
        out.plot_z_gcm2 = d.plot_z_cm.*2.77;
        for a = 1:length(out.plot_z_gcm2)
            out.plot_N14(a) = integral(@(t) PofZofT(t,d.t1.times,d.t1.surface_elvs,(mass_depth_gcm2_v+out.plot_z_gcm2(a)),d.P,14).*exp(-d.l14.*t),0,maxt);
            out.plot_N10(a) = d.N10inh + integral(@(t) PofZofT(t,d.t1.times,d.t1.surface_elvs,(mass_depth_gcm2_v+out.plot_z_gcm2(a)),d.P,10).*exp(-d.l10.*t),0,maxt);
            out.plot_N26(a) = d.N26inh + integral(@(t) PofZofT(t,d.t1.times,d.t1.surface_elvs,(mass_depth_gcm2_v+out.plot_z_gcm2(a)),d.P,26).*exp(-d.l26.*t),0,maxt);
        end
    end  
else
    % Case used as objective function, just send M
    out = M;
end








