function [age,varargout]=Nye_interp(z,ACP,H,varargin)
% Interplates ice-ages at depths z, given the Age Control Points (ACP) and core
% depth (H)

% [age,varargout]=Nye_interp(z,ACP,H,...)

% Inputs
%     z: depths to interpolate ages for
%     ACP: N x 2 matrix of Age Control Points ([Depths, Age])
%     H: Full Length of core under strain conditions

% Varargin
%     .delta_age: if this argument is used, ACPs are input as gas ages and
%         the argument out becomes gas age.
%     .ref_date: Define 'Present' (default = 2000 CE)
%     .date0: Age of ice at 0 m depth (date core drilled?)

% Outputs
%     age: age of ice at depth 'z'
%     b_dot: accumulation history (ice age);


% Written by James E. Lee, College of Earth, Ocean, Atmospheric Sciences,
% Oregon State University
% Last Update: 03/16/15
%% Deal with input format
if ~isnumeric(z) || ~isnumeric(ACP) || ~isnumeric(H)
    error('Missing required input (z,ACP,H)')
end
  
%% Formating varargin fields
if rem(size(varargin,2),2)==0
    fields=reshape(varargin,2,size(varargin,2)/2);
    fields=cell2struct(fields(2,:),fields(1,:),2);
else
    fprintf('Input fields must enter as pairs')
end

%% Preset Values for optional inputs
delta_age=0;
ref_date=2000; %sets reference date
date0=2013; %age of depth at 0 depth

if isfield(fields,'delta_age');delta_age=fields.delta_age;end
if isfield(fields,'ref_date');ref_date=fields.ref_date;end
if isfield(fields,'date0');date0=fields.date0;end
  
%% Deal with z formatting
if ~isvector(z)
    error('Error: z must be vector')
end
if max(z)>=H;
    fprintf('Input z abbreviated to depths less than H');
    z=z(z<H);
end
z=reshape(z,[],1); %Reshapes z into column vector;

%% Deal with ACP formatting

ACP(:,2)=ACP(:,2)+delta_age; %convert to ice age from ACP
if ACP(1,1); %adds ACP at 0 m depth
    ACP=[0 ref_date-date0;ACP];
end

id=isnan(sum(ACP,2));ACP=ACP(~id,:);
if ~issorted(ACP(:,1));ACP=sortrows(ACP,1);end
if ~issorted(ACP(:,2));error('Error: ACP not in stratigraphic order');end

Z=ACP(:,1); %Age control depths
Age=ACP(:,2); %Age control ages
%% Calculate stuff

% mean accumulation rate between ACP's;
b_dot=H*log((H-Z(1:end-1))./(H-Z(2:end)))./diff(Age);
b_dot(end+1)=b_dot(end); %Assume accumulation rate remains constant deeper than last ACP

% solve for ages
j=sum(z*ones(size(Z'))>ones(size(z))*Z',2); %Index of ACP's shallower/younger than elements in z
age=Age(j)+(H./b_dot(j)).*log((H-Z(j))./(H-z)); 

age=age-delta_age;

%% Varargout
varargout={b_dot};