%BAYES_SEDIMENT.m
%Edwin Kite
%
%Purpose:      Determine Mars sedimentation rate using frequency of
%              embedded craters.
%Description: 
%Inputs:       A list or lists of crater diameters. 
%Outputs:     
%Assumptions: (i)   The landscape is initially horizontal, and aggrades vertically at a uniform rate Ws
%             (ii)	Cratering is a Poisson process with initial crater
%             profiles (including rim) given by Melosh, and the impact flux
%             doesn't change with time.
%             (iii)	When the level of the landscape exceeds the level of the rim, the crater is immediately filled to the brim.
%             (iv)	Erosion does not preferentially expose craters (reword but retain this point)
%          
%Notes:       

obliterationdepthfraction = 0.2; %From Melosh, 1989 (approximately)

%**********OBSERVATIONS TO BE FIT ***************

countarea = (2302.143-178.01)*1e12; %count area in m^2 <-- area of channel belts minus area of holes

%From circle fitting to preserved part of arc: 
sizes = [  387.7016
   56.4344
   70.4752
   68.8069
   81.3375
  118.8421
  131.6533
  118.5345
  146.2326
  154.3903
  176.5247
  175.1765
  157.8368
  495.5084
  346.2190
  237.3461
  187.5713]*2/1e3; %diameter in km

%From circle fitting to preserved part of arc - candidates: 
csizes = 1/1e3*2*[  229.2447
  557.9057
  348.4017
  339.3956
  477.2053
  163.3322
  209.0119
  117.7217
   98.7535
  212.4212
  244.8331
  269.7765
   95.8815
   38.7171
   57.2853
  197.3651
   50.3188
  105.9558
   98.9631
   63.5207
   90.4993
  233.6493
  104.4690
  203.0253
  189.7922
  177.6827
  248.5467
  164.6746
  316.4349
  178.1513
  219.7324
   48.8328
  125.3092
   73.4154
   77.2019
  214.2061
   76.0797
   89.8379
  195.8929
  141.1061
  163.1332
  103.0620
  152.9041];


%*****************CHRONOLOGY AND CRATER-PRODUCTION MODELS**********************

wtt = ... %werner_tanaka_icarus2011_table1
[11  -0.00070349 0.00005805492 0;  ...
10 0.0023664     0.0006232845 0; ...
9 0.0108278     -0.004753462 0; ...
8 -0.04034      -0.01180639 0; ...
7 -0.0568509     0.06755923 0; ...
6 0.256181       0.1015683 0;  ...
5 0.109792      -0.3630098 0;  ...
4 -0.7735       -0.4860814 0;  ...
3 -0.0111926     0.7915374 0;  ...
2 1.25281        1.256814 0; ...
1 -2.65407      -3.197453  -2; ...
0 -3.1665       -3.383677 -3.383677]'; %Production functions. 1 Hartmann, 2 Ivanov, 3 'minus two'

%Martian chronology functions 
%1 Ivanov, 2001; 
%2 Hartmann & Neukum, 2001; 
%3 Hartmann, 2005]

%NOTE CHANGE PER FEB 2013 EMAIL FROM CALEB FASSETT:
%"[...] Hartmann 2005 chronology function.  It turns out that Greg Michael and I have recently realized
%that the chronology function in WT2011 is mistaken, and that the correct Hartmann chronology function 
%is something like:
%N(1)=3.79E-14*(exp(6.93t)-1)+5.84E-04t  (Greg's calculation)
%N(1)=3.86E-14*(exp(6.93t)-1)+5.92E-04t   (What I've been using since 2006)

mcf = [  2.68e-14,   4.1e-4; ... 
         3.22e-14,   4.875e-4; ...
         3.79e-14,   5.84e-4];
   
 % 4.4246e-14,   6.8158e-4]; %DEPRECATED per Caleb Fassett email
%  
   
   
   clear ncum
   D = logspace(-3,1,61);%logspace(-3,3,22);
   ncum = zeros(3,length(D));
for mpf = 1:3
     for coeff = 1:12
         ncum(mpf,:) = ncum(mpf,:) + wtt(mpf+1,coeff).*log10(D).^wtt(1,coeff);
     end
         %ncum(mpf,:) = 10.^ncum(mpf,:);%Leave as linear for now
end

age = [0:0.01:4.6];

for l = 1:3 %mcf %These ncr agree with WernerTanaka2011 Table 2 (sometimes disagree in last decimal place)
    ncr(l,:) = mcf(l,1).*(exp(6.93.*age) - 1) + mcf(l,2).*age;
end

dncr = diff(ncr')'; 

%Make isochron diagram (for testing)

%Deprecated
ats = [1e-3,1,3,3.5,4.0,4.3];%ages to show


aa = 0;

for aoi = [3.7, 3.9]; %[3.6,3.7,3.9,4.0]; %Age of interest
    aa = aa + 1;

    for timearound = 1:2 %The two different cratering models used by WernerTanaka2011
    
        if timearound == 1
            PF = 2; %Production function
            CM = 2; %Chronology model
        elseif timearound == 2
            PF = 1; %Production function
            CM = 3; %Chronology model
        end
    
    
%Ivanov checked against WT11Icarus Fig 1d, looks good
%Hartmann checked against WT11Icarus Figure 1c, looks good (cursory check)
a0corr = log10(interp1(age,ncr(CM,:),ats)) - interp1(D,ncum(PF,:),1); %Absolute number of additional craters
for l =1:length(ats)
    tsi(l,:) = ncum(PF,:) + a0corr(l);  %To-show isochrons
end 

%Contour the accumulation rate:
%NOT NECCESSARY, TESTING IS DONE

%For a range of sizes of interest
sthr = D; %Do not change this line km
ND = logspace(-10,10,401);%logspace(-10,10,41); Altered in revisions to improve Bayesian error bar precision
    
cakm = countarea./1e12;

small_dt = 0.025; %Gyr    


a0corr = log10(interp1(age,ncr(CM,:),[aoi-small_dt aoi+small_dt])) - interp1(D,ncum(PF,:),1); %Absolute number of additional craters
    n1 = ncum(PF,:) + a0corr(1);  %To-show isochrons
    n2 = ncum(PF,:) + a0corr(2);  %To-show isochrons


        for l = 1:length(sthr) %A range of minimum-crater-size threshholds

            ns(l) = sum(sizes>sthr(l))./cakm; %size threshold = sthr. Convert to 'per km2'
            cns(l) = ns(l) + sum(csizes>sthr(l))./cakm; %size threshold = sthr. Convert to 'per km2'

            soffset(l) = interp1(D,ncum(PF,:),sthr(l)) - interp1(D,ncum(PF,:),1); %offset for this size

            %Find the flux *per year* at this size range *at D>1km*

            %n2 = interp1(age,ncr(CM,:),aoi+small_dt) + a0corr + soffset(l);%log10
            %n1 = interp1(age,ncr(CM,:),aoi-small_dt) + a0corr + soffset(l);%log10

            %Diff at age of interest; convert to rate in N(>)/yr *for 1km*
            daoi(l) = ((10.^n2(l) - 10.^n1(l))./(2.*small_dt))./1e9; %per year

            %(fluxatsr is per 10^6 km^2)
            fluxatsr(l) = daoi(l);%(Rate of new craters per year at this D(>) size range)
            %within the study area
            fluxatsr(l) = fluxatsr(l) .* cakm;

            accumfor1(l) = fluxatsr(l)*obliterationdepthfraction*sthr(l);
            car(aa,l,:) = accumfor1(l)./ND;


        end

        
%*****************MAXIMUM LIKELIHOOD ESTIMATE PART OF LOOP (LEGACY SQRT(N) ERRORBARS) **********************

%Find maximum-likelihood accumulation rate for each timearound
accr = logspace(-8,0,5001);


size_range_to_fit = find(sthr>0.29 & ns > 0); %300m is about the breakpoint. Include size class 292m.

    for datatype = 1:2
        
        if datatype == 1;%confirmed only
            ctf = ns;%Count to fit
        elseif datatype == 2; %include 'candidate' RCIs as well as confirmed
            ctf = cns; %Count to fit
        end

            for tshr = size_range_to_fit;
                %Note error bars are just sqrt(N)
                sig = (sqrt(ctf(tshr).*cakm)./cakm); %sigma   
                model(datatype,timearound,aa,tshr,:) = interp1(ND,squeeze(car(aa,tshr,:)),accr);

                L(datatype,timearound,aa,tshr,:) = ...
                      exp(- (model(datatype,timearound,aa,tshr,:) - ctf(tshr)).^2 / (2*sig.^2)); %Wall & Jenkins Practical Stat. for Astronomers, 2nd edition, p.129
            end

            for taccr = 1:length(accr) %log-likelihoods
                mle(datatype,timearound,aa,taccr) = sum(log([squeeze(L(datatype,timearound,aa,size_range_to_fit,taccr))]));
            end

    end
    
    [C,I] = max(mle,[],4);
    %plot(model(datatype,timearound,aa,I(datatype,timearound,aa)),[0.8 0.8 0.8],'LineWidth',2)

%***************** ADDED IN REVISION: BAYESIAN ERRORBARS **********************

%Note: The reviewers didn't object to the sqrt(N) erorrbars in our
%submitted ms., and Bayesian errorbars are unlikely to affect our conclusions, but
%it's bad practice not to do it properly when we can.

%implied flat log prior
accr = logspace(-8,0,5001); %vector of trial accumulation rates

size_range_to_fit = find(sthr>0.29 & ns > 0); %300m is about the breakpoint. Include size class 292m.


    for datatype = 1:2
        
        if datatype == 1;%confirmed only
            ctf = ns;%Count to fit
        elseif datatype == 2; %include 'candidate' RCIs as well as confirmed
            ctf = cns; %Count to fit
        end

            for tshr = find(sthr>0.1); %size_range_to_fit
                %psuedo for each trial accum
                %pseudo what is the probability 

                k = round(ctf(tshr).*cakm);%number of events actually observed
                if abs(k-ctf(tshr).*cakm)>1e-10;error('Non-integer number of impacts?!');end
                
                if k > 0
                
                    lambda = cakm.*interp1(ND,squeeze(car(aa,tshr,:)),accr);%expected value of events given my trail accum rate
                    poissonlog = (log10(lambda.^k.* exp(-lambda)) - log10(factorial(k))); %Sivia 2nd edition equation 5.49.             
                    poisson = 10.^poissonlog;
                    poisson(isnan(poisson)==1)=0;
                    poisson = cumsum(poisson)./sum(poisson);

               model(datatype,timearound,aa,tshr,:) = poisson;
               lambdam(datatype,timearound,aa,tshr,:) = lambda;
               [BP,IP,JP] = unique(poisson);

               lowaccum =  interp1(BP,log10(accr(IP)),  0.5*0.3173); %log10(1 sigma low accumulation rate)
               highaccum = interp1(BP,log10(accr(IP)),1-0.5*0.3173); %log10(1 sigma high accumulation rate)

               %accumulation rates:
               errorlow(datatype,timearound,aa,tshr) = 10.^lowaccum;
              errorhigh(datatype,timearound,aa,tshr) = 10.^highaccum;
            bayesmedian(datatype,timearound,aa,tshr) = 10.^interp1(BP,log10(accr(IP)),         0.5);

                %convert to craters-per-km2 (craterdensity):
                cdlow(datatype,timearound,aa,tshr) = interp1(log10(ND),squeeze(car(aa,tshr,:)), lowaccum); %;lowaccum is in log form
               cdhigh(datatype,timearound,aa,tshr) = interp1(log10(ND),squeeze(car(aa,tshr,:)),highaccum);
                    end
                
            end
    end
    
%*****************DISPLAY**********************

    
if timearound == 2
    %The 0.9999999 is to allow for an error bar that doesn't touch zero
    [C,h] = contour(logspace(log10(min(sthr)),log10(max(sthr)),481), ...
                    logspace(log10(min(ND)),log10(max(ND)), 3201), ... %changed in revision from 321
                    interp2(log10(squeeze(car(aa,:,:))'),3),[-6:1:0],'Color',[0.2 0.2 0.2]);
    
    %Add the best fitting accumulation rate as a line on the plot:
    [C,h] = contour(logspace(log10(min(sthr)),log10(max(sthr)),481), ...
                    logspace(log10(min(ND)),log10(max(ND)), 3201), ... %changed in revision from 321
                    interp2(log10(squeeze(car(aa,:,:))'),3),[log10(accr(I(datatype,timearound,aa)))], ...
                    'LineWidth',2,'Color','b'); %Dark gray ([0.8 0.8 0.8]) replaced by blue on direct reviewer suggestion
            
    %These are totally lame sqrt(N) legacy errorbars; deprecated            
    %errorbar(sthr(sthr>0.1),ns(sthr>0.1),(sqrt(ns(sthr>0.1).*cakm)./cakm)*0.999999999,'r*','LineStyle','-','LineWidth',1)
    %hold on
    %errorbar(sthr(sthr>0.1),cns(sthr>0.1),(sqrt(cns(sthr>0.1).*cakm)./cakm)*0.999999999,'k*','LineStyle','-','LineWidth',1)
   
    %Here are the proper Bayesian errorbars. There is almost no difference
    %from sqrt(N) except for the very largest bins (fewest craters)
    line(sthr(sthr>0.1),ns(sthr>0.1),'Color','r','LineWidth',1);
    hold on
    line(sthr(sthr>0.1),cns(sthr>0.1),'Color','k','LineWidth',1);
    
    for tshr = find(sthr>0.1) 
        if tshr <= size(cdlow,4)
             errorbar(sthr(tshr),ns(tshr),ns(tshr)-cdlow(1,timearound,aa,tshr), ...
                 cdhigh(1,timearound,aa,tshr)-ns(tshr),'r*','LineWidth',1)
             errorbar(sthr(tshr),cns(tshr),cns(tshr)-cdlow(2,timearound,aa,tshr), ...
                 cdhigh(2,timearound,aa,tshr)-cns(tshr),'k*','LineWidth',1)
 
        end
    end
    %errorbar(sthr(sth r>0.1),ns(sthr>0.1),interp(),'b*','LineStyle','-','LineWidth',1)
    hold on
    %errorbar(sthr(sthr>0.1),cns(sthr>0.1),interp(),'m*','LineStyle','-','LineWidth',1)
    
    
    
    %end of revised section doing the proper errorbars.
    
    set(gca,'yscale','log')
    set(gca,'xscale','log')
    xlim([10^-(1.5),1e1]) %xlim([1e-2,1e1])
    ylim([1e-9,1e2]) %ylim([1e-10,1e4])
    %clabel(C,h);
    box on
    set(gca,'FontSize',20)
    xlabel('Crater Diameter (km)','FontSize',20)
    ylabel('NCraters (>Diameter)/km^2','FontSize',20)
    labelname = strcat('Assumed Age = ',strcat(' ',num2str(aoi)),' Gyr');
    text(0.35,1e1,labelname, ...
    'FontSize',18)
    for lindex = -6:0
            label_location = interp1(ND,squeeze(car(aa,end,:)),10^lindex);
            label_location = 10^(log10(label_location) +0.7); %graphical shift, for clarity
            labelname = strcat('10^{',num2str(lindex),'} m/yr');
            text(3,label_location,labelname, ...
            'FontSize',18)
    end
   % title('Aeolis-Zephyria Planum Crater-River Interaction Frequency','FontSize',16)
   % colorbar
    axis square
    set(gca,'ytick',logspace(-8,2,6))
   % legend({'Ivanov PF, Hartmann+Neukum CF';'Hartmann PF, Hartmann05 CF';'Crater-River Interactions';'Crater-River Interactions including Candidates';'Best fit'},'FontSize',12)
elseif timearound == 1
    figure
    [C,h] = contour(logspace(log10(min(sthr)),log10(max(sthr)),481), ...
                    logspace(log10(min(ND)),log10(max(ND)), 3201), ... %changed in revision to '3201' from '321'
                    interp2(log10(squeeze(car(aa,:,:))'),3),[-6:1:0], ...
                    'LineStyle','-.','LineWidth',1,'Color','m');%'Color',[0.2 0.2 0.2]);
                hold on
end

end
    
end


%*****************SUMMARY GRAPHICS**********************
figure
plot(accr,squeeze(mle(1,1,:,:))','LineStyle',':')
hold on
plot(accr,squeeze(mle(2,1,:,:))')
plot(accr,squeeze(mle(1,2,:,:))','LineStyle',':','LineWidth',2)
hold on
plot(accr,squeeze(mle(2,2,:,:))','LineWidth',2)
set(gca,'xscale','log'); set(gca,'yscale','log')
grid on
xlabel('Accumulation rate (m/yr)')
ylabel('Natural Log of Likelihood (Yes, Two Logs)')


%plot(age(1:end-1),dncr(:,1:end)./min(dncr'))

save Bayes_sediment_output_REVISIONS
