function CratModelFig7and8 %Reads .diam file and compares it to modeled surfaces from different ages clear all; close all; SECONDARY = true; %apply Smith's secondary crater model or not--true for fig 7 false for fig 8 EROSION = true; %whether or not to use erosion model for surface. should be on CUMULATIVE = false; %cumulative crater plots--false for fig 7 true for fig 8 SKIP_EMPTY = true; %skip empty bins with cumulative max_beta = 400; %max beta value to test, nm/a num_betas = 50; %number of betas to test, between 0 and max_beta prec = 0.25; %precision to consider crater count densities the same. 0.25 seems good. iters = 100; %iterations min_age = .1; %min age to be considered max_age = 4; %max age to be considered age_dev = 0.1; %space between ages to be considered num_ages = (max_age - min_age)/age_dev + 1; %number of ages to be considered area_surface = 10000; % large area with known age (i.e., true_age_surface), in km2 USE_COUNT_MIN_MAX = true; %whether or not to use the max/min bins from input .diam file %true for fig 7 false for fig 8 diameter_min = .088; %used if USE_COUNT_MIN_MAX is false diameter_max = 1000; %used if USE_COUNT_MIN_MAX is false %Creates age and beta matrices age_matrix = linspace(min_age,max_age,num_ages); beta_matrix = linspace(0,max_beta,num_betas); rng('shuffle'); %seeds random number generator with current time %% Opens diam file and reads craters [in_file,fpath]=uigetfile('*.diam','Pick a diam file!'); fid = fopen(strcat(fpath,in_file)); line = fgetl(fid); while line~=123&isnan(str2double(line)) %identifies start of crater list by '{' if strfind(line,'age') %finds age of parent surface surf_age = line((strfind(line,':')+1):end) end if strfind(line,'B') %finds beta surface was made with surf_beta = str2double(line((strfind(line,'=sur')+1):end)); end if strfind(line,'Large') %finds parent surface area surf_area = str2double(line((strfind(line,':')+1):end)); end if length(strfind(line,'area ='))>0||length(strfind(line,'Area'))>0 %finds subsurface area sub_surf = str2double(line((strfind(line,'=')+1):end)); end line = fgetl(fid); end if line==123 %my convention for start of craters, not needed if just start with numbers line = fgetl(fid); end i=1; while line~=125&line~=-1 %collects crater diameters until end of list '}' crater_str{i} = line; line = fgetl(fid); i=i+1; end fclose(fid); craters=str2double(crater_str); %converts file text to doubles diameters=craterfreq(1,0,1200); %sets diameter bins, from 0.0039-724 km %% Bins craters for i = 1:length(diameters)-1 %bins craters crater_count(i) = length(craters(craters>=diameters(i)&craters=diameters(end))); %last bin min_diam_idx = min(find(crater_count~=0)); %finds smallest bin max_diam_idx = max(find(crater_count~=0)); %finds largeest bin crater_count = crater_count(min_diam_idx:max_diam_idx); %takes smallest to largest present bins diameters = diameters(min_diam_idx:max_diam_idx); %takes smallest to largest present bins if ~USE_COUNT_MIN_MAX %if using user-defined crater sizes [~, idx_min] = min(abs(diameter_min-diameters)); %finds smallest and largest bins and defines diameters accordingly [~, idx_max] = min(abs(diameter_max-diameters)); diameters = diameters(idx_min:idx_max); crater_count = crater_count(idx_min:idx_max); end if CUMULATIVE %cumulates crates if needed compare_count = cumul(crater_count, SKIP_EMPTY); else compare_count = crater_count; end diameter_min = min(diameters); diameter_max = max(diameters); %% Creates and subsamples surfaces, and compares to inputted crater count fit_count = 0; for h = 1:num_ages disp('Percent done'); perc = (1-(num_ages-h)/num_ages)*100 %percent finished for k = 1:num_betas %crater density age = age_matrix(h); beta = beta_matrix(k); [test_diameters, test_density_craters] = craterfreq(age,diameter_min,... diameter_max,EROSION,SECONDARY,beta); %finds diameters and # of craters area_crater = pi*(test_diameters/2).^2; %computes crater areas %Calculates total number of craters after rounding num_craters_pre_round = area_surface.*test_density_craters; num_craters = round(num_craters_pre_round); num_bins = length(test_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)); %assigns random locations to craters 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 for b=1:iters test_crater_count = zeros(1,num_bins); %Defines a random box for this iteration box(b,1) = rand*(x_max-sqrt(sub_surf)); box(b,2) = rand*(y_max-sqrt(sub_surf)); %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)= box(b,1) && address_book(c,d,1)<= (box(b,1)+sqrt(sub_surf))... && address_book(c,d,2)>= box(b,2) && address_book(c,d,2)<= (box(b,2)+sqrt(sub_surf)) test_crater_count(c) = test_crater_count(c)+1; end end end end if CUMULATIVE %cumulates synthetic count if needed test_count = cumul(test_crater_count, SKIP_EMPTY); else test_count = test_crater_count; end perc_diff(h,k,b,:) = abs(test_count-compare_count)./((test_count+compare_count)./2); %calculates difference between .diam count and synthetic count perc_diff(h,k,b,isnan(perc_diff(h,k,b,:))) = 0; %fixes divide by zero errors when bins are empty if perc_diff(h,k,b,:) < prec %checks that all bins meet prec criteria fit_count = fit_count + 1; %counts number of fits fits(fit_count,1) = age; %stores age of fit fits(fit_count,2) = beta; %stores beta of fit fit_craters(fit_count,:) = test_crater_count; %stores synthetic count for displaying later end end clr_vars = {'test_diameters', 'test_density_craters','area_crater','num_craters_pre_round',... 'num_craters','num_bins','total_craters','max_craters_per_bin','address_book',... 'test_crater_count','RMSE_density'}; clear clr_vars{:}; end end %% Plots histogram f1 = figure(1); axes1 = axes('Parent',f1); edge{1} = [0,0.5:.5:max_age-.5,4.1]; %histogram edges edge{2} = [0:50:401]; hist3(fits,'edges',edge) xlim([0 4.1]) ylim([0 400.5]) xlabel('Ages (Ga)'); ylabel('\beta (nm a^{-1})'); zlabel('Modeled surfaces'); if CUMULATIVE title(strcat(in_file,'; ',num2str(sub_surf),' km^{2}; D_{min}:',... num2str(diameter_min), ' km; cumulative; iters: ',num2str(iters),'; prec:',num2str(prec))); else title(strcat(in_file,'; ',num2str(sub_surf),' km^{2}; D_{min}:',... num2str(diameter_min), ' km; not cumulative; iters: ',num2str(iters),'; prec:',num2str(prec))); end set(axes1,'FontSize',10,'FontWeight','bold'); %% Plots percent of good fits for each age bin f2 = figure(2); axes2 = axes('Parent',f2); h=hist3(fits,'edges',edge); %bins data bar([.25:.5:max_age-.25],sum(h(1:8,1:8)')./fit_count.*100,1) %plots data as sumsof beta bins, may need to adjust depending on bin sizes axis([0 4 0 100]); xlabel('Age (Ga)'); ylabel('Percent of matching surfaces'); yticks = get(gca,'ytick')'; % Gets the axis tick labels. percentsy = repmat('%', length(yticks),1); % makes a matrix yticklabel = [num2str(yticks) percentsy]; set(gca,'yticklabel',yticklabel) if CUMULATIVE title(strcat(in_file,'; ',num2str(sub_surf),' km^{2}; D_{min}:',... num2str(diameter_min), ' km; cumulative; iters: ',num2str(iters),... '; prec:',num2str(prec))); else title(strcat(in_file,'; ',num2str(sub_surf),' km^{2}; D_{min}:',... num2str(diameter_min), ' km; not cumulative; iters: ',num2str(iters)... ,'; prec:',num2str(prec))); end set(axes2,'FontSize',10,'FontWeight','bold'); %% Shows loaded and modeled crater counts. User can scroll throug them w/left and right keys. Space ends. buttona = 1; plot_num = 1; %which fitted count to display figure(3) while buttona ~= 32 %spacebar quits cla if CUMULATIVE scatter(diameters,cumul(fit_craters(plot_num,:)./sub_surf,SKIP_EMPTY),'DisplayName','Modeled') hold on scatter(diameters,compare_count./sub_surf,'d','DisplayName','Measured') ylabel('Cumulative crater count density (cumulative craters/km^2') else scatter(diameters,fit_craters(plot_num,:)./sub_surf,'DisplayName','Modeled') hold on scatter(diameters,compare_count./sub_surf,'d','DisplayName','Measured') ylabel('Crater count density (craters/km^2') end set(gca, 'YScale', 'log') set(gca, 'XScale', 'log') xlabel('Crater Diameter (km)') axis([10^-3 10^3 10^-6 10^3]) legend('show') title(strcat(num2str(fits(plot_num,1)),' Ga; \beta=',... num2str(fits(plot_num,2)),'; Fit ',num2str(plot_num),'/',num2str(fit_count))); [x,y,buttona] = ginput(1); if buttona == 29 %right key if plot_num1 plot_num = plot_num-1; end end end 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