function [t_array,nN_array] = naturalDRC_GOK(final_time,D_dot,rhoPrime,D_0,cVal)
%   N. Brown, 4/22/21
%   Simulates n/N (r') buildup over geologic time following GOK

% final_time in s
% D_dot in Gy/ka

%dimensionless lattice distances considered, m
rP_min=0.01;
rP_step=0.01;
rP_max=2.5;
rP_array=[rP_min:rP_step:rP_max]';    
    
%empty traps to begin with
nN_rP_0=zeros(cast(((rP_max-rP_min)/rP_step)+1,'uint8'),1);

    s=3e15;     

    %terms
    alpha=cVal+1;
    
    Ds=(D_dot/(1e3*365.25*24*3600))/D_0;%(Gy/ka/Gy to Gy/s/Gy) to 1/s
    P_g=s*exp(-rhoPrime^(-1/3)*rP_array);
        
    %remove dimensionality for GOK approach
    function [dnN_array] = dnN(t,nN_rP)
        %integrate for every r val
        term1=Ds*(1-nN_rP).^alpha;
        term2=P_g.*nN_rP;         
        dnN_array = term1 - term2;
    end

%execute solver
[t_array,nN_rP_array]=ode23tb(@dnN,[0 final_time],nN_rP_0);
    
nN_array=rP_step*sum(nN_rP_array.*[3*rP_array.^2.*exp(-(rP_array).^3)]',2);
end