---
title: "Titan Interior Evolution Visualization Toolkit"
subtitle: "Interactive analysis of thermodynamic models for Titan's interior"
author: "M. Melwani Daswani (Jet Propulsion Laboratory, California Institute of Technology)"
date: "May 31, 2025"
output: 
  html_notebook:
    toc: true
    toc_float: true
    code_folding: show
    theme: flatly
    highlight: tango
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE, fig.width = 10, fig.height = 6)
```

# Overview

This notebook was constructed for the paper by K. E. Miller*, M. Melwani Daswani, C. Sotin, C. S. Cockell, C. Neish, M. J. Malaska, K. K. Farnsworth, C. Nixon, K. M. Soderlund, P. M. Higgins, A. Affholder, K. Kalousová, and R. M. C. Lopes "Titan’s refractory core evolution: implications for organics in its subsurface ocean".

For correspondence about the paper, please direct enquiries to K. E. Miller at kelly.miller@swri.org.
For enquiries about this visualization toolkit R notebook, please direct enquiries to M. Melwani Daswani at daswani@jpl.nasa.gov.

This notebook provides interactive visualization tools for analyzing results from thermodynamic models of Titan's interior evolution. The models have already been pre-run, and the data are presented here. Please include the uncompressed "titan_data" folder in the same folder as this notebook, or change the paths in the code below to direct to the uncompressed "titan_data" folder. The analysis covers three compositional cases:

- **Case A**: Chondritic composition (100% CI chondrite)
- **Case B**: Mixed composition (26% chondritic, 74% cometary) 
- **Case C**: Cometary composition (100% cometary)

## What's Included

- **Time Series Analysis**: Element extraction, density evolution, hydrosphere growth
- **Phase Assemblage Diagrams**: Mineral stability at different pressures and temperatures
- **Density Evolution**: Interior structure changes over time
- **Fluid Chemistry**: Aqueous solute concentrations (Cases A & C)
- **Mass Fractions**: Composition of exsolved species
- **Cross-Case Comparisons**: Side-by-side analysis of all three cases

# Setup and Data Loading

## Load Required Libraries

```{r load_libraries}
# Core data manipulation and visualization
library(tidyverse)
library(viridis)
library(ggpubr)
library(scales)
library(patchwork)

# Set default theme for all plots
theme_set(theme_bw())

# Define consistent color scheme for cases
case_colors <- c("Case A" = "#648FFF", "Case B" = "#DC267F", "Case C" = "#FFB000")
case_shapes <- c("Case A" = 16, "Case B" = 17, "Case C" = 18)

# Physical constants
TITAN_HYDROSPHERE_MASS <- 3.74e22  # kg (reference)

cat("✓ Libraries loaded successfully\n")
cat("✓ Visualization theme and color scheme set\n")
```

## Data Loading Functions

```{r data_loading_functions}
#' Load and validate CSV data with error handling
#' 
#' @param file_path Path to CSV file
#' @param data_type Description of data type for error messages
#' @return Data frame or NULL if loading fails
load_data_safe <- function(file_path, data_type = "data") {
  if (!file.exists(file_path)) {
    cat(sprintf("⚠️  Warning: %s not found at %s\n", data_type, file_path))
    return(NULL)
  }
  
  tryCatch({
    data <- read.csv(file_path, stringsAsFactors = FALSE)
    cat(sprintf("✓ Loaded %s (%d rows, %d columns)\n", data_type, nrow(data), ncol(data)))
    return(data)
  }, error = function(e) {
    cat(sprintf("❌ Error loading %s: %s\n", data_type, e$message))
    return(NULL)
  })
}

#' Load time series data for a specific case and element
#' 
#' @param case_name Case identifier ("case_a", "case_b", "case_c")
#' @param element Element name ("c", "h", "total_mass", etc.)
#' @param data_type "extracted", "retained", or "density_evolution"
#' @return Data frame with time series data
load_timeseries <- function(case_name, element, data_type = "extracted") {
  if (data_type == "density_evolution") {
    file_path <- sprintf("titan_data/timeseries/%s_density_evolution.csv", case_name)
  } else {
    file_path <- sprintf("titan_data/timeseries/%s_%s_%s.csv", case_name, element, data_type)
  }
  
  data_description <- sprintf("%s %s %s", case_name, element, data_type)
  return(load_data_safe(file_path, data_description))
}

#' Load phase assemblage data
#' 
#' @param case_name Case identifier
#' @param pressure Pressure level ("1GPa", "2500MPa", "3240MPa", "5GPa")
#' @return Data frame with phase assemblage data
load_phases <- function(case_name, pressure) {
  file_path <- sprintf("titan_data/phase_assemblages/%s_phases_%s.csv", case_name, pressure)
  data_description <- sprintf("%s phases at %s", case_name, pressure)
  return(load_data_safe(file_path, data_description))
}

cat("✓ Data loading functions defined\n")
```

## Load Essential Datasets

```{r load_essential_data}
cat("Loading essential datasets...\n")

# Hydrosphere comparison data
hydrosphere_comparison <- load_data_safe("titan_data/timeseries/hydrosphere_comparison_all_cases.csv", 
                                        "Hydrosphere comparison")

# Thermal evolution data
thermal_evolution <- load_data_safe("titan_data/thermal_evolution/sotin_thermal_evolution.csv", 
                                   "Sotin thermal evolution")

# Fluid chemistry data (Cases A & C only)
case_a_fluid <- load_data_safe("titan_data/fluid_chemistry/case_a_aqueous_solutes.csv", 
                               "Case A fluid chemistry")
case_c_fluid <- load_data_safe("titan_data/fluid_chemistry/case_c_aqueous_solutes.csv", 
                               "Case C fluid chemistry")

# Mass fractions (Cases A & C only)
case_a_mass_fractions <- load_data_safe("titan_data/mass_fractions/case_a_exsolved_species.csv", 
                                       "Case A mass fractions")
case_c_mass_fractions <- load_data_safe("titan_data/mass_fractions/case_c_exsolved_species.csv", 
                                       "Case C mass fractions")

cat("\n✓ Essential datasets loaded\n")
```

# Time Series Analysis

## Hydrosphere Comparison Across All Cases

```{r hydrosphere_comparison_plot}
if (!is.null(hydrosphere_comparison)) {
  
  # Create the main comparison plot
  hydrosphere_plot <- ggplot(hydrosphere_comparison, aes(x = Gyr, y = sum, color = source, shape = source)) +
    geom_point(size = 3) +
    geom_line(linewidth = 1) +
    scale_color_manual(values = case_colors) +
    scale_shape_manual(values = case_shapes) +
    scale_y_continuous(
      breaks = seq(0, 4e22, by = 0.5e22),
      labels = function(x) sprintf("%.1f", x/1e22),
      name = expression(paste("Total extracted mass (×10"^"22", " kg)")),
      sec.axis = sec_axis(~ . / TITAN_HYDROSPHERE_MASS * 100, name = "% Titan hydrosphere")
    ) +
    coord_cartesian(ylim = c(1e22, 4e22)) +
    labs(
      title = "Hydrosphere Mass Evolution: Comparison of All Cases",
      subtitle = "Total extracted mass over time for different starting compositions",
      x = "Time (Gyr)",
      caption = "Case A: Chondritic | Case B: Mixed | Case C: Cometary"
    ) +
    theme(
      legend.position = c(0.8, 0.3),
      legend.title = element_blank(),
      legend.background = element_rect(fill = "white"),
      plot.title = element_text(size = 14, face = "bold"),
      plot.subtitle = element_text(size = 12)
    )
  
  print(hydrosphere_plot)
  
  # Summary statistics
  final_masses <- hydrosphere_comparison %>%
    group_by(source) %>%
    filter(Gyr == max(Gyr)) %>%
    mutate(percent_hydrosphere = sum / TITAN_HYDROSPHERE_MASS * 100)
  
  cat("\n📊 Final Hydrosphere Masses:\n")
  for (i in 1:nrow(final_masses)) {
    cat(sprintf("   %s: %.2f × 10²² kg (%.1f%% of reference hydrosphere)\n", 
                final_masses$source[i], 
                final_masses$sum[i]/1e22, 
                final_masses$percent_hydrosphere[i]))
  }
  
} else {
  cat("❌ Cannot create hydrosphere comparison plot - data not available\n")
}
```

## Individual Element Extraction Time Series

```{r element_timeseries_functions}
#' Create time series plot for a specific element across all cases
#' 
#' @param element Element name
#' @param data_type "extracted" or "retained"
#' @param log_scale Use log scale for y-axis
#' @return ggplot object
plot_element_timeseries <- function(element, data_type = "extracted", log_scale = TRUE) {
  
  # Load data for all cases
  case_a_data <- load_timeseries("case_a", element, data_type)
  case_b_data <- load_timeseries("case_b", element, data_type)
  case_c_data <- load_timeseries("case_c", element, data_type)
  
  # Combine data
  combined_data <- bind_rows(
    if (!is.null(case_a_data)) mutate(case_a_data, case = "Case A"),
    if (!is.null(case_b_data)) mutate(case_b_data, case = "Case B"),
    if (!is.null(case_c_data)) mutate(case_c_data, case = "Case C")
  )
  
  if (nrow(combined_data) == 0) {
    cat(sprintf("⚠️  No data available for %s %s\n", element, data_type))
    return(NULL)
  }
  
  # Create plot
  p <- ggplot(combined_data, aes(x = Gyr, y = value, color = case, shape = case)) +
    geom_point(size = 2) +
    geom_line(linewidth = 1) +
    scale_color_manual(values = case_colors) +
    scale_shape_manual(values = case_shapes) +
    labs(
      title = sprintf("%s %s Over Time", str_to_title(element), str_to_title(data_type)),
      subtitle = "Comparison across all compositional cases",
      x = "Time (Gyr)",
      y = ifelse(data_type == "density_evolution", "Density (kg/m³)", "Mass (kg)")
    ) +
    theme(
      legend.position = "bottom",
      legend.title = element_blank()
    )
  
  # Apply log scale if requested
  if (log_scale && data_type != "density_evolution") {
    p <- p + scale_y_log10(labels = scales::scientific)
  }
  
  return(p)
}

#' Create a grid of element time series plots
#' 
#' @param elements Vector of element names
#' @param data_type "extracted" or "retained"
#' @param ncol Number of columns in grid
#' @return Combined ggplot object
plot_element_grid <- function(elements, data_type = "extracted", ncol = 3) {
  
  plots <- map(elements, ~ plot_element_timeseries(.x, data_type, log_scale = TRUE))
  
  # Remove NULL plots
  plots <- plots[!map_lgl(plots, is.null)]
  
  if (length(plots) == 0) {
    cat(sprintf("❌ No plots could be created for %s elements\n", data_type))
    return(NULL)
  }
  
  # Combine plots
  combined_plot <- wrap_plots(plots, ncol = ncol)
  
  return(combined_plot)
}

cat("✓ Element time series plotting functions defined\n")
```

```{r major_elements_timeseries}
# Plot major elements extraction
major_elements <- c("total_mass", "h", "c", "o", "n")

major_elements_plot <- plot_element_grid(major_elements, "extracted", ncol = 2)

if (!is.null(major_elements_plot)) {
  print(major_elements_plot)
} else {
  cat("❌ Could not create major elements plot\n")
}
```

```{r trace_elements_timeseries}
# Plot trace elements extraction
trace_elements <- c("na", "mg", "al", "si", "s", "k", "ca", "fe")

trace_elements_plot <- plot_element_grid(trace_elements, "extracted", ncol = 3)

if (!is.null(trace_elements_plot)) {
  print(trace_elements_plot)
} else {
  cat("❌ Could not create trace elements plot\n")
}
```

## Density Evolution

```{r density_evolution}
# Create density evolution plot
density_plots <- map(c("case_a", "case_b", "case_c"), function(case_name) {
  density_data <- load_timeseries(case_name, "density", "density_evolution")
  
  if (!is.null(density_data)) {
    case_label <- str_to_title(str_replace(case_name, "_", " "))
    
    ggplot(density_data, aes(x = Gyr, y = value)) +
      geom_point(size = 2, color = case_colors[case_label]) +
      geom_line(linewidth = 1, color = case_colors[case_label]) +
      labs(
        title = case_label,
        x = "Time (Gyr)",
        y = "Density (kg/m³)"
      ) +
      theme(
        plot.title = element_text(hjust = 0.5)
      )
  }
})

# Remove NULL plots
density_plots <- density_plots[!map_lgl(density_plots, is.null)]

if (length(density_plots) > 0) {
  combined_density_plot <- wrap_plots(density_plots, ncol = 3) +
    plot_annotation(
      title = "Mean Mantle Density Evolution",
      subtitle = "Average density changes over time for each compositional case"
    )
  
  print(combined_density_plot)
} else {
  cat("❌ Could not create density evolution plots\n")
}
```

# Phase Assemblage Analysis

## Phase Diagram Functions

```{r phase_diagram_functions}
#' Create phase assemblage plot for a specific case and pressure
#' 
#' @param case_name Case identifier
#' @param pressure Pressure level
#' @param thermal_profiles Optional thermal profiles to overlay
#' @return ggplot object
plot_phase_assemblage <- function(case_name, pressure, thermal_profiles = NULL) {
  
  # Load phase data
  phase_data <- load_phases(case_name, pressure)
  
  if (is.null(phase_data)) {
    return(NULL)
  }
  
  # Calculate pressure value for title
  pressure_value <- case_when(
    pressure == "1GPa" ~ "1 GPa",
    pressure == "2500MPa" ~ "2.5 GPa", 
    pressure == "3240MPa" ~ "3.24 GPa",
    pressure == "5GPa" ~ "5 GPa",
    TRUE ~ pressure
  )
  
  # Calculate radius estimate (approximate)
  radius_km <- case_when(
    pressure == "1GPa" ~ 1964,
    pressure == "2500MPa" ~ 1554,
    pressure == "3240MPa" ~ 1304, 
    pressure == "5GPa" ~ 83,
    TRUE ~ NA_real_
  )
  
  # Define phase colors - expanded palette
  phase_colors <- c(
    "Amph" = "#89C5DA", "Atg" = "#DA5724", "Bt" = "#599861", "Carb" = "#74D944", 
    "Chl" = "#CE50CA", "Cpx" = "#3F4921", "Deer" = "#7FDCC0", "Dia" = "#C0717C", 
    "Gth" = "#CBD588", "Gph" = "#5F7FC7", "Liz" = "#673770", "Mica" = "#C84248", 
    "Ol" = "#D3D93E", "Opx" = "#38333E", "Pu" = "#508578", "Py" = "#D7C1B1", 
    "Po" = "#8569D5", "Sp" = "#AD6F3B", "Tlc" = "#CD9BCD", "Tro" = "#6DDE88",
    "Rgw" = "#5E738F", "Gt" = "#D1A33D", "Grn" = "#8A7C64", "Law" = "#FF5733"
  )
  
  # Calculate maximum mass percent for scaling
  max_mass_percent <- phase_data %>%
    group_by(T_K) %>%
    summarise(total = sum(mass_percent, na.rm = TRUE), .groups = 'drop') %>%
    pull(total) %>%
    max(na.rm = TRUE)
  
  # Create base plot
  p <- ggplot(phase_data, aes(x = T_K, y = mass_percent, fill = phase)) +
    geom_area() +
    scale_fill_manual(values = phase_colors, na.value = "grey50") +
    scale_x_continuous(expand = c(0, 0)) +
    scale_y_continuous(expand = c(0, 0)) +
    coord_cartesian(ylim = c(0, max_mass_percent)) +
    labs(
      title = sprintf("%s at %s", str_to_title(str_replace(case_name, "_", " ")), pressure_value),
      subtitle = ifelse(!is.na(radius_km), sprintf("r ≈ %d km", radius_km), ""),
      x = "Temperature (K)",
      y = "Weight %",
      fill = "Phases"
    ) +
    theme(
      legend.position = "right",
      plot.title = element_text(hjust = 0.5),
      plot.subtitle = element_text(hjust = 0.5)
    ) +
    guides(fill = guide_legend(ncol = 3))
  
  # Add thermal profile lines if provided
  if (!is.null(thermal_profiles)) {
    # Add thermal evolution lines
    for (i in seq_along(thermal_profiles)) {
      profile <- thermal_profiles[[i]]
      # This would need the thermal profile data at the specific pressure
      # For now, we'll skip this feature
    }
  }
  
  return(p)
}

#' Create phase assemblage comparison across cases for a specific pressure
#' 
#' @param pressure Pressure level
#' @return Combined ggplot object
plot_phase_comparison <- function(pressure) {
  
  cases <- c("case_a", "case_b", "case_c")
  
  plots <- map(cases, ~ plot_phase_assemblage(.x, pressure))
  
  # Remove NULL plots
  plots <- plots[!map_lgl(plots, is.null)]
  
  if (length(plots) == 0) {
    cat(sprintf("❌ No phase assemblage data available for %s\n", pressure))
    return(NULL)
  }
  
  # Combine plots
  combined_plot <- wrap_plots(plots, ncol = length(plots)) +
    plot_annotation(
      title = sprintf("Phase Assemblage Comparison at %s", 
                     case_when(
                       pressure == "1GPa" ~ "1 GPa",
                       pressure == "2500MPa" ~ "2.5 GPa",
                       pressure == "3240MPa" ~ "3.24 GPa", 
                       pressure == "5GPa" ~ "5 GPa",
                       TRUE ~ pressure
                     )),
      subtitle = "Temperature vs. weight% phase diagrams across compositional cases"
    )
  
  return(combined_plot)
}

cat("✓ Phase diagram functions defined\n")
```

## Phase Assemblages at Different Pressures

```{r phase_assemblages_1GPa}
# Phase assemblages at 1 GPa
phase_1GPa <- plot_phase_comparison("1GPa")

if (!is.null(phase_1GPa)) {
  print(phase_1GPa)
} else {
  cat("❌ Could not create 1 GPa phase comparison\n")
}
```

```{r phase_assemblages_2500MPa}
# Phase assemblages at 2.5 GPa
phase_2500MPa <- plot_phase_comparison("2500MPa")

if (!is.null(phase_2500MPa)) {
  print(phase_2500MPa)
} else {
  cat("❌ Could not create 2.5 GPa phase comparison\n")
}
```

```{r phase_assemblages_5GPa}
# Phase assemblages at 5 GPa
phase_5GPa <- plot_phase_comparison("5GPa")

if (!is.null(phase_5GPa)) {
  print(phase_5GPa)
} else {
  cat("❌ Could not create 5 GPa phase comparison\n")
}
```

# Density-Temperature-Radius-Time Profiles

```{r density_radius_time_functions}
#' Create density vs radius plots over time
#' 
#' @param case_name Case identifier
#' @return ggplot object
plot_density_radius_time <- function(case_name) {
  
  file_path <- sprintf("titan_data/density_evolution/%s_density_radius_time.csv", case_name)
  density_data <- load_data_safe(file_path, sprintf("%s density-radius-time", case_name))
  
  if (is.null(density_data)) {
    return(NULL)
  }
  
  # Create the plot
  case_label <- str_to_title(str_replace(case_name, "_", " "))
  
  p <- ggplot(density_data, aes(x = radius_km, y = Bulk_rs_density, color = factor(Gyr))) +
    geom_line(linewidth = 1) +
    scale_color_viridis_d(name = "Time (Gyr)") +
    labs(
      title = sprintf("%s: Density vs Radius Over Time", case_label),
      subtitle = "Interior density evolution with time",
      x = "Radius (km)",
      y = "Density (kg/m³)"
    ) +
    theme(
      legend.position = "right"
    )
  
  return(p)
}

#' Create temperature vs radius plots over time
#' 
#' @param case_name Case identifier 
#' @return ggplot object
plot_temperature_radius_time <- function(case_name) {
  
  file_path <- sprintf("titan_data/density_evolution/%s_density_radius_time.csv", case_name)
  density_data <- load_data_safe(file_path, sprintf("%s temperature-radius-time", case_name))
  
  if (is.null(density_data)) {
    return(NULL)
  }
  
  # Create the plot
  case_label <- str_to_title(str_replace(case_name, "_", " "))
  
  p <- ggplot(density_data, aes(x = radius_km, y = T_K, color = factor(Gyr))) +
    geom_line(linewidth = 1) +
    scale_color_viridis_d(name = "Time (Gyr)") +
    labs(
      title = sprintf("%s: Temperature vs Radius Over Time", case_label),
      subtitle = "Interior temperature evolution with time",
      x = "Radius (km)", 
      y = "Temperature (K)"
    ) +
    theme(
      legend.position = "right"
    )
  
  return(p)
}

cat("✓ Density-radius-time plotting functions defined\n")
```

```{r density_radius_profiles}
# Create density-radius-time profiles for all cases
cases <- c("case_a", "case_b", "case_c")

density_radius_plots <- map(cases, plot_density_radius_time)
names(density_radius_plots) <- cases

# Remove NULL plots
density_radius_plots <- density_radius_plots[!map_lgl(density_radius_plots, is.null)]

if (length(density_radius_plots) > 0) {
  combined_density_radius <- wrap_plots(density_radius_plots, ncol = 1)
  print(combined_density_radius)
} else {
  cat("❌ Could not create density-radius-time profiles\n")
}
```

```{r temperature_radius_profiles}
# Create temperature-radius-time profiles for all cases
temperature_radius_plots <- map(cases, plot_temperature_radius_time)
names(temperature_radius_plots) <- cases

# Remove NULL plots  
temperature_radius_plots <- temperature_radius_plots[!map_lgl(temperature_radius_plots, is.null)]

if (length(temperature_radius_plots) > 0) {
  combined_temperature_radius <- wrap_plots(temperature_radius_plots, ncol = 1)
  print(combined_temperature_radius)
} else {
  cat("❌ Could not create temperature-radius-time profiles\n")
}
```

# Fluid Chemistry Analysis

## Aqueous Solute Concentrations

```{r fluid_chemistry_functions}
#' Plot aqueous solute concentrations vs pressure
#' 
#' @param fluid_data Fluid chemistry dataframe
#' @param case_label Case label for plot title
#' @return ggplot object
plot_fluid_chemistry <- function(fluid_data, case_label) {
  
  if (is.null(fluid_data)) {
    return(NULL)
  }
  
  # Define custom colors for elements
  element_colors <- c(
    "Sum_C" = "#332288", "Sum_Si" = "#117733", "Sum_Ca" = "#44AA99", 
    "Sum_Mg" = "#88CCEE", "Sum_K" = "#DDCC77", "Sum_Na" = "#CC6677", 
    "Sum_S" = "#AA4499", "Sum_NH4." = "#882255"
  )
  
  # Get element columns (exclude Pressure_bar and pH)
  element_cols <- names(fluid_data)[!names(fluid_data) %in% c("Pressure_bar", "pH")]
  
  # Create long format data
  fluid_long <- fluid_data %>%
    select(Pressure_bar, pH, all_of(element_cols)) %>%
    pivot_longer(cols = all_of(element_cols), names_to = "Element", values_to = "Concentration") %>%
    filter(Concentration >= 1e-10)  # Filter out very low concentrations
  
  # Create element labels with proper formatting
  element_labels <- c(
    "Sum_C" = "ΣC", "Sum_Si" = "ΣSi", "Sum_Ca" = "ΣCa", 
    "Sum_Mg" = "ΣMg", "Sum_K" = "ΣK", "Sum_Na" = "ΣNa",
    "Sum_S" = "ΣS", "Sum_NH4." = "ΣN"
  )
  
  # Get concentration range for pH scaling
  log_conc_range <- range(log10(fluid_long$Concentration))
  ph_range <- range(fluid_data$pH, na.rm = TRUE)
  
  # Create the main plot
  p <- ggplot() +
    # Plot element concentrations
    geom_line(data = fluid_long, 
              aes(x = Pressure_bar, y = Concentration, color = Element, linetype = Element),
              linewidth = 1) +
    # Plot pH on secondary axis
    geom_line(data = fluid_data,
              aes(x = Pressure_bar, 
                  y = 10^(log_conc_range[1] + (log_conc_range[2] - log_conc_range[1]) * 
                           (pH - ph_range[1]) / (ph_range[2] - ph_range[1]))),
              color = "#882255", linewidth = 1.2) +
    scale_y_log10(
      name = "Concentration (mol/kg)",
      limits = c(1e-10, NA),
      sec.axis = sec_axis(~ ph_range[1] + (ph_range[2] - ph_range[1]) * 
                           (log10(.) - log_conc_range[1]) / (log_conc_range[2] - log_conc_range[1]), 
                         name = "pH")
    ) +
    scale_x_log10(name = "Pressure (bar)") +
    scale_color_manual(values = element_colors, labels = element_labels) +
    scale_linetype_discrete(labels = element_labels) +
    labs(
      title = sprintf("%s: Aqueous Solute Concentrations", case_label),
      subtitle = "Elemental concentrations and pH vs pressure",
      color = "Element",
      linetype = "Element"
    ) +
    theme(
      legend.position = "bottom",
      panel.grid.minor = element_blank()
    )
  
  return(p)
}

cat("✓ Fluid chemistry plotting functions defined\n")
```

```{r fluid_chemistry_plots}
# Plot fluid chemistry for available cases
fluid_plots <- list()

if (!is.null(case_a_fluid)) {
  fluid_plots[["Case A"]] <- plot_fluid_chemistry(case_a_fluid, "Case A")
}

if (!is.null(case_c_fluid)) {
  fluid_plots[["Case C"]] <- plot_fluid_chemistry(case_c_fluid, "Case C")
}

# Display plots
if (length(fluid_plots) > 0) {
  for (case_name in names(fluid_plots)) {
    if (!is.null(fluid_plots[[case_name]])) {
      print(fluid_plots[[case_name]])
    }
  }
} else {
  cat("❌ No fluid chemistry data available\n")
}
```

# Mass Fraction Analysis

## Exsolved Species Composition

```{r mass_fraction_functions}
#' Create lollipop chart of exsolved species mass fractions
#' 
#' @param mass_frac_data Mass fraction dataframe
#' @param case_label Case label for plot title
#' @return ggplot object
plot_mass_fractions <- function(mass_frac_data, case_label) {
  
  if (is.null(mass_frac_data) || nrow(mass_frac_data) == 0) {
    return(NULL)
  }
  
  # Filter and order data
  filtered_data <- mass_frac_data %>%
    filter(mass_frac >= 1e-10) %>%
    arrange(desc(mass_frac))
  
  # Create formatted labels for species
  format_species_label <- function(species) {
    case_when(
      species == "H2O_mass_frac" ~ "H₂O (g)",
      species == "CH4_mass_frac" ~ "CH₄ (g)",
      species == "H2_mass_frac" ~ "H₂ (g)",
      species == "H2S_mass_frac" ~ "H₂S (g)",
      species == "NH3_mass_frac" ~ "NH₃ (g)",
      species == "CO2_mass_frac" ~ "CO₂ (g)",
      species == "CO_mass_frac" ~ "CO (g)",
      species == "N2_mass_frac" ~ "N₂ (g)",
      species == "S2_mass_frac" ~ "S₂ (g)",
      species == "SO2_mass_frac" ~ "SO₂ (g)",
      species == "C2H6_mass_frac" ~ "C₂H₆ (g)",
      species == "lizardite_kg" ~ "Lizardite",
      species == "lizardite_mass_frac" ~ "Lizardite",
      species == "diopside_mass_frac" ~ "Diopside",
      species == "wollastonite_mass_frac" ~ "Wollastonite",
      species == "aragonite_mass_frac" ~ "Aragonite",
      species == "calcite_mass_frac" ~ "Calcite",
      TRUE ~ str_replace_all(species, "_mass_frac|_kg", "")
    )
  }
  
  # Apply formatting
  filtered_data$species_formatted <- format_species_label(filtered_data$species_mass_frac)
  filtered_data$species_formatted <- factor(filtered_data$species_formatted, 
                                            levels = rev(filtered_data$species_formatted))
  
  # Determine colors (gases vs solids)
  filtered_data$type <- ifelse(grepl("\\(g\\)", filtered_data$species_formatted), "Gas", "Solid")
  
  type_colors <- c("Gas" = "#1E88E5", "Solid" = "#FFC107")
  
  # Create lollipop chart
  p <- ggplot(filtered_data, aes(x = mass_frac, y = species_formatted, color = type)) +
    geom_segment(aes(x = 1e-10, xend = mass_frac, yend = species_formatted), 
                 color = "grey70", linewidth = 0.5) +
    geom_point(size = 4) +
    scale_x_log10(
      breaks = scales::trans_breaks("log10", function(x) 10^x),
      labels = scales::trans_format("log10", scales::math_format(10^.x))
    ) +
    scale_color_manual(values = type_colors) +
    labs(
      title = sprintf("%s: Exsolved Species Mass Fractions", case_label),
      subtitle = "Mass fraction of extracted hydrosphere",
      x = "Mass Fraction",
      y = "",
      color = "Type"
    ) +
    theme(
      legend.position = "bottom",
      panel.grid.minor = element_blank(),
      axis.text.y = element_text(size = 10)
    )
  
  return(p)
}

cat("✓ Mass fraction plotting functions defined\n")
```

```{r mass_fraction_plots}
# Plot mass fractions for available cases
mass_frac_plots <- list()

if (!is.null(case_a_mass_fractions)) {
  mass_frac_plots[["Case A"]] <- plot_mass_fractions(case_a_mass_fractions, "Case A")
}

if (!is.null(case_c_mass_fractions)) {
  mass_frac_plots[["Case C"]] <- plot_mass_fractions(case_c_mass_fractions, "Case C")
}

# Display plots
if (length(mass_frac_plots) > 0) {
  for (case_name in names(mass_frac_plots)) {
    if (!is.null(mass_frac_plots[[case_name]])) {
      print(mass_frac_plots[[case_name]])
    }
  }
} else {
  cat("❌ No mass fraction data available\n")
}
```

# Custom Analysis Functions

## Retained vs Extracted Mass Ratios (Case B & C)

```{r retention_analysis}
#' Calculate and plot retention ratios for specific elements
#' 
#' @param case_name Case identifier
#' @param elements Vector of element names
#' @return ggplot object
plot_retention_ratios <- function(case_name, elements = c("c", "h", "k", "si")) {
  
  # Load retained and extracted data for each element
  retention_data <- map_dfr(elements, function(element) {
    retained_data <- load_timeseries(case_name, element, "retained")
    extracted_data <- load_timeseries(case_name, element, "extracted")
    
    if (!is.null(retained_data) && !is.null(extracted_data)) {
      # Merge by time
      combined <- inner_join(
        retained_data %>% select(Gyr, retained = value),
        extracted_data %>% select(Gyr, extracted = value),
        by = "Gyr"
      ) %>%
        mutate(
          ratio = retained / (retained + extracted),
          element = toupper(element)
        )
      
      return(combined)
    }
    
    return(NULL)
  })
  
  if (nrow(retention_data) == 0) {
    cat(sprintf("⚠️  No retention data available for %s\n", case_name))
    return(NULL)
  }
  
  # Create plot
  case_label <- str_to_title(str_replace(case_name, "_", " "))
  
  p <- ggplot(retention_data, aes(x = Gyr, y = ratio, color = element)) +
    geom_point(size = 2) +
    geom_line(linewidth = 1) +
    scale_y_continuous(labels = scales::percent, limits = c(0, 1)) +
    scale_color_viridis_d() +
    labs(
      title = sprintf("%s: Element Retention Ratios", case_label),
      subtitle = "Fraction of each element retained in the solid residue",
      x = "Time (Gyr)",
      y = "Retention Ratio",
      color = "Element"
    ) +
    theme(
      legend.position = "bottom"
    )
  
  return(p)
}

# Plot retention ratios for Cases B and C
retention_plots <- map(c("case_b", "case_c"), plot_retention_ratios)
names(retention_plots) <- c("Case B", "Case C")

# Remove NULL plots
retention_plots <- retention_plots[!map_lgl(retention_plots, is.null)]

if (length(retention_plots) > 0) {
  combined_retention <- wrap_plots(retention_plots, ncol = 2) +
    plot_annotation(
      title = "Element Retention Analysis",
      subtitle = "Comparison of retention vs extraction for key elements"
    )
  
  print(combined_retention)
} else {
  cat("❌ Could not create retention ratio plots\n")
}
```

# Summary Analysis

## Key Findings Summary

```{r summary_analysis}
cat("📋 SUMMARY OF KEY FINDINGS\n")
cat("========================\n\n")

# Hydrosphere mass summary
if (!is.null(hydrosphere_comparison)) {
  final_masses <- hydrosphere_comparison %>%
    group_by(source) %>%
    filter(Gyr == max(Gyr)) %>%
    mutate(percent_hydrosphere = sum / TITAN_HYDROSPHERE_MASS * 100)
  
  cat("🌊 FINAL HYDROSPHERE MASSES (at 4.5 Gyr):\n")
  for (i in 1:nrow(final_masses)) {
    cat(sprintf("   %s: %.2f × 10²² kg (%.1f%% of reference)\n", 
                final_masses$source[i], 
                final_masses$sum[i]/1e22, 
                final_masses$percent_hydrosphere[i]))
  }
  cat("\n")
}

# Phase assemblage summary
cat("🔬 PHASE ASSEMBLAGE COVERAGE:\n")
pressure_levels <- c("1GPa", "2500MPa", "3240MPa", "5GPa")
cases <- c("case_a", "case_b", "case_c")

for (pressure in pressure_levels) {
  available_cases <- map_lgl(cases, ~ !is.null(load_phases(.x, pressure)))
  n_available <- sum(available_cases)
  cat(sprintf("   %s: %d/3 cases available\n", pressure, n_available))
}
cat("\n")

# Fluid chemistry summary
cat("💧 FLUID CHEMISTRY DATA:\n")
cat(sprintf("   Case A: %s\n", ifelse(!is.null(case_a_fluid), "✓ Available", "❌ Not available")))
cat(sprintf("   Case B: %s\n", "❌ Not available (expected)"))
cat(sprintf("   Case C: %s\n", ifelse(!is.null(case_c_fluid), "✓ Available", "❌ Not available")))
cat("\n")

# Data completeness summary
cat("📊 DATASET COMPLETENESS:\n")

# Count available timeseries
timeseries_elements <- c("total_mass", "h", "c", "o", "n", "na", "mg", "al", "si", "s", "k", "ca", "fe")
extraction_coverage <- map_dfr(cases, function(case) {
  available <- map_lgl(timeseries_elements, ~ !is.null(load_timeseries(case, .x, "extracted")))
  tibble(
    case = str_to_title(str_replace(case, "_", " ")),
    available_elements = sum(available),
    total_elements = length(timeseries_elements),
    coverage_percent = round(sum(available) / length(timeseries_elements) * 100, 1)
  )
})

for (i in 1:nrow(extraction_coverage)) {
  cat(sprintf("   %s: %d/%d elements (%.1f%%)\n", 
              extraction_coverage$case[i],
              extraction_coverage$available_elements[i],
              extraction_coverage$total_elements[i],
              extraction_coverage$coverage_percent[i]))
}

cat("\n")
cat("✅ Analysis complete! All available data has been visualized.\n")
```

# Interactive Exploration Guide

## How to Use This Notebook

```{r usage_guide, echo=FALSE, results='asis'}
cat("## 🚀 Getting Started\n\n")

cat("This notebook provides comprehensive visualization tools for analyzing Titan interior evolution models. Here's how to explore the data:\n\n")

cat("### 📁 Data Organization\n")
cat("- **Time Series**: Element extraction/retention over 4.5 billion years\n")
cat("- **Phase Assemblages**: Mineral stability diagrams at different pressures\n")
cat("- **Density Evolution**: Interior structure changes over time\n")
cat("- **Fluid Chemistry**: Aqueous solute concentrations (Cases A & C)\n")
cat("- **Mass Fractions**: Composition of exsolved gases and minerals\n\n")

cat("### 🎯 Key Comparisons\n")
cat("- **Case A (Chondritic)**: Pure CI chondrite composition\n")
cat("- **Case B (Mixed)**: 26% chondritic + 74% cometary\n") 
cat("- **Case C (Cometary)**: Pure cometary composition\n\n")

cat("### 🔧 Customization Options\n")
cat("You can modify the analysis by:\n")
cat("1. **Changing elements**: Edit the `elements` vectors in plotting functions\n")
cat("2. **Adjusting time ranges**: Filter data by `Gyr` column\n")
cat("3. **Modifying pressure levels**: Select different pressure values for phase diagrams\n")
cat("4. **Custom visualizations**: Use the provided functions as templates\n\n")

cat("### 📊 Reproducing Key Figures\n")
cat("- **Main hydrosphere plot**: Already generated in the hydrosphere comparison section\n")
cat("- **Phase assemblage grids**: Use `plot_phase_comparison()` function\n")
cat("- **Element time series**: Use `plot_element_timeseries()` function\n")
cat("- **Density evolution**: Use `plot_density_radius_time()` function\n\n")

cat("### 💡 Pro Tips\n")
cat("- Check the console output for data loading status\n")
cat("- Missing plots usually indicate missing source data\n")
cat("- Use `View()` function to inspect data frames directly\n")
cat("- Modify color schemes by editing the `case_colors` vector\n\n")
```

## Available Data Files Reference

```{r data_reference, echo=FALSE, results='asis'}
if (file.exists("titan_data/data_manifest.csv")) {
  manifest <- read.csv("titan_data/data_manifest.csv", stringsAsFactors = FALSE)
  
  cat("## 📚 Complete Data Reference\n\n")
  
  # Group by dataset type
  by_type <- split(manifest, manifest$dataset_type)
  
  for (type in names(by_type)) {
    cat(sprintf("### %s\n", type))
    type_data <- by_type[[type]]
    
    for (i in 1:nrow(type_data)) {
      cat(sprintf("- **%s** (%s)\n", basename(type_data$file_path[i]), type_data$case[i]))
      cat(sprintf("  - Columns: %s\n", type_data$columns[i]))
      cat(sprintf("  - Rows: %s\n", type_data$rows_approx[i]))
    }
    cat("\n")
  }
} else {
  cat("⚠️ Data manifest not found. Run the data extraction script first.\n")
}
```

---

## Session Information

```{r session_info}
cat("R Session Information:\n")
cat("=====================\n")
sessionInfo()
```

---

**📝 Citation**: If you use this toolkit for research, please cite the original Titan interior evolution study and acknowledge the data sources.

**🔗 Source Code**: This visualization toolkit is designed to be modular and extensible. Feel free to adapt the functions for your own analysis needs.

**📧 Support**: For questions about the data or visualization methods, please refer to the original research publication.