function [age_MCarray,age2D0]=main_MC(DRCdata,fadeData,labDoseRate,geolDoseRate,geolDoseRateErr,titleString,numMCruns)
% N. Brown, 4/22/21
% INPUTS
%   DRCdata: m x 4 matrix. Columns are: dose (s beta), Li/Ti, Li/Ti errors,
%       effective delay times of dose points (s)
%   fadeData: m x 3 matrix. Columns are: Li/Ti, Li/Ti errors, effective
%      delay time of fade points (s)
%   labDoseRate: at grain (Gy/s)
%   geolDoseRate: (Gy/ka)
%   titleString: used to label plots
%   numMCruns: for error analysis, number of Monte Carlo iterations
%
% OUTPUTS
%   age_MCarray: from which mean and std dev are calculated, based on
%      'numMCruns' iterations
%   age2D0: resolvable upper age limit (2*D0/Ddot)

    %%  Fit fading data to Kars et al. (2008) expression to get rho' value
    %       Also plots result.

    % First, find which is shortest duration fade point(s). Sometimes there are multiple points
    %    at approximately or exactly the same fade duration (e.g., the recycled
    %    point within a DRC) which should be averaged to get the initial
    %    intensity. Assume these effective fade times are within 2% of each
    %    other.
    figure('Position',[679,79,414,650])
    subplot(2,1,1)
    
    sortedFadeData=sortrows(fadeData,3);
    L0T0inds=find(sortedFadeData(:,3)<1.02*sortedFadeData(1,3));

    s=3e15;%s^-1, Huntley (2006)
    I0_guess=3;%starting guess for I_0 within Eq. 5 of K+08
    rhoP_guess=1e-6;
    [~,rhoPrime,fitParamSigmas]=rhoPrimeFitter(I0_guess,rhoP_guess,fadeData(:,1)./mean(sortedFadeData(L0T0inds,1)),fadeData(:,2)./fadeData(:,1),fadeData(:,3),s);
    rhoPrimeErr=fitParamSigmas(2);

    title(titleString)
    
    %%  Plot dose response data, both as measured and 'unfaded'; rho' errors not propogated
    subplot(2,1,2)
    
    z=1.8;%e.g., Jain et al., (2012) p. 6

    % unfade the non-zero-dose points; errors are simply borrowed
    DRC_uf=DRCdata(:,1:3);%copy 'as measured' DRC to start with
    nonzeroInds=find(DRCdata(:,1)~=0);

    %%  Fit the unfaded data to a GOK growth curve (Guralnik et al., 2015) to get unfaded D0 value
    kineticOrderTrapping=2;%first-order ~ 1.001 (1 undefined), second-order = 2; serves as upper limit within solver
    D0_MCarray=zeros(numMCruns,1);
    LTmax_MCarray=zeros(numMCruns,1);
    cuf_MCarray=zeros(numMCruns,1);
    
    doseArrayDRCuf=logspace(-1,log10(max(labDoseRate*DRC_uf(:,1))),1e3);
    
    h=waitbar(0,'Finding unfaded D_0...');
    for i=1:numMCruns
        waitbar(i/numMCruns)
        rhoPrime_i=normrnd(rhoPrime,rhoPrimeErr);
        if rhoPrime_i<0
            rhoPrime_i=1e-8;
        end
        DRC_uf(nonzeroInds,2)=DRCdata(nonzeroInds,2)./exp(-rhoPrime_i*log(z*s*DRCdata(nonzeroInds,4)).^3);
        [paramsGOK,~]=DRCfitter_GOK(DRC_uf(:,2),DRC_uf(:,3),labDoseRate*DRCdata(:,1),kineticOrderTrapping);
        LTmax_MCarray(i)=paramsGOK(1);
        D0_MCarray(i)=paramsGOK(2);
        cuf_MCarray(i)=paramsGOK(3);
        
        %plot each curve
        fitArrayDRCuf=paramsGOK(1)*(1-(1+paramsGOK(3)*doseArrayDRCuf./paramsGOK(2)).^(-1/paramsGOK(3)));    
        plot(doseArrayDRCuf,fitArrayDRCuf,'color',[0.9,0.9,0.9,0.05])        
        hold on
    end
    close(h)
    
    D0uf=mean(D0_MCarray);
    D0ufErr=std(D0_MCarray);
    LTmaxuf=mean(LTmax_MCarray);
    LTmaxErr=std(LTmax_MCarray);
    cuf=mean(cuf_MCarray);
    cufErr=std(cuf_MCarray);
    
    age2D0=2*D0uf/geolDoseRate;
    
    % plot as measured data
    xLL=3e1;
    xUL=3e3;
    errorbar(labDoseRate*DRCdata(2:end,1),DRCdata(2:end,2),DRCdata(2:end,3),'ob')
    hold on
    plot([xLL,xUL],[DRCdata(1,2),DRCdata(1,2)],'r')
    plot([xLL,xUL],[DRCdata(1,2)+DRCdata(1,3),DRCdata(1,2)+DRCdata(1,3)],'--r')    
    plot([xLL,xUL],[DRCdata(1,2)-DRCdata(1,3),DRCdata(1,2)-DRCdata(1,3)],'--r')    
    
    text(min(labDoseRate*DRC_uf(nonzeroInds,1)),0.9*max(DRC_uf(:,2)),['D_0: ' num2str(D0uf,'%.f') char(177) num2str(D0ufErr,'%.f') ' Gy' char(10)...
        'c: ' num2str(cuf,'%.1f') char(177) num2str(cufErr,'%.1f') char(10)...
        'LT_{max}: ' num2str(LTmaxuf,'%.1f') char(177) num2str(LTmaxErr,'%.1f')])
       
    set(gca,'yscale','log')
    set(gca,'xscale','log')
    xlim([xLL xUL])
    ylim([0.1 1e3])

    %%  Simulate n(r') evolution to produce natural dose response points ('fully faded')
    final_time_s=(max(labDoseRate*DRCdata(:,1))/geolDoseRate)*1e3*365.25*24*3600;
    
    De_MCarray=zeros(numMCruns,1);
    age_MCarray=zeros(numMCruns,1);
    
    h=waitbar(0,'Simulating natural dose response...');
    i=1;
    while i<=numMCruns
        waitbar(i/numMCruns)
        
        %get randomized inputs
        rhoPrime_i=normrnd(rhoPrime,rhoPrimeErr);
        if rhoPrime_i<0
            rhoPrime_i=1e-8;
        end
        D0ufi=normrnd(D0uf,D0ufErr);
        cufi=normrnd(cuf,cufErr);
        if cufi>(kineticOrderTrapping-1)
            cufi=kineticOrderTrapping-1;
        elseif cufi<=0
            cufi=0.001;            
        end       
        geolDoseRatei=normrnd(geolDoseRate,geolDoseRateErr);
        LTmaxufi=normrnd(LTmaxuf,LTmaxErr);
        LNTNi=normrnd(DRCdata(1,2),DRCdata(1,3));
        
        %simulate natural history
        [natDRCt_arrayi,natDRC_nN_arrayi] = naturalDRC_GOK(final_time_s,geolDoseRatei,rhoPrime_i,D0ufi,cufi);
        
%         %plot results
%         plot(geolDoseRatei/(1e3*365.25*24*3600)*natDRCt_arrayi,LTmaxufi*natDRC_nN_arrayi,'color',[1,0,0,0.05],'.')
         hold on
        
        %fit results to GOK expression
        [paramsGOKnati,~]=DRCfitter_GOK(LTmaxufi*natDRC_nN_arrayi,0.1*ones(size(LTmaxufi*natDRC_nN_arrayi)),geolDoseRatei/(1e3*365.25*24*3600)*natDRCt_arrayi,kineticOrderTrapping);
        LTmaxNati=paramsGOKnati(1);
        D0nati=paramsGOKnati(2);
        cNati=paramsGOKnati(3);
        
        fitArrayDRCnati=LTmaxNati*(1-(1+cNati*doseArrayDRCuf./D0nati).^(-1/cNati));
        plot(doseArrayDRCuf,fitArrayDRCnati,'color',[1,0,0,0.05]) 
        
        % get De value
        Dei=((-LNTNi/LTmaxNati + 1)^-cNati - 1)*(D0nati/cNati);
        if isreal(Dei)&& Dei > 0
            De_MCarray(i)=Dei;
            age_MCarray(i)=De_MCarray(i)/geolDoseRatei;
            i=i+1;
        end
    end
    close(h)
    
    De=mean(De_MCarray);
    DeErr=std(De_MCarray);
    age=mean(age_MCarray);
    ageErr=std(age_MCarray);
    
    ylim([1e-1 1.05*max(DRC_uf(:,2))])
    xlabel('Dose (Gy)')
    ylabel('L_x/T_x')   
    
    text(0.5*labDoseRate*max(DRC_uf(:,1)),0.2*LTmaxuf,['D_e: ' num2str(De,'%.f') char(177) num2str(DeErr,'%.f') ' Gy' char(10)...
        'Age: ' num2str(age,'%.1f') char(177) num2str(ageErr,'%.1f') ' ka'])    
end