diff --git a/.gitignore b/.gitignore
index f7f3d9d..7f7e849 100644
--- a/.gitignore
+++ b/.gitignore
@@ -18,6 +18,8 @@ Manifest.toml
 examples/DenverUPbExampleData/results.csv
 results.csv
 distresults.csv
+AgeDepthModel.csv
+stratonly_example_height_results.csv
 
 .ipynb_check*
 /etc
diff --git a/examples/Chron1.0Coupled.jl b/examples/Chron1.0Coupled.jl
index 7aad51c..c2b3a6b 100644
--- a/examples/Chron1.0Coupled.jl
+++ b/examples/Chron1.0Coupled.jl
@@ -111,7 +111,7 @@
 
     # Run the stratigraphic MCMC model
     @time (mdl, agedist, lldist) = StratMetropolisDist(smpl, config)
-
+    exportdataset(NamedTuple(mdl), "AgeDepthModel.csv")
 
 ## --- Plot stratigraphic model - - - - - - - - - - - - - - - - - - - - - - - -
 
@@ -178,7 +178,13 @@
     dhdt_84p = nanpctile(dhdt_dist,84.135,dim=2) # Plus 1-sigma (84.135th percentile)
 
     # Plot results
-    hdl = plot(agebincenters,dhdt, label="Mean", color=:black, linewidth=2)
+    hdl = plot(
+        xlabel="Age ($(smpl.Age_Unit))", 
+        ylabel="Depositional Rate ($(smpl.Height_Unit) / $(smpl.Age_Unit) over $binwidth $(smpl.Age_Unit))", 
+        fg_color_legend=:white,
+        framestyle=:box,
+    )
+    plot!(hdl, agebincenters,dhdt, label="Mean", color=:black, linewidth=2)
     plot!(hdl,[agebincenters; reverse(agebincenters)],[dhdt_16p; reverse(dhdt_84p)], fill=(0,0.2,:darkblue), linealpha=0, label="68% CI")
     for lci in 20:5:45
         dhdt_lp = nanpctile(dhdt_dist,lci,dim=2)
@@ -186,8 +192,7 @@
         plot!(hdl,[agebincenters; reverse(agebincenters)],[dhdt_lp; reverse(dhdt_up)], fill=(0,0.2,:darkblue), linealpha=0, label="")
     end
     plot!(hdl, agebincenters,dhdt_50p, label="Median", color=:grey, linewidth=1)
-    plot!(hdl, xlabel="Age ($(smpl.Age_Unit))", ylabel="Depositional Rate ($(smpl.Height_Unit) / $(smpl.Age_Unit) over $binwidth $(smpl.Age_Unit))", fg_color_legend=:white)
-    # savefig(hdl,"DepositionRateModelCI.pdf")
+    savefig(hdl,"DepositionRateModelCI.pdf")
     display(hdl)
 
 ## --- Make heatmap - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
@@ -250,12 +255,15 @@
 
     # Run the model. Note the additional `hiatus` arguments
     @time (mdl, agedist, hiatusdist, lldist) = StratMetropolisDist(smpl, hiatus, config)
+    exportdataset(NamedTuple(mdl), "AgeDepthModel.csv")
 
     # Plot results (mean and 95% confidence interval for both model and data)
     hdl = plot([mdl.Age_025CI; reverse(mdl.Age_975CI)],[mdl.Height; reverse(mdl.Height)], fill=(minimum(mdl.Height),0.5,:blue), label="model")
     plot!(hdl, mdl.Age, mdl.Height, linecolor=:blue, label="")
     plot!(hdl, smpl.Age, smpl.Height, xerror=(smpl.Age-smpl.Age_025CI,smpl.Age_975CI-smpl.Age),label="data",seriestype=:scatter,color=:black)
     plot!(hdl, xlabel="Age (Ma)", ylabel="Height (cm)", framestyle=:box)
+    savefig(hdl,"AgeDepthModel.pdf")
+    display(hdl)
 
 ## --- (Optional) Add systematic uncertainties for U-Pb data
 
diff --git a/examples/Chron1.0CoupledConcordia.jl b/examples/Chron1.0CoupledConcordia.jl
index 9c7619a..6355725 100644
--- a/examples/Chron1.0CoupledConcordia.jl
+++ b/examples/Chron1.0CoupledConcordia.jl
@@ -110,7 +110,7 @@
 
     # Run the stratigraphic MCMC model
     @time (mdl, agedist, lldist) = StratMetropolisDist(smpl, config)
-
+    exportdataset(NamedTuple(mdl), "AgeDepthModel.csv")
 
 ## --- Plot stratigraphic model - - - - - - - - - - - - - - - - - - - - - - - -
 
@@ -178,7 +178,13 @@
     dhdt_84p = nanpctile(dhdt_dist,84.135,dim=2) # Plus 1-sigma (84.135th percentile)
 
     # Plot results
-    hdl = plot(agebincenters,dhdt, label="Mean", color=:black, linewidth=2)
+    hdl = plot(
+        xlabel="Age ($(smpl.Age_Unit))", 
+        ylabel="Depositional Rate ($(smpl.Height_Unit) / $(smpl.Age_Unit) over $binwidth $(smpl.Age_Unit))", 
+        fg_color_legend=:white,
+        framestyle=:box,
+    )
+    plot!(hdl, agebincenters,dhdt, label="Mean", color=:black, linewidth=2)
     plot!(hdl,[agebincenters; reverse(agebincenters)],[dhdt_16p; reverse(dhdt_84p)], fill=(0,0.2,:darkblue), linealpha=0, label="68% CI")
     for lci in 20:5:45
         dhdt_lp = nanpctile(dhdt_dist,lci,dim=2)
@@ -186,8 +192,7 @@
         plot!(hdl,[agebincenters; reverse(agebincenters)],[dhdt_lp; reverse(dhdt_up)], fill=(0,0.2,:darkblue), linealpha=0, label="")
     end
     plot!(hdl, agebincenters,dhdt_50p, label="Median", color=:grey, linewidth=1)
-    plot!(hdl, xlabel="Age ($(smpl.Age_Unit))", ylabel="Depositional Rate ($(smpl.Height_Unit) / $(smpl.Age_Unit) over $binwidth $(smpl.Age_Unit))", fg_color_legend=:white)
-    # savefig(hdl,"DepositionRateModelCI.pdf")
+    savefig(hdl,"DepositionRateModelCI.pdf")
     display(hdl)
 
 ## --- End of File
diff --git a/examples/Chron1.0CoupledSystematic.jl b/examples/Chron1.0CoupledSystematic.jl
index 4e875e3..86e85a3 100644
--- a/examples/Chron1.0CoupledSystematic.jl
+++ b/examples/Chron1.0CoupledSystematic.jl
@@ -124,6 +124,7 @@
 
     # Run the stratigraphic MCMC model
     @time (mdl, agedist, lldist) = StratMetropolisDist(smpl, config, systematic)
+    exportdataset(NamedTuple(mdl), "AgeDepthModel.csv")
 
     # # Other options:
     # # Youngest Zircon
@@ -213,7 +214,13 @@
     dhdt_84p = nanpctile(dhdt_dist,84.135,dim=2) # Plus 1-sigma (84.135th percentile)
 
     # Plot results
-    hdl = plot(agebincenters,dhdt, label="Mean", color=:black, linewidth=2)
+    hdl = plot(
+        xlabel="Age ($(smpl.Age_Unit))", 
+        ylabel="Depositional Rate ($(smpl.Height_Unit) / $(smpl.Age_Unit) over $binwidth $(smpl.Age_Unit))", 
+        fg_color_legend=:white,
+        framestyle=:box,
+    )
+    plot!(hdl, agebincenters,dhdt, label="Mean", color=:black, linewidth=2)
     plot!(hdl,[agebincenters; reverse(agebincenters)],[dhdt_16p; reverse(dhdt_84p)], fill=(0,0.2,:darkblue), linealpha=0, label="68% CI")
     for lci in 20:5:45
         dhdt_lp = nanpctile(dhdt_dist,lci,dim=2)
@@ -221,8 +228,7 @@
         plot!(hdl,[agebincenters; reverse(agebincenters)],[dhdt_lp; reverse(dhdt_up)], fill=(0,0.2,:darkblue), linealpha=0, label="")
     end
     plot!(hdl, agebincenters,dhdt_50p, label="Median", color=:grey, linewidth=1)
-    plot!(hdl, xlabel="Age ($(smpl.Age_Unit))", ylabel="Depositional Rate ($(smpl.Height_Unit) / $(smpl.Age_Unit) over $binwidth $(smpl.Age_Unit))", fg_color_legend=:white)
-    # savefig(hdl,"DepositionRateModelCI.pdf")
+    savefig(hdl,"DepositionRateModelCI.pdf")
     display(hdl)
 
 ## --- Probability that a given interval of stratigraphy was deposited entirely before/after a given time
diff --git a/examples/Chron1.0CoupledTabularInput.jl b/examples/Chron1.0CoupledTabularInput.jl
index 2418798..e3b09da 100644
--- a/examples/Chron1.0CoupledTabularInput.jl
+++ b/examples/Chron1.0CoupledTabularInput.jl
@@ -12,18 +12,19 @@
 
     using Chron
     using Plots
+    cd(@__DIR__) # Move to script directory
 
 ## --- Define sample properties - - - - - - - - - - - - - - - - - - - - - - - -
 
     # As in Chron1.0Coupled, but here we read in the information from a csv table:
-    ds = importdataset(joinpath(@__DIR__, "coupled_example_input.csv"), importas=:Tuple)
+    ds = importdataset("coupled_example_input.csv", importas=:Tuple)
 
     smpl = ChronAgeData(length(ds.Name))
     smpl.Name      =  Tuple(ds.Name)
     smpl.Height   .=  ds.Height
     smpl.Height_sigma .= ds.Height_sigma
     smpl.Age_Sidedness .= ds.Age_Sidedness              # Sidedness (zeros by default: geochron constraints are two-sided). Use -1 for a maximum age and +1 for a minimum age, 0 for two-sided
-    smpl.Path = joinpath(@__DIR__, first(ds.Path))      # Where are the data files?
+    smpl.Path = first(ds.Path)                          # Where are the data files?
     smpl.inputSigmaLevel = first(ds.inputSigmaLevel)    # i.e., are the data files 1-sigma or 2-sigma. Integer.
     smpl.Age_Unit = first(ds.Age_Unit)                  # Unit of measurement for ages and errors in the data files
     smpl.Height_Unit = first(ds.Height_Unit)            # Unit of measurement for Height and Height_sigma
@@ -103,7 +104,7 @@
 
     # Run the stratigraphic MCMC model
     @time (mdl, agedist, lldist) = StratMetropolisDist(smpl, config)
-
+    exportdataset(NamedTuple(mdl), "AgeDepthModel.csv")
 
 ## --- Plot stratigraphic model - - - - - - - - - - - - - - - - - - - - - - - -
 
@@ -126,25 +127,19 @@
     display(hdl)
 
 
-## --- Interpolate model age at a specific stratigraphic height - - - - - - - -
-
-    # Stratigraphic height at which to interpolate
-    interp_height = 0
+## --- Interpolate results at a specific heights read from file  - - - - - - - -
 
-    age_at_height = linterp1s(mdl.Height,mdl.Age,interp_height)
-    age_at_height_min = linterp1s(mdl.Height,mdl.Age_025CI,interp_height)
-    age_at_height_max = linterp1s(mdl.Height,mdl.Age_975CI,interp_height)
-    print("Interpolated age at height=$interp_height: $age_at_height +$(age_at_height_max-age_at_height)/-$(age_at_height-age_at_height_min) $(smpl.Age_Unit)")
+    # Stratigraphic heights at which to interpolate
+    hds = importdataset("coupled_example_height_query.csv", importas=:Tuple)
 
-    # Optional: interpolate full age distribution
-    interpolated_distribution = Array{Float64}(undef,size(agedist,2))
-    for i=1:size(agedist,2)
-        interpolated_distribution[i] = linterp1s(mdl.Height,agedist[:,i],interp_height)
-    end
-    hdl = histogram(interpolated_distribution, nbins=50, label="", framestyle=:box)
-    plot!(hdl, xlabel="Age ($(smpl.Age_Unit)) at height=$interp_height", ylabel="Likelihood (unnormalized)")
-    savefig(hdl, "Interpolated age distribution.pdf")
-    display(hdl)
+    hds = (;
+        Height = hds.Height,
+        Age = linterp1s(mdl.Height, mdl.Age, hds.Height),
+        Age_sigma = linterp1s(mdl.Height, mdl.Age_sigma, hds.Height),
+        Age_025CI = linterp1s(mdl.Height, mdl.Age_025CI, hds.Height),
+        Age_975CI = linterp1s(mdl.Height, mdl.Age_975CI, hds.Height),
+    )
+    exportdataset(hds, "coupled_example_height_results.csv")
 
 ## --- Calculate deposition rate binned by age  - - - - - - - - - - - - - - - -
 
@@ -170,7 +165,13 @@
     dhdt_84p = nanpctile(dhdt_dist,84.135,dim=2) # Plus 1-sigma (84.135th percentile)
 
     # Plot results
-    hdl = plot(agebincenters,dhdt, label="Mean", color=:black, linewidth=2)
+    hdl = plot(
+        xlabel="Age ($(smpl.Age_Unit))", 
+        ylabel="Depositional Rate ($(smpl.Height_Unit) / $(smpl.Age_Unit) over $binwidth $(smpl.Age_Unit))", 
+        fg_color_legend=:white,
+        framestyle=:box,
+    )
+    plot!(hdl, agebincenters,dhdt, label="Mean", color=:black, linewidth=2)
     plot!(hdl,[agebincenters; reverse(agebincenters)],[dhdt_16p; reverse(dhdt_84p)], fill=(0,0.2,:darkblue), linealpha=0, label="68% CI")
     for lci in 20:5:45
         dhdt_lp = nanpctile(dhdt_dist,lci,dim=2)
@@ -178,8 +179,7 @@
         plot!(hdl,[agebincenters; reverse(agebincenters)],[dhdt_lp; reverse(dhdt_up)], fill=(0,0.2,:darkblue), linealpha=0, label="")
     end
     plot!(hdl, agebincenters,dhdt_50p, label="Median", color=:grey, linewidth=1)
-    plot!(hdl, xlabel="Age ($(smpl.Age_Unit))", ylabel="Depositional Rate ($(smpl.Height_Unit) / $(smpl.Age_Unit) over $binwidth $(smpl.Age_Unit))", fg_color_legend=:white)
-    # savefig(hdl,"DepositionRateModelCI.pdf")
+    savefig(hdl,"DepositionRateModelCI.pdf")
     display(hdl)
 
 ## --- Make heatmap - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
@@ -216,12 +216,15 @@
 
     # Run the model. Note the additional `hiatus` arguments
     @time (mdl, agedist, hiatusdist, lldist) = StratMetropolisDist(smpl, hiatus, config)
+    exportdataset(NamedTuple(mdl), "AgeDepthModel.csv")
 
     # Plot results (mean and 95% confidence interval for both model and data)
     hdl = plot([mdl.Age_025CI; reverse(mdl.Age_975CI)],[mdl.Height; reverse(mdl.Height)], fill=(minimum(mdl.Height),0.5,:blue), label="model")
     plot!(hdl, mdl.Age, mdl.Height, linecolor=:blue, label="")
     plot!(hdl, smpl.Age, smpl.Height, xerror=(smpl.Age-smpl.Age_025CI,smpl.Age_975CI-smpl.Age),label="data",seriestype=:scatter,color=:black)
     plot!(hdl, xlabel="Age (Ma)", ylabel="Height (cm)", framestyle=:box)
+    savefig(hdl,"AgeDepthModel.pdf")
+    display(hdl)
 
 ## --- (Optional) Add systematic uncertainties for U-Pb data
 
diff --git a/examples/Chron1.0Radiocarbon.jl b/examples/Chron1.0Radiocarbon.jl
index a0fb10c..9bde75b 100644
--- a/examples/Chron1.0Radiocarbon.jl
+++ b/examples/Chron1.0Radiocarbon.jl
@@ -84,6 +84,7 @@
 
     # Run the stratigraphic MCMC model
     @time (mdl, agedist, lldist) = StratMetropolis14C(smpl, config)
+    exportdataset(NamedTuple(mdl), "AgeDepthModel.csv")
 
     # Plot results (mean and 95% confidence interval for both model and data)
     hdl = plot(framestyle=:box,
@@ -140,7 +141,13 @@
     dhdt_84p = nanpctile(dhdt_dist,84.135,dim=2) # Plus 1-sigma (84.135th percentile)
 
     # Plot results
-    hdl = plot(agebincenters,dhdt, label="Mean", color=:black, linewidth=2)
+    hdl = plot(
+        xlabel="Age ($(smpl.Age_Unit))", 
+        ylabel="Depositional Rate ($(smpl.Height_Unit) / $(smpl.Age_Unit) over $binwidth $(smpl.Age_Unit))", 
+        fg_color_legend=:white,
+        framestyle=:box,
+    )
+    plot!(hdl, agebincenters,dhdt, label="Mean", color=:black, linewidth=2)
     plot!(hdl,[agebincenters; reverse(agebincenters)],[dhdt_16p; reverse(dhdt_84p)], fill=(0,0.2,:darkblue), linealpha=0, label="68% CI")
     for lci in 20:5:45
         dhdt_lp = nanpctile(dhdt_dist,lci,dim=2)
@@ -148,8 +155,7 @@
         plot!(hdl,[agebincenters; reverse(agebincenters)],[dhdt_lp; reverse(dhdt_up)], fill=(0,0.2,:darkblue), linealpha=0, label="")
     end
     plot!(hdl, agebincenters,dhdt_50p, label="Median", color=:grey, linewidth=1)
-    plot!(hdl, xlabel="Age ($(smpl.Age_Unit))", ylabel="Depositional Rate ($(smpl.Height_Unit) / $(smpl.Age_Unit) over $binwidth $(smpl.Age_Unit))", fg_color_legend=:white)
-    # savefig(hdl,"DepositionRateModelCI.pdf")
+    savefig(hdl,"DepositionRateModelCI.pdf")
     display(hdl)
 
 ## --- Optional: Stratigraphic model including hiatuses - - - - - - - - - - - -
@@ -164,6 +170,7 @@
 
     # Run the model. Note the additional `hiatus` arguments
     @time (mdl, agedist, hiatusdist, lldist) = StratMetropolis14C(smpl, hiatus, config)
+    exportdataset(NamedTuple(mdl), "AgeDepthModel.csv")
 
     # Plot results (mean and 95% confidence interval for both model and data)
     hdl = plot(framestyle=:box,
@@ -173,5 +180,7 @@
     plot!(hdl, [mdl.Age_025CI; reverse(mdl.Age_975CI)],[mdl.Height; reverse(mdl.Height)], fill=(minimum(mdl.Height),0.5,:blue), label="model")
     plot!(hdl, mdl.Age, mdl.Height, linecolor=:blue, label="", fg_color_legend=:white)
     plot!(hdl, smpl.Age, smpl.Height, xerror=(smpl.Age-smpl.Age_025CI,smpl.Age_975CI-smpl.Age),label="data",seriestype=:scatter,color=:black)
+    savefig(hdl, "Interpolated age distribution.pdf")
+    display(hdl)
 
 ## --- End of File - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
diff --git a/examples/Chron1.0StratOnly.jl b/examples/Chron1.0StratOnly.jl
index 56c24da..bf71a77 100644
--- a/examples/Chron1.0StratOnly.jl
+++ b/examples/Chron1.0StratOnly.jl
@@ -50,6 +50,7 @@
 
     # Run the stratigraphic MCMC model
     @time (mdl, agedist, lldist) = StratMetropolis(smpl, config)
+    exportdataset(NamedTuple(mdl), "AgeDepthModel.csv")
 
     # Plot results (mean and 95% CI for model / 2-sigma for data)
     hdl = plot(framestyle=:box,
@@ -114,7 +115,13 @@
     dhdt_84p = nanpctile(dhdt_dist,84.135,dim=2) # Plus 1-sigma (84.135th percentile)
 
     # Plot results
-    hdl = plot(agebincenters,dhdt, label="Mean", color=:black, linewidth=2)
+    hdl = plot(
+        xlabel="Age ($(smpl.Age_Unit))", 
+        ylabel="Depositional Rate ($(smpl.Height_Unit) / $(smpl.Age_Unit) over $binwidth $(smpl.Age_Unit))", 
+        fg_color_legend=:white,
+        framestyle=:box,
+    )
+    plot!(hdl, agebincenters,dhdt, label="Mean", color=:black, linewidth=2)
     plot!(hdl,[agebincenters; reverse(agebincenters)],[dhdt_16p; reverse(dhdt_84p)], fill=(0,0.2,:darkblue), linealpha=0, label="68% CI")
     for lci in 20:5:45
         dhdt_lp = nanpctile(dhdt_dist,lci,dim=2)
@@ -122,10 +129,8 @@
         plot!(hdl,[agebincenters; reverse(agebincenters)],[dhdt_lp; reverse(dhdt_up)], fill=(0,0.2,:darkblue), linealpha=0, label="")
     end
     plot!(hdl, agebincenters,dhdt_50p, label="Median", color=:grey, linewidth=1)
-    plot!(hdl, xlabel="Age ($(smpl.Age_Unit))", ylabel="Depositional Rate ($(smpl.Height_Unit) / $(smpl.Age_Unit) over $binwidth $(smpl.Age_Unit))", fg_color_legend=:white)
-    # savefig(hdl,"DepositionRateModelCI.pdf")
+    savefig(hdl,"DepositionRateModelCI.pdf")
     display(hdl)
-
 ## --- Optional: Stratigraphic model including hiatuses - - - - - - - - - - - -
 
     # Data about hiatuses
@@ -138,11 +143,14 @@
 
     # Run the model. Note the additional `hiatus` arguments
     @time (mdl, agedist, hiatusdist, lldist) = StratMetropolis(smpl, hiatus, config); sleep(0.5)
+    exportdataset(NamedTuple(mdl), "AgeDepthModel.csv")
 
     # Plot results (mean and 95% confidence interval for both model and data)
     hdl = plot([mdl.Age_025CI; reverse(mdl.Age_975CI)],[mdl.Height; reverse(mdl.Height)], fill=(minimum(mdl.Height),0.5,:blue), label="model")
     plot!(hdl, mdl.Age, mdl.Height, linecolor=:blue, label="", fg_color_legend=:white)
     plot!(hdl, smpl.Age, smpl.Height, xerror=smpl.Age_sigma*2,label="data",seriestype=:scatter,color=:black)
     plot!(hdl, xlabel="Age ($(smpl.Age_Unit))", ylabel="Height ($(smpl.Height_Unit))")
+    savefig(hdl, "Interpolated age distribution.pdf")
+    display(hdl)
 
 ## --- End of File - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
diff --git a/examples/Chron1.0StratOnlyTabularInput.jl b/examples/Chron1.0StratOnlyTabularInput.jl
index ea09e61..b8d61c3 100644
--- a/examples/Chron1.0StratOnlyTabularInput.jl
+++ b/examples/Chron1.0StratOnlyTabularInput.jl
@@ -9,13 +9,14 @@
 
     using Chron
     using Plots
+    cd(@__DIR__) # Move to script directory
 
 ## --- Define sample properties - - - - - - - - - - - - - - - - - - - - - - - -
 
     # As in Chron1.0Coupled, but here we read in the information from a csv table:
 
     # # # # # # # # # # # Enter sample information here # # # # # # # # # # # # #
-    ds = importdataset(joinpath(@__DIR__, "stratonly_example_input.csv"), importas=:Tuple)
+    ds = importdataset("stratonly_example_input.csv", importas=:Tuple)
     smpl = ChronAgeData(length(ds.Name))
     smpl.Name       = Tuple(ds.Name)
     smpl.Age       .= ds.Age
@@ -50,6 +51,7 @@
 
     # Run the stratigraphic MCMC model
     @time (mdl, agedist, lldist) = StratMetropolis(smpl, config)
+    exportdataset(NamedTuple(mdl), "AgeDepthModel.csv")
 
     # Plot results (mean and 95% CI for model / 2-sigma for data)
     hdl = plot(framestyle=:box,
@@ -70,25 +72,19 @@
     savefig(hdl,"AgeDepthModel.pdf")
     display(hdl)
 
-## --- Interpolate results at a specific height - - - - - - - - - - - - - - - -
+## --- Interpolate results at a specific heights read from file  - - - - - - - -
 
-    # Stratigraphic height at which to interpolate
-    height = -404
+    # Stratigraphic heights at which to interpolate
+    hds = importdataset("stratonly_example_height_query.csv", importas=:Tuple)
 
-    age_at_height = linterp1s(mdl.Height,mdl.Age,height)
-    age_at_height_min = linterp1s(mdl.Height,mdl.Age_025CI,height)
-    age_at_height_max = linterp1s(mdl.Height,mdl.Age_975CI,height)
-    print("Interpolated age at height=$height: $age_at_height +$(age_at_height_max-age_at_height)/-$(age_at_height-age_at_height_min) $(smpl.Age_Unit)")
-
-    # Optional: interpolate full age distribution
-    interpolated_distribution = Array{Float64}(undef,size(agedist,2))
-    for i=1:size(agedist,2)
-        interpolated_distribution[i] = linterp1s(mdl.Height,agedist[:,i],height)
-    end
-    hdl = histogram(interpolated_distribution, nbins=50, label="")
-    plot!(hdl, xlabel="Age ($(smpl.Age_Unit)) at height=$height", ylabel="Likelihood (unnormalized)")
-    savefig(hdl, "Interpolated age distribution.pdf")
-    display(hdl)
+    hds = (;
+        Height = hds.Height,
+        Age = linterp1s(mdl.Height, mdl.Age, hds.Height),
+        Age_sigma = linterp1s(mdl.Height, mdl.Age_sigma, hds.Height),
+        Age_025CI = linterp1s(mdl.Height, mdl.Age_025CI, hds.Height),
+        Age_975CI = linterp1s(mdl.Height, mdl.Age_975CI, hds.Height),
+    )
+    exportdataset(hds, "stratonly_example_height_results.csv")
 
 ## --- Calculate deposition rate binned by age - - - - - - - - - - - - - - - - -
 
@@ -114,7 +110,13 @@
     dhdt_84p = nanpctile(dhdt_dist,84.135,dim=2) # Plus 1-sigma (84.135th percentile)
 
     # Plot results
-    hdl = plot(agebincenters,dhdt, label="Mean", color=:black, linewidth=2)
+    hdl = plot(
+        xlabel="Age ($(smpl.Age_Unit))", 
+        ylabel="Depositional Rate ($(smpl.Height_Unit) / $(smpl.Age_Unit) over $binwidth $(smpl.Age_Unit))", 
+        fg_color_legend=:white,
+        framestyle=:box,
+    )
+    plot!(hdl, agebincenters,dhdt, label="Mean", color=:black, linewidth=2)
     plot!(hdl,[agebincenters; reverse(agebincenters)],[dhdt_16p; reverse(dhdt_84p)], fill=(0,0.2,:darkblue), linealpha=0, label="68% CI")
     for lci in 20:5:45
         dhdt_lp = nanpctile(dhdt_dist,lci,dim=2)
@@ -122,8 +124,7 @@
         plot!(hdl,[agebincenters; reverse(agebincenters)],[dhdt_lp; reverse(dhdt_up)], fill=(0,0.2,:darkblue), linealpha=0, label="")
     end
     plot!(hdl, agebincenters,dhdt_50p, label="Median", color=:grey, linewidth=1)
-    plot!(hdl, xlabel="Age ($(smpl.Age_Unit))", ylabel="Depositional Rate ($(smpl.Height_Unit) / $(smpl.Age_Unit) over $binwidth $(smpl.Age_Unit))", fg_color_legend=:white)
-    # savefig(hdl,"DepositionRateModelCI.pdf")
+    savefig(hdl,"DepositionRateModelCI.pdf")
     display(hdl)
 
 ## --- Optional: Stratigraphic model including hiatuses - - - - - - - - - - - -
@@ -138,11 +139,14 @@
 
     # Run the model. Note the additional `hiatus` arguments
     @time (mdl, agedist, hiatusdist, lldist) = StratMetropolis(smpl, hiatus, config); sleep(0.5)
+    exportdataset(NamedTuple(mdl), "AgeDepthModel.csv")
 
     # Plot results (mean and 95% confidence interval for both model and data)
     hdl = plot([mdl.Age_025CI; reverse(mdl.Age_975CI)],[mdl.Height; reverse(mdl.Height)], fill=(minimum(mdl.Height),0.5,:blue), label="model")
     plot!(hdl, mdl.Age, mdl.Height, linecolor=:blue, label="", fg_color_legend=:white)
     plot!(hdl, smpl.Age, smpl.Height, xerror=smpl.Age_sigma*2,label="data",seriestype=:scatter,color=:black)
     plot!(hdl, xlabel="Age ($(smpl.Age_Unit))", ylabel="Height ($(smpl.Height_Unit))")
+    savefig(hdl,"AgeDepthModel.pdf")
+    display(hdl)
 
 ## --- End of File - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
diff --git a/examples/coupled_example_height_query.csv b/examples/coupled_example_height_query.csv
new file mode 100644
index 0000000..75405d0
--- /dev/null
+++ b/examples/coupled_example_height_query.csv
@@ -0,0 +1,147 @@
+Height
+-52
+-51
+-50
+-49
+-48
+-47
+-46
+-45
+-44
+-43
+-42
+-41
+-40
+-39
+-38
+-37
+-36
+-35
+-34
+-33
+-32
+-31
+-30
+-29
+-28
+-27
+-26
+-25
+-24
+-23
+-22
+-21
+-20
+-19
+-18
+-17
+-16
+-15
+-14
+-13
+-12
+-11
+-10
+-9
+-8
+-7
+-6
+-5
+-4
+-3
+-2
+-1
+0
+1
+2
+3
+4
+5
+6
+7
+8
+9
+10
+11
+12
+13
+14
+15
+16
+17
+18
+19
+20
+21
+22
+23
+24
+25
+26
+27
+28
+29
+30
+31
+32
+33
+34
+35
+36
+37
+38
+39
+40
+41
+42
+43
+44
+45
+46
+47
+48
+49
+50
+51
+52
+53
+54
+55
+56
+57
+58
+59
+60
+61
+62
+63
+64
+65
+66
+67
+68
+69
+70
+71
+72
+73
+74
+75
+76
+77
+78
+79
+80
+81
+82
+83
+84
+85
+86
+87
+88
+89
+90
+91
+92
+93
diff --git a/examples/stratonly_example_height_query.csv b/examples/stratonly_example_height_query.csv
new file mode 100644
index 0000000..7ad7c5f
--- /dev/null
+++ b/examples/stratonly_example_height_query.csv
@@ -0,0 +1,58 @@
+Height
+-355
+-356
+-357
+-358
+-359
+-360
+-361
+-362
+-363
+-364
+-365
+-366
+-367
+-368
+-369
+-370
+-371
+-372
+-373
+-374
+-375
+-376
+-377
+-378
+-379
+-380
+-381
+-382
+-383
+-384
+-385
+-386
+-387
+-388
+-389
+-390
+-391
+-392
+-393
+-394
+-395
+-396
+-397
+-398
+-399
+-400
+-401
+-402
+-403
+-404
+-405
+-406
+-407
+-408
+-409
+-410
+-411