function CratModelFig2and5
    %Plots calculated subsurface ages given a parent age, surface area, and Beta value

    %% initialize variables and input parameters
    clear all; close all; 

    min_age = 1; %min age to be considered
    max_age = 3;  %max age to be considered
    num_ages = 3; %number of ages to be considered
    area_surface = 20000; % large area with known age (i.e., true_age_surface), in km2. 200000 used for subsurface areas greater than 10 km
    sub_area = 100;       %supsampled area, km
    iters = 100;  %number of iterations
    diameter_max = 1000; % largest crater diameter on area, in km (max = 1024 km)
    EROSION = true;  %if erosion is modeled Based on Smith
    SECONDARY = false;  %if secondary craters are modeled based on SMith
    CUMULATIVE = true; %counts and plots cumulatively
    SKIP_EMPTY = true; %skips bins with no craters when using cumulative
    USE_QUANTILE = true;  %use interquartile range or mean/std dev to invert?
    beta = 10; %erosion rate in nm/a
    %how to weight RMSE (0=no weight, 1=/#craters, 2=/#craters^2.
    %4=logarithmically
    weight = 2;
    FIT_EROSION = true;  %If fitting counted crater densities uses Smith model to include erosion
    AVG_NO_CRATERS = false;  %if should average subareas with no craters or ignore them (have 0.1 Ga age)

    %Creates age and min crater size matrices
    age_matrix = linspace(min_age,max_age,num_ages);
    num_ages = length(age_matrix);
    min_crat_matrix = [.016 .063 .125 .178 .25 .375 .5 1]; %minimum crater diamters, km, smaller sizes omitted for subsurface areas greater than 10 km
    
    rng('shuffle');             %seeds random number generator with current time

    %% Calculates average ages
    for h = 1:num_ages
        disp('Percent done');
        perc = (1-(num_ages-h)/num_ages)*100  %percent finished
        %crater density
        age = age_matrix(h);
        [diameters, density_craters] = craterfreq(age,min_crat_matrix(1),...
            diameter_max,EROSION,SECONDARY,beta);  %finds diameters and # of craters

        area_crater = pi*(diameters/2).^2; %computes crater areas

        %Calculates total number of craters after rounding
        num_craters_pre_round = area_surface.*density_craters;
        num_craters = round(num_craters_pre_round);
        num_bins = length(diameters);
        total_craters = sum(num_craters);
        max_craters_per_bin = max(num_craters);

        % create matrix of (crater diameters bins, each crater in the bin, x&y) 
        % create coordinates for edges of surface
        x_length = sqrt(area_surface);
        y_length = sqrt(area_surface);

        x_min = 0;
        y_min = 0;
        x_max = sqrt(area_surface);
        y_max = sqrt(area_surface);

        address_book = nan(num_bins,max_craters_per_bin,2);

        for i=1:num_bins
            if num_craters(i)>0   %checks if any craters are in this diameter bin
                address_book(i,1:num_craters(i),1) = x_length.*rand(1,num_craters(i));
                address_book(i,1:num_craters(i),2) = y_length.*rand(1,num_craters(i));
            end
        end

        %Counts craters located within sub-sample areas (each sub-sampled area
        % placed on the board randomly and number of sub-samples determined by
        % "iters"

        %counts number of craters for each sub area, for each iteration, in each bin
        crater_count = zeros(length(sub_area),iters,num_bins);

        for a = 1:length(sub_area)
            for b = 1:iters
                %Defines a random box for this iteration
                box(a,b,1) = rand*(x_max-sqrt(sub_area(a)));
                box(a,b,2) = rand*(y_max-sqrt(sub_area(a))); %coordinates of subsurface a

                %Checks each crater to see if it's in this interation's box
                
                for c = 1:num_bins      %each diameter bin size
                    if num_craters(c)>0 && area_crater(c)<sub_area(a)    %If there are craters in bin smaller than subsurface
                        for d=1:num_craters(c)     %for each crater in the bin
                            if address_book(c,d,1)>= box(a,b,1) && address_book(c,d,1)<= (box(a,b,1)+sqrt(sub_area(a)))...
                                  && address_book(c,d,2)>= box(a,b,2) && address_book(c,d,2)<= (box(a,b,2)+sqrt(sub_area(a)))
                                crater_count(a,b,c) = crater_count(a,b,c)+1;
                            end
                        end
                    end
                    crater_count_density(a,b,c) = crater_count(a,b,c)/sub_area(a);
                end
            end
        end
        %Calculates RMSE compared to isochrons with no erosion
        for k = 1:length(min_crat_matrix)    
            [~,idx] = min(abs(diameters-min_crat_matrix(k))); %bin for min crater diameter
            for i = 1:length(sub_area)   %subsurface areas 
                for j = 1:iters  %iterations
                    if CUMULATIVE %cumulates crater densities if desired
                        RMSE_density = cumul(crater_count_density(i,j,:), SKIP_EMPTY);
                    else
                        RMSE_density = crater_count_density(i,j,:);
                    end

                    [age_calc(h,i,j,k), beta_calc(h,i,j,k), RMSE_min(h,i,j,k)] =...
                        craterfit(diameters(idx:end),RMSE_density(idx:end),...
                        weight,CUMULATIVE,SKIP_EMPTY,FIT_EROSION,SECONDARY); 
                end
                %calculates averages and standard deviations over all iterations
                if AVG_NO_CRATERS            
                    age_to_avg = age_calc(h,i,:,k);  %includes all ages, even if no craters (0.1 Ga)
                    beta_to_avg = beta_calc(h,i,:,k);
                else
                    age_to_avg = age_calc(h,i,age_calc(h,i,:,k)>0.1,k);  %only includes ages with craters (>0.1 Ga)
                    beta_to_avg = beta_calc(h,i,beta_calc(h,i,:,k)>0.1,k);
                end
                age_avg(h,i,k) = mean(age_to_avg);  %calculates mean and std dev to plot
                age_std(h,i,k) = std(age_to_avg);
                beta_avg(h,i,k) = mean(beta_to_avg);
                beta_std(h,i,k) = std(beta_to_avg);
                age_med(h,i,k) = median(age_to_avg);  %calculates median and interquantile to plot
                age_25(h,i,k) = quantile(age_to_avg,.25);
                age_75(h,i,k) = quantile(age_to_avg,.75);
            end    
        end
        clr_vars = {'diameters', 'density_craters','area_crater','num_craters_pre_round',...
            'num_craters','num_bins','total_craters','max_craters_per_bin','address_book',...
            'crater_count','RMSE_density'};
        clear clr_vars{:};
    end
    
    %% plots lines
    SymbOrb={'o', 'x','h','s','d', '^','<','>','x', 'o', '+', '.', 'p','h','s','d', '^','<','>','x', 'o', '+', '.', 'p','h','s','d', '^','<','>','x', 'o', '+', '.', 'p','h','s','d', '^','<','>','x', 'o', '+', '.', 'p','h','s','d', '^','<','>','x', 'o', '+', '.', 'p','h','s','d', '^','<','>','x', 'o', '+', '.', 'p','h','s','d', '^','<','>','x', 'o', '+', '.', 'p','h','s','d', '^','<','>','x', 'o', '+', '.', 'p','h','s','d', '^','<','>','x', 'o', '+', '.', 'p','h','s','d', '^','<','>','x', 'o', '+', '.', 'p','h','s','d', '^','<','>','x', 'o', '+', '.', 'p','h','s','d', '^','<','>','x', 'o', '+', '.', 'p','h','s','d', '^','<','>','x', 'o', '+', '.', 'p','h','s','d', '^','<','>','x', 'o', '+', '.', 'p','h','s','d', '^','<','>','x', 'o', '+', '.', 'p','h','s','d', '^','<','>','x', 'o', '+', '.', 'p','h','s','d', '^','<','>','x', 'o', '+', '.', 'p','h','s','d', '^','<','>','x', 'o', '+', '.', 'p','h','s','d', '^','<','>','x', 'o', '+', '.', 'p','h','s','d', '^','<','>','x', 'o', '+', '.', 'p','h','s','d', '^','<','>','x', 'o', '+', '.', 'p','h','s','d', '^','<','>','x', 'o', '+', '.', 'p','h','s','d', '^','<','>','x', 'o', '+', '.', 'p','h','s','d', '^','<','>','x', 'o', '+', '.', 'p','h','s','d', '^','<','>','x', 'o', '+', '.', 'p','h','s','d', '^','<','>','x', 'o', '+', '.', 'p','h','s','d', '^','<','>','x', 'o', '+', '.', 'p','h','s','d', '^','<','>','x', 'o', '+', '.', 'p','h','s','d', '^','<','>','x', 'o', '+', '.', 'p','h','s','d', '^','<','>','x', 'o', '+', '.', 'p','h','s','d', '^','<','>','x', 'o', '+', '.', 'p','h','s','d', '^','<','>','x', 'o', '+', '.', 'p','h','s','d', '^','<','>'};
    
    f = figure(1);

    for i = 1:length(age_matrix)
        med=squeeze(age_med(i,1,:))';
        quant25=med-squeeze(age_25(i,1,:))';
        quant75=squeeze(age_75(i,1,:))'-med;
        errorbar(min_crat_matrix.*1000,med,quant25,quant75,...
            SymbOrb{i})
        set(gca,'xscale','log')
        hold on
    end
    ylim([0 4]);   %formats graph
    xlim([10 1200]);
    xlabel('Minimum crater diameter (m)');
    ylabel('Median calculated age (Ga)');
    title(strcat(num2str(sub_area(1)),' km^2'));
    ax=gca;
    ax.XTick=[16 63 125 250 500 1000];
    set(gca,'YMinorTick','on');
    hold off
end

function [CF_crater_D, CF_crater_N] = craterfreq(CF_age, CF_diameter_min,...
    CF_diameter_max, CF_EROSION, CF_SECONDARY, CF_Beta)
    %craterfreq: Based on range of craters chosen and surface age (Ga)
    %returns list of crater frequencies for 1 Ga surface age (no. craters/km2).
    %If EROSION is true uses Smith et al 2008 model for erosion and B is
    %erosion rate (nm/a)

    %if false, uses production function from Smith 2008
    %if true, uses hartmann 1 Ga density, times ratio, divided by age
    CF_HARTMANN_PROD = true; 


    switch nargin   %allows for fewer inputs to function, giving default values below
        case 0
            CF_age = 1;
            CF_diamter_min = 0.0039;
            CF_diameter_max = 1024;
            CF_EROSION = false;  
        case 1
            CF_diamter_min = 0.0039;
            CF_diameter_max = 1024;
            CF_EROSION = false;
        case 2
            CF_diameter_max = 1024;
            CF_EROSION = false;
        case 3
            CF_EROSION = false;
        case 4
            CF_SECONDARY = false;
            CF_Beta = 1e-5;
        case 5
            CF_Beta = 1e-5;
    end
    %% Calculates (with erosion) or loads crater sizes and frequencies
     %crater diameter bins, Michael Icarus 2013
    CF_D = [0.00391 0.00553 0.00782 0.0111 0.01565 0.0221 0.0313 0.0442 0.06251...
            0.0887 0.125 0.177 0.25 0.354 0.5 0.7075 1 1.415 2 2.83...
            4 5.66 8.05 11.32 16.05 22.63 32.05 45.3 64.05 90.6 128.05 181.1...
            256.05 362.1 512.05 724.1];  %km   
    % Crater frequencies for 1 Ga surface, Michael Icarus 2013
    CF_Nh = [4.04E+03 2.33E+03 1.14E+03 4.58E+02 1.91E+02 6.66E+01 2.40E+01...
            9.44E+00 3.30E+00 1.22E+00 4.37E-01 1.47E-01 4.70E-02 1.38E-02 ...
            4.02E-03 1.15E-03 3.08E-04 1.28E-04 6.85E-05 3.67E-05 1.98E-05 ...
            1.06E-05 5.68E-06 3.04E-06 1.62E-06 8.71E-07 4.67E-07 2.40E-07 ...
            1.12E-07 5.21E-08 2.43E-08 1.13E-08 5.28E-09 2.47E-09 1.15E-09 ...
            5.37E-10];

    CF_N_1Ga = 3.79e-14*(exp(6.93)-1)+5.84e-4;   %Michael 2013 Icarus, pg 889, eqn 3
    CF_ratio = (3.79e-14*(exp(6.93*CF_age)-1)+5.84e-4*CF_age)/CF_N_1Ga;  %ratio to 1 Ga

    if CF_EROSION
        if CF_Beta == 0
            CF_Beta = 1e-5;  %minimum erosion rate
        end
        for i = 1:length(CF_D)
            CF_p=0;, CF_L=0;, CF_Xi=0;, CF_Psi=0;

            if ~CF_HARTMANN_PROD
               if CF_D(i) < 1.4  %sets production fnctn p
                    CF_p = 0.0035*(0.13*log(CF_D(i))+0.83)/(CF_D(i)^3.3);
                else if CF_D(i) <= 48.1
                        CF_p = 10^(-1.8*log10(CF_D(i))-2.59);
                    else
                        CF_p = 10^(-2.2*log10(CF_D(i))-1.89);
                    end
                end
                CF_p = CF_p * 0.29;    
            else
                CF_p = CF_Nh(i)*CF_ratio/CF_age;
            end


            if CF_D(i) < 5.8 %sets Xi   
                CF_Xi = 0.2*CF_D(i);
            else
                CF_Xi = 0.42*log(CF_D(i))-0.01;
            end

            if CF_D(i) < 1.2 && CF_SECONDARY  %sets Psi
                CF_Psi = 1 + 0.1225/(0.1225+CF_D(i)^2);  %yields eqn 10
            else
                CF_Psi = 1;
            end

            CF_L = CF_Psi*CF_Beta/(1000*CF_Xi);  %lambda: crater loss fnctn

            CF_N(i) = CF_p/CF_L*(1-exp(-CF_L*CF_age));
        end
    else
        CF_N = CF_Nh .* CF_ratio;
    end
    %% Crater vectors of crater sizes and frequencies based on desired crater size range

    % find closest value in D to diameter_min
    [~, CF_idx_min] = min(abs(CF_D-CF_diameter_min));

    % find closest value in D to diameter_max
    [~, CF_idx_max] = min(abs(CF_D-CF_diameter_max));

    % create arrays for diameter range and crater frequency range

    CF_crater_D = CF_D(CF_idx_min:CF_idx_max);
    CF_crater_N = CF_N(CF_idx_min:CF_idx_max);
end

function [C_cumulate_count_density] = cumul(C_crater_count_density, C_IGNORE_EMPTY)
%%Cumulates crater count densities

    switch nargin
        case 1
            C_IGNORE_EMPTY = true;  %skips bins with no craters
    end

    C_num_bins = length(C_crater_count_density);    
    for C_k = 1:C_num_bins   %size bins
        C_cumulate_count_density(C_k) = 0;   
        if C_crater_count_density(C_k) > 0 || ~C_IGNORE_EMPTY
            for C_m = 1:(C_num_bins - C_k + 1)
                C_m_rev = C_num_bins-C_m+1;
                C_cumulate_count_density(C_k) = C_cumulate_count_density(C_k)...
                  + C_crater_count_density(C_m_rev);
            end
        end    
    end
end

function [FIT_age_calc, FIT_beta_calc, FIT_RMSE_min] = craterfit(FIT_diameters,...
    FIT_count_dens, FIT_weight, FIT_CUMULATIVE,FIT_SKIP_EMPTY,FIT_FIT_EROSION,FIT_SECONDARY)
    
    %% Provides best fit age and beta (if FIT_EROSION) and RMSE for a crater count
    switch nargin
        case 5
            FIT_FIT_EROSION = false;
            FIT_SECONDARY = false;
        case 6
            FIT_SECONDARY = false;
    end

    FIT_prec = 10;  %precision for ages--higher slows runtime but gives more precise ages

    if FIT_FIT_EROSION  %Test betas if FIT_EROSION is on
        FIT_test_b = [1 10:10:400];
    else
        FIT_test_b = [1];
    end
    
    if ~FIT_CUMULATIVE
        FIT_max_d = FIT_diameters(end);   %just the diameters measured
    else
        FIT_max_d = 1000;        %all craters larger than measured, since need to cumulate them
    end

    for FIT_i = 1:length(FIT_test_b)  %betas
        for FIT_k = 1:4*FIT_prec %ages
            FIT_RMSE(FIT_i,FIT_k)=0;
            if FIT_FIT_EROSION   %gets theoretical crater frequency
                [~, FIT_d_theory] = craterfreq(FIT_k/FIT_prec,FIT_diameters(1),...
                    FIT_max_d,true,FIT_SECONDARY,FIT_test_b(FIT_i));
            else
                [~, FIT_d_theory] = craterfreq(FIT_k/FIT_prec,FIT_diameters(1),FIT_max_d);
            end
            if FIT_CUMULATIVE  %Cumulates crater frequnecy if desired
                FIT_d_theory = cumul(FIT_d_theory, FIT_SKIP_EMPTY);
            end
            for FIT_m = 1:length(FIT_diameters)
                if FIT_weight==4  %gets RMSE for log weighting
                    if FIT_count_dens(FIT_m)>0
                        FIT_RMSE(FIT_i,FIT_k) = FIT_RMSE(FIT_i,FIT_k)+...
                            (log10(FIT_count_dens(FIT_m))...
                            -log10(FIT_d_theory(FIT_m))).^2;
                    end
                else
                    if ~isnan(FIT_count_dens(FIT_m))  %gets RMSE for non-log weighting
                        if ~FIT_SKIP_EMPTY||FIT_count_dens(FIT_m)~=0
                            FIT_RMSE(FIT_i,FIT_k) = FIT_RMSE(FIT_i,FIT_k)+(FIT_count_dens(FIT_m)...
                            - FIT_d_theory(FIT_m)).^2/FIT_d_theory(FIT_m).^(FIT_weight-1);
                        end
                    end
                end
            end
            FIT_RMSE(FIT_i,FIT_k) = FIT_RMSE(FIT_i,FIT_k)^.5;    
        end
    end
    if FIT_FIT_EROSION  %finds minimum RMSE, depending on whether FIT_EROSION, which requires checking multiple columns
        FIT_RMSE_min = min(min(FIT_RMSE));
        [~, FIT_idx_k] = min(min(FIT_RMSE));
        [~, FIT_idx_i] = min(FIT_RMSE(:,FIT_idx_k));
        FIT_age_calc = FIT_idx_k/FIT_prec;
        FIT_beta_calc = FIT_test_b(FIT_idx_i);
    else
        FIT_RMSE_min = min(FIT_RMSE);
        [~, FIT_idx] = min(FIT_RMSE);
        FIT_age_calc = FIT_idx/FIT_prec;
        FIT_beta_calc = 0;
    end
end