---
title: "Heckman_in_time"
author: "Arthur Vienne"
date: "2024-10-29"
output: html_document
---

#Heckman measurement 1: 130 days
```{r}
setwd ("C:/Users/AVienne/OneDrive - Universiteit Antwerpen/Labexp 3 (start Oct 2022)/PAPER_")
library (readxl)
library (dplyr)
data_allpools <- data.frame()
pools <- c("pyrophosphate","NH2OH","Dithionite")
empty_list <- list () 

for (j in (1:length (pools))) { 
  pool <- pools[j]
  data <- read_excel("SequentialExtraction_dataset_heckman.xlsx", sheet = pool, col_names = T)
  data$pool <- pool 
  data <- data [-which (data$Treatment=="FB1Bac"|data$Treatment=="C1Bac"),]
  blankdata <- data[which (data$Treatment=="BLANK"),]
  data <- data[-which (data$Treatment=="BLANK"),]
  
  #Blank normalization
  data$Fe_fraction2  <- data$Fe_fraction2 - mean (blankdata$Fe_fraction2)
  data$K_fraction2  <- data$K_fraction2 - mean (blankdata$Fe_fraction2)
  data$Na_fraction2  <- data$Na_fraction2 - mean (blankdata$Na_fraction2)
   data$Mg_fraction2 <- data$Mg_fraction2 - mean (blankdata$Fe_fraction2)
  data$Ca_fraction2  <- data$Ca_fraction2 - mean (blankdata$Ca_fraction2)
  data$Al_fraction2  <- data$Al_fraction2 - mean (blankdata$Al_fraction2)
  data$Si_fraction2  <- data$Si_fraction2 - mean (blankdata$Si_fraction2)
  
  data$Fe <- data$Fe_fraction2* data$V_mL_fraction2/data$sampleMass_g #µg/g
  data$K <- data$K_fraction2* data$V_mL_fraction2/data$sampleMass_g #µg/g
  data$Mg<- data$Mg_fraction2* data$V_mL_fraction2/data$sampleMass_g #µg/g
  data$Ca<- data$Ca_fraction2* data$V_mL_fraction2/data$sampleMass_g #µg/g
  data$Al<- data$Al_fraction2* data$V_mL_fraction2/data$sampleMass_g #µg/g
  data$Si<- data$Si_fraction2* data$V_mL_fraction2/data$sampleMass_g #µg/g
  data$Na<- data$Na_fraction2* data$V_mL_fraction2/data$sampleMass_g #µg/g
  
  
  Puresilicatesdata <- data[is.na(data$soilmass), ]
  data <- data[!is.na(data$soilmass), ]
  Treatments <- unique (data$Treatment)
  Treatments_letters <- unique(gsub("[0-9]", "", Treatments))
  Treatments <- Treatments_letters [-which (Treatments_letters=="C")]
  
  #rock/soil ratio? 
  combined_data <- data.frame()
  
  for (i in (1: length (Treatments))) {
    # data$ Dunite_soil <- 0.11 / data$soilmass #g/g
    data$Basalt_soil <- 1.1 / data$soilmass #g/g
    treatment <- Treatments[i]
    corrected_data <- data[grep(paste0("^",treatment), data$Treatment), ]
    
    #to do add if to use dunite or basalt/soil ratio if the rock is dunite or basalt 
    
    #correction for all elements:
    if (treatment == "FB" |treatment == "CB" ){  
      corrected_data$Ca <- corrected_data$Ca - corrected_data$Basalt_soil * mean (Puresilicatesdata [which (Puresilicatesdata$Treatment==treatment),]$Ca)
      corrected_data$Fe <- corrected_data$Fe - corrected_data$Basalt_soil * mean (Puresilicatesdata [which (Puresilicatesdata$Treatment==treatment),]$Fe)
      corrected_data$K <- corrected_data$K - corrected_data$Basalt_soil * mean (Puresilicatesdata [which (Puresilicatesdata$Treatment==treatment),]$K)
      corrected_data$Mg <- corrected_data$Mg - corrected_data$Basalt_soil * mean (Puresilicatesdata [which (Puresilicatesdata$Treatment==treatment),]$Mg)
      corrected_data$Na <- corrected_data$Na - corrected_data$Basalt_soil * mean (Puresilicatesdata [which (Puresilicatesdata$Treatment==treatment),]$Na)
      corrected_data$Al <- corrected_data$Al - corrected_data$Basalt_soil * mean (Puresilicatesdata [which (Puresilicatesdata$Treatment==treatment),]$Al)
      corrected_data$Si <- corrected_data$Si - corrected_data$Basalt_soil * mean (Puresilicatesdata [which (Puresilicatesdata$Treatment==treatment),]$Si)
    }
    combined_data <- rbind(combined_data, corrected_data)
  }
  
  data_control <- data[which (data$Treatment=="C3"|data$Treatment=="C1" ),]
  combined_data <- rbind(combined_data, data_control)
  
  #elements in combined_data are in µg/g --> µmol/g
  combined_data$Fe_mol <- combined_data$Fe * combined_data$soilmass * 1000/ 55.845 /10^(6)
  combined_data$K_mol <- combined_data$K   * combined_data$soilmass * 1000/ 39.098/10^(6)
  combined_data$Ca_mol <- combined_data$Ca  * combined_data$soilmass * 1000/ 40.078/10^(6)
  combined_data$Mg_mol <- combined_data$Mg  * combined_data$soilmass * 1000/ 24.305/10^(6)
  combined_data$Na_mol <- combined_data$Na  * combined_data$soilmass * 1000/ 22.98976/10^(6)
  combined_data$Al_mol <- combined_data$Al  * combined_data$soilmass * 1000/ 26.98153/10^(6)
  combined_data$Si_mol <- combined_data$Si  * combined_data$soilmass * 1000/ 26.98153/10^(6)
    #units in mol element/mesocosm 
  data_allpools <- rbind(data_allpools,combined_data)
}

#replace <0 values (smaller than blank) by 0
data_allpools <- data_allpools %>%
  mutate(across(c(Fe_mol, Al_mol, Ca_mol, Mg_mol, K_mol, Na_mol, Si_mol), ~ ifelse(. < 0, 0, .)))
#OK, now normalize for controls
elements <- c("Fe","Al","Ca","Mg","Si")
plot_list <- list()

for (j in (1: length(elements))) {
  element <- elements[j]
  
  #subtract controls for each element
  library (Rmisc)
  library (ggplot2)
  
  measurevar <- paste0(element, "_mol")
  
  # Apply the summarySE function using the dynamically created measurevar
  sum <- summarySE(data_allpools[complete.cases(data_allpools[[measurevar]]), ], 
                   measurevar = measurevar, 
                   groupvars = c("Treatment", "pool"))
  
  
  # Initialize a list to store results
  results_list <- list()
  
  # Define a function to calculate difference and propagate standard error
  calculate_difference <- function(SIC1, SE1, SIC2, SE2) {
    diff <- SIC1 - SIC2
    propagated_se <- sqrt(SE1^2 + SE2^2)
    return(c(diff, propagated_se))
  }
  
  # Subtract values for FB1 by C1 for each pool
  for (pool in unique(sum$pool)) {
    results_list[[paste("FB1 - C1", pool)]] <- calculate_difference(
      sum[[measurevar]][sum$Treatment == "FB1" & sum$pool == pool], 
      sum$se[sum$Treatment == "FB1" & sum$pool == pool], 
      sum[[measurevar]][sum$Treatment == "C1" & sum$pool == pool], 
      sum$se[sum$Treatment == "C1" & sum$pool == pool]
    )
  }
  
  # Loop through treatments ending with '3' and subtract from C3, for each pool
  for (treatment in sum$Treatment[grepl("3$", sum$Treatment) & sum$Treatment != "C3"]) {
    for (pool in unique(sum$pool)) {
      comparison_name <- paste(treatment, "- C3", pool)
      results_list[[comparison_name]] <- calculate_difference(
        sum[[measurevar]][sum$Treatment == treatment & sum$pool == pool] , 
        sum$se[sum$Treatment == treatment & sum$pool == pool], 
        sum[[measurevar]][sum$Treatment == "C3" & sum$pool == pool], 
        sum$se[sum$Treatment == "C3" & sum$pool == pool]
      )
    }
  }
  
  # Initialize an empty list to store the final rows
  results_data <- list()
  
  # Loop through the names of the results_list to construct a dataframe
  for (name in names(results_list)) {
    # Split the comparison name into treatment and pool
    name_split <- unlist(strsplit(name, " "))
    
    # Extract the treatment name (first part before subtraction)
    treatment <- name_split[1]  # e.g., "FB1", "CB3", etc.
    
    # Extract the pool name (third element in the split string)
    name_spli2 <- unlist(strsplit(name, " "))
    
    pool <- name_split[4]  # e.g., "Exchangeable", "Carbonate", etc.
    
    # Extract the difference and propagated SE
    diff <- results_list[[name]][1]
    se <- results_list[[name]][2]
    
    # Append the row to the results_data list
    results_data[[name]] <- data.frame(
      diff = diff,
      se = se,
      treatment = treatment,
      pool = pool
    )
  }
  
  # Combine all rows into a single dataframe
  results_df <- do.call(rbind, results_data)
  
  
  library(dplyr)
  # Create a summary dataframe by treatment, summing the 'diff' and calculating propagated 'se'
  total_rows <- results_df %>%
    group_by(treatment) %>%
    summarize(
      diff = sum(diff, na.rm = TRUE),  # Sum the 'diff' for each treatment
      se = sqrt(sum(se^2, na.rm = TRUE))  # Calculate the propagated SE
    ) %>%
    mutate(pool = "total")  # Add a column for pool as "total"
  # Bind the original dataframe with the new 'total' rows
  results_df_with_total <- bind_rows(results_df, total_rows)
  
  
  #make a ggplot 
  # Assign the plot to the dynamically created name
  treatment_order <- c( "FB1", "CB3","FB3","CD3","MD3","FD3")
  # Convert Treatment column to factor with the specified order
  results_df_with_total$treatment <- factor(results_df_with_total$treatment, levels = treatment_order)
    results_df_with_total1 <- results_df_with_total #Store under different name
    results_df_with_total1$element <- element
empty_list[[j]] <- results_df_with_total1

  
  library (RColorBrewer)
  cols = brewer.pal(n = 8, name = 'Set2')

  #now add simulation data!!
  setwd ("C:/Users/AVienne/OneDrive - Universiteit Antwerpen/Labexp 3 (start Oct 2022)/PAPER_")
  treatmentvector <- c("FB1","FB3","CB3")
  #replace by for loop for i once all the data is there
 
  plot_name <- paste(element, "plot", sep = "")  # Create the plot name dynamically
  #make the combined plot now
  assign(plot_name, ggplot(results_df_with_total, aes(x = treatment, y = diff, color = pool)) +
           geom_point(size = 3,position = position_dodge(width = 0.5)) +  # Plot points for 'diff' values
           geom_errorbar(aes(ymin = diff - se, ymax = diff + se), width = 0.2, position = position_dodge(width = 0.5)) +  # Add error bars
           theme_bw() +  # Minimal theme for cleaner look
           theme(axis.title.y = element_text(size= 10), 
                 axis.title.x = element_text(size = 10), 
                 axis.title = element_text(face = "bold"),
                 panel.grid.major = element_blank(), 
                 panel.grid.minor = element_blank(), 
                 panel.background = element_blank(), 
                 axis.line = element_line(colour = "black")) +
           scale_color_manual(values = cols) +  # Set colors from cols
           labs(
             x = "Treatment", 
             y = paste("\u0394", "mol", element), 
             color = "Pool"
           ) +  theme(legend.position = "none") )
  
  plot_list[[j]]  <- get (plot_name)
  
  
  
  for (i in (1: length(treatmentvector))){  
    treatment <- treatmentvector[i]
    
    data <- read_excel(paste0("Simulationdata_",treatment,".xlsx"), sheet = "Sheet1", col_names = TRUE)
    #Data of day 130
    data130 <- data[which (data$days==130),]
    secondaryminerals <- c("gibbsite","goethite","calcite","Nesquehonite","CaMontmorillonite","Chlorite","Fe.OH.3.a.","Al.OH.3.a.","SiO2.a.")
    Alvector <- c(1,0,0,0,2.33,2,0,1,0)
    Cavector <- c(0,0,1,0,0.165,0,0,0,0)
    Mgvector <- c(0,0,0,1,0,5,0,0,0)
    Fevector <- c(0,1,0,0,0,0,1,0,0)
    Sivector <- c(0,0,0,0,3.67,3,0,0,1) #mol/mol
    
    for (k in (1:length (secondaryminerals))) {
      mineral <- secondaryminerals[k]
      diff <- data130[data130$type == "rock", mineral] - data130[data130$type == "control", mineral]
      #mol difference after 130 days
      
      Alvector[k] <- diff * Alvector[k]
      Cavector[k] <- diff * Cavector[k]
      Mgvector[k] <- diff * Mgvector[k]
      Fevector[k] <- diff * Fevector[k]
      Sivector[k] <- diff * Sivector[k] }
    
    results_df <- data.frame(Mineral = secondaryminerals,
                             Al = unlist(Alvector),
                             Ca = unlist(Cavector),
                             Mg = unlist(Mgvector),
                             Fe = unlist(Fevector),
                             Si = unlist(Sivector), stringsAsFactors = FALSE , treatment)
    
    if (!(element %in% c("Na", "K"))) {
      results_df_element <- results_df[, c(element, "Mineral","treatment"), drop = FALSE]
      results_df_element <-results_df_element[-which(results_df_element[,1]==0),]
      #calculate total simulated value
      total_value <- sum(results_df_element[[element]], na.rm = TRUE)
      total_row <- data.frame(
        element = total_value,      # Use dynamic column name
        Mineral = "total",             # Set the mineral as "total"
        treatment = unique(results_df_element$treatment))
      colnames (total_row)[1] <- element
      
      # Get the treatment from the existing data
      results_df_element <- rbind(results_df_element, total_row)
      available_shapes <- c(17, 18, 15, 19, 3, 4, 7, 20, 9)
      shape_mapping <- setNames(available_shapes[1:length(secondaryminerals)], secondaryminerals)
      shape_mapping["total"] <- 8
    }
    
    secondaryminerals <- c("gibbsite","goethite","calcite","Nesquehonite","CaMontmorillonite","Chlorite","Fe.OH.3.a.","Al.OH.3.a.","SiO2.a.")
    
    
    
    #trick for legend below:
    secondaryminerals <- c("gibbsite","goethite","calcite","Nesquehonite","CaMontmorillonite","Chlorite","Fe.OH.3.a.","Al.OH.3.a.","SiO2.a.")
    mineral_assignments <- rep(secondaryminerals, length.out = nrow(results_df_with_total))
    # Shuffle the mineral_assignments to create random order
    set.seed(123)  # Set seed for reproducibility
    mineral_assignments <- sample(mineral_assignments)
    mineral_assignments[1]<- "total"
    # Add the Minerals column to the dataframe
    results_df_with_total$Mineral <- mineral_assignments
    
    name <- "forlegend"
    available_shapes <- c(17, 18, 15, 19, 3, 4, 7, 20, 9)
    shape_mapping <- setNames(available_shapes[1:length(secondaryminerals)], secondaryminerals)
    shape_mapping["total"] <- 8
    
    assign(name, ggplot(results_df_with_total, aes(x = treatment, y = diff, color = pool)) +
             geom_point(size = 3,position = position_dodge(width = 0.5)) +  # Plot points for 'diff' values
             geom_errorbar(aes(ymin = diff - se, ymax = diff + se), width = 0.2, position = position_dodge(width = 0.5)) +  # Add error bars
             theme_bw() +  # Minimal theme for cleaner look
             theme(axis.title.y = element_text(size= 10), 
                   axis.title.x = element_text(size = 10), 
                   axis.title = element_text(face = "bold"),
                   panel.grid.major = element_blank(), 
                   panel.grid.minor = element_blank(), 
                   panel.background = element_blank(), 
                   axis.line = element_line(colour = "black")) +
             scale_color_manual(values = cols, name = "Pool \n (extracted)") +  # Set colors from cols
             labs(
               x = "Treatment", 
               y = paste("\u0394", "mol", element), 
               color = "Pool"
             )       + 
             geom_point(data = results_df_with_total, 
                        aes(x = treatment, y = diff, shape = Mineral),  # Correct dynamic column access
                        size = 3, fill = NA, col = "black") + 
             scale_shape_manual(values = shape_mapping, name = "Secondary minerals \n (simulated)"))
   
        # Access the plot from the list by index
      # Modify the plot by adding points
      plot_list[[j]] <- plot_list[[j]] + 
        geom_point(data = results_df_element, 
                   aes(x = treatment, y = .data[[element]], shape = Mineral),  # Correct dynamic column access
                   size = 3, fill = NA, col = "black") + 
        scale_shape_manual(values = shape_mapping)
  }
  
}

  dataframe_1 <- do.call(rbind, empty_list)

library (cowplot)
legend <- get_legend (forlegend)
final_plot <- plot_grid (plot_grid(plotlist = plot_list, labels = "AUTO"), legend, ncol = 2,  rel_widths  = c(1, 0.2))

output_path <- "C:/Users/AVienne/OneDrive - Universiteit Antwerpen/Labexp 3 (start Oct 2022)/PAPER_/Sec_minerals_heckman.jpeg"
ggsave(filename = output_path, plot = final_plot, width = 10, height = 8, units = "in", dpi = 600)

#to do 

  
```

#Heckman measurement 2: 386 days
```{r}
setwd ("C:/Users/AVienne/OneDrive - Universiteit Antwerpen/Labexp 3 (start Oct 2022)/PAPER_")
library (readxl)
library (dplyr)
data_allpools <- data.frame()
pools <- c("pyrophosphate","NH2OH","Dithionite")
empty_list <- list () 

for (j in (1:length (pools))) { 
  pool <- pools[j]
  data <- read_excel("SequentialExtraction_dataset_heckman2.xlsx", sheet = pool, col_names = T)
  data$pool <- pool 
  blankdata <- data[which (data$Treatment=="BLANK"),]
  data <- data[-which (data$Treatment=="BLANK"),]
  data <- data[-which (data$Treatment=="FD3" |data$Treatment=="MD3" | data$Treatment=="CD3"),]
  
  #Blank normalization
  data$Fe_fraction2  <- data$Fe_fraction2 - mean (blankdata$Fe_fraction2)
  data$K_fraction2  <- data$K_fraction2 - mean (blankdata$Fe_fraction2)
  data$Na_fraction2  <- data$Na_fraction2 - mean (blankdata$Na_fraction2)
   data$Mg_fraction2 <- data$Mg_fraction2 - mean (blankdata$Fe_fraction2)
  data$Ca_fraction2  <- data$Ca_fraction2 - mean (blankdata$Ca_fraction2)
  data$Al_fraction2  <- data$Al_fraction2 - mean (blankdata$Al_fraction2)
  data$Si_fraction2  <- data$Si_fraction2 - mean (blankdata$Si_fraction2)
  
  data$Fe <- data$Fe_fraction2* data$V_mL_fraction2/data$sampleMass_g #µg/g
  data$K <- data$K_fraction2* data$V_mL_fraction2/data$sampleMass_g #µg/g
  data$Mg<- data$Mg_fraction2* data$V_mL_fraction2/data$sampleMass_g #µg/g
  data$Ca<- data$Ca_fraction2* data$V_mL_fraction2/data$sampleMass_g #µg/g
  data$Al<- data$Al_fraction2* data$V_mL_fraction2/data$sampleMass_g #µg/g
  data$Si<- data$Si_fraction2* data$V_mL_fraction2/data$sampleMass_g #µg/g
  data$Na<- data$Na_fraction2* data$V_mL_fraction2/data$sampleMass_g #µg/g
  
  
  Puresilicatesdata <- data[is.na(data$soilmass), ]
  data <- data[!is.na(data$soilmass), ]
  Treatments <- unique (data$Treatment)
  Treatments_letters <- unique(gsub("[0-9]", "", Treatments))
  Treatments <- Treatments_letters [-which (Treatments_letters=="C")]
  
  #rock/soil ratio? 
  combined_data <- data.frame()
  
  for (i in (1: length (Treatments))) {
    # data$ Dunite_soil <- 0.11 / data$soilmass #g/g
    data$Basalt_soil <- 1.1 / data$soilmass #g/g
    treatment <- Treatments[i]
    corrected_data <- data[grep(paste0("^",treatment), data$Treatment), ]
    
    #to do add if to use dunite or basalt/soil ratio if the rock is dunite or basalt 
    
    #correction for all elements:
    if (treatment == "FB" |treatment == "CB" ){  
      corrected_data$Ca <- corrected_data$Ca - corrected_data$Basalt_soil * mean (Puresilicatesdata [which (Puresilicatesdata$Treatment==treatment),]$Ca)
      corrected_data$Fe <- corrected_data$Fe - corrected_data$Basalt_soil * mean (Puresilicatesdata [which (Puresilicatesdata$Treatment==treatment),]$Fe)
      corrected_data$K <- corrected_data$K - corrected_data$Basalt_soil * mean (Puresilicatesdata [which (Puresilicatesdata$Treatment==treatment),]$K)
      corrected_data$Mg <- corrected_data$Mg - corrected_data$Basalt_soil * mean (Puresilicatesdata [which (Puresilicatesdata$Treatment==treatment),]$Mg)
      corrected_data$Na <- corrected_data$Na - corrected_data$Basalt_soil * mean (Puresilicatesdata [which (Puresilicatesdata$Treatment==treatment),]$Na)
      corrected_data$Al <- corrected_data$Al - corrected_data$Basalt_soil * mean (Puresilicatesdata [which (Puresilicatesdata$Treatment==treatment),]$Al)
      corrected_data$Si <- corrected_data$Si - corrected_data$Basalt_soil * mean (Puresilicatesdata [which (Puresilicatesdata$Treatment==treatment),]$Si)
    }
    combined_data <- rbind(combined_data, corrected_data)
  }
  
  data_control <- data[which (data$Treatment=="C3"|data$Treatment=="C1" ),]
  combined_data <- rbind(combined_data, data_control)
  
  #elements in combined_data are in µg/g --> µmol/g
  combined_data$Fe_mol <- combined_data$Fe * combined_data$soilmass * 1000/ 55.845 /10^(6)
  combined_data$K_mol <- combined_data$K   * combined_data$soilmass * 1000/ 39.098/10^(6)
  combined_data$Ca_mol <- combined_data$Ca  * combined_data$soilmass * 1000/ 40.078/10^(6)
  combined_data$Mg_mol <- combined_data$Mg  * combined_data$soilmass * 1000/ 24.305/10^(6)
  combined_data$Na_mol <- combined_data$Na  * combined_data$soilmass * 1000/ 22.98976/10^(6)
  combined_data$Al_mol <- combined_data$Al  * combined_data$soilmass * 1000/ 26.98153/10^(6)
  combined_data$Si_mol <- combined_data$Si  * combined_data$soilmass * 1000/ 26.98153/10^(6)
    #units in mol element/mesocosm 
  data_allpools <- rbind(data_allpools,combined_data)
}
#replace <0 values (smaller than blank) by 0
data_allpools <- data_allpools %>%
  mutate(across(c(Fe_mol, Al_mol, Ca_mol, Mg_mol, K_mol, Na_mol, Si_mol), ~ ifelse(. < 0, 0, .)))
#OK, now normalize for controls
elements <- c("Fe","Al","Ca","Mg","Si")
plot_list <- list()

for (j in (1: length(elements))) {
  element <- elements[j]
  
  #subtract controls for each element
  library (Rmisc)
  library (ggplot2)
  
  measurevar <- paste0(element, "_mol")
  
  # Apply the summarySE function using the dynamically created measurevar
  sum <- summarySE(data_allpools[complete.cases(data_allpools[[measurevar]]), ], 
                   measurevar = measurevar, 
                   groupvars = c("Treatment", "pool"))
  
  
  # Initialize a list to store results
  results_list <- list()
  
  # Define a function to calculate difference and propagate standard error
  calculate_difference <- function(SIC1, SE1, SIC2, SE2) {
    diff <- SIC1 - SIC2
    propagated_se <- sqrt(SE1^2 + SE2^2)
    return(c(diff, propagated_se))
  }
  
  # Subtract values for FB1 by C1 for each pool
  for (pool in unique(sum$pool)) {
    results_list[[paste("FB1 - C1", pool)]] <- calculate_difference(
      sum[[measurevar]][sum$Treatment == "FB1" & sum$pool == pool], 
      sum$se[sum$Treatment == "FB1" & sum$pool == pool], 
      sum[[measurevar]][sum$Treatment == "C1" & sum$pool == pool], 
      sum$se[sum$Treatment == "C1" & sum$pool == pool]
    )
  }
  
  # Loop through treatments ending with '3' and subtract from C3, for each pool
  for (treatment in sum$Treatment[grepl("3$", sum$Treatment) & sum$Treatment != "C3"]) {
    for (pool in unique(sum$pool)) {
      comparison_name <- paste(treatment, "- C3", pool)
      results_list[[comparison_name]] <- calculate_difference(
        sum[[measurevar]][sum$Treatment == treatment & sum$pool == pool] , 
        sum$se[sum$Treatment == treatment & sum$pool == pool], 
        sum[[measurevar]][sum$Treatment == "C3" & sum$pool == pool], 
        sum$se[sum$Treatment == "C3" & sum$pool == pool]
      )
    }
  }
  
  # Initialize an empty list to store the final rows
  results_data <- list()
  
  # Loop through the names of the results_list to construct a dataframe
  for (name in names(results_list)) {
    # Split the comparison name into treatment and pool
    name_split <- unlist(strsplit(name, " "))
    
    # Extract the treatment name (first part before subtraction)
    treatment <- name_split[1]  # e.g., "FB1", "CB3", etc.
    
    # Extract the pool name (third element in the split string)
    name_spli2 <- unlist(strsplit(name, " "))
    
    pool <- name_split[4]  # e.g., "Exchangeable", "Carbonate", etc.
    
    # Extract the difference and propagated SE
    diff <- results_list[[name]][1]
    se <- results_list[[name]][2]
    
    # Append the row to the results_data list
    results_data[[name]] <- data.frame(
      diff = diff,
      se = se,
      treatment = treatment,
      pool = pool
    )
  }
  
  # Combine all rows into a single dataframe
  results_df <- do.call(rbind, results_data)
  
  
  
  library(dplyr)
  # Create a summary dataframe by treatment, summing the 'diff' and calculating propagated 'se'
  total_rows <- results_df %>%
    group_by(treatment) %>%
    summarize(
      diff = sum(diff, na.rm = TRUE),  # Sum the 'diff' for each treatment
      se = sqrt(sum(se^2, na.rm = TRUE))  # Calculate the propagated SE
    ) %>%
    mutate(pool = "total")  # Add a column for pool as "total"
  # Bind the original dataframe with the new 'total' rows
  results_df_with_total <- bind_rows(results_df, total_rows)
  
  
  #make a ggplot 
  # Assign the plot to the dynamically created name
  treatment_order <- c( "FB1", "CB3","FB3","CD3","MD3","FD3")
  # Convert Treatment column to factor with the specified order
  results_df_with_total$treatment <- factor(results_df_with_total$treatment, levels = treatment_order)
  
  
  results_df_with_total2 <- results_df_with_total #Store under different name
    results_df_with_total2$element <- element
empty_list[[j]] <- results_df_with_total2

  library (RColorBrewer)
  cols = brewer.pal(n = 8, name = 'Set2')

  #now add simulation data!!
  setwd ("C:/Users/AVienne/OneDrive - Universiteit Antwerpen/Labexp 3 (start Oct 2022)/PAPER_")
  treatmentvector <- c("FB1","FB3","CB3")
  #replace by for loop for i once all the data is there
 
  plot_name <- paste(element, "plot", sep = "")  # Create the plot name dynamically
  #make the combined plot now
  assign(plot_name, ggplot(results_df_with_total, aes(x = treatment, y = diff, color = pool)) +
           geom_point(size = 3,position = position_dodge(width = 0.5)) +  # Plot points for 'diff' values
           geom_errorbar(aes(ymin = diff - se, ymax = diff + se), width = 0.2, position = position_dodge(width = 0.5)) +  # Add error bars
           theme_bw() +  # Minimal theme for cleaner look
           theme(axis.title.y = element_text(size= 10), 
                 axis.title.x = element_text(size = 10), 
                 axis.title = element_text(face = "bold"),
                 panel.grid.major = element_blank(), 
                 panel.grid.minor = element_blank(), 
                 panel.background = element_blank(), 
                 axis.line = element_line(colour = "black")) +
           scale_color_manual(values = cols) +  # Set colors from cols
           labs(
             x = "Treatment", 
             y = paste("\u0394", "mol", element), 
             color = "Pool"
           ) +  theme(legend.position = "none") )
  
  plot_list[[j]]  <- get (plot_name)
  
  
  
  for (i in (1: length(treatmentvector))){  
    treatment <- treatmentvector[i]
    
    data <- read_excel(paste0("Simulationdata_",treatment,".xlsx"), sheet = "Sheet1", col_names = TRUE)
    #Data of day 130
    data130 <- data[which (data$days==386),]
    secondaryminerals <- c("gibbsite","goethite","calcite","Nesquehonite","CaMontmorillonite","Chlorite","Fe.OH.3.a.","Al.OH.3.a.","SiO2.a.")
    Alvector <- c(1,0,0,0,2.33,2,0,1,0)
    Cavector <- c(0,0,1,0,0.165,0,0,0,0)
    Mgvector <- c(0,0,0,1,0,5,0,0,0)
    Fevector <- c(0,1,0,0,0,0,1,0,0)
    Sivector <- c(0,0,0,0,3.67,3,0,0,1) #mol/mol
    
    for (k in (1:length (secondaryminerals))) {
      mineral <- secondaryminerals[k]
      diff <- data130[data130$type == "rock", mineral] - data130[data130$type == "control", mineral]
      #mol difference after 130 days
      
      Alvector[k] <- diff * Alvector[k]
      Cavector[k] <- diff * Cavector[k]
      Mgvector[k] <- diff * Mgvector[k]
      Fevector[k] <- diff * Fevector[k]
      Sivector[k] <- diff * Sivector[k] }
    
    results_df <- data.frame(Mineral = secondaryminerals,
                             Al = unlist(Alvector),
                             Ca = unlist(Cavector),
                             Mg = unlist(Mgvector),
                             Fe = unlist(Fevector),
                             Si = unlist(Sivector), stringsAsFactors = FALSE , treatment)
    
    if (!(element %in% c("Na", "K"))) {
      results_df_element <- results_df[, c(element, "Mineral","treatment"), drop = FALSE]
      results_df_element <-results_df_element[-which(results_df_element[,1]==0),]
      #calculate total simulated value
      total_value <- sum(results_df_element[[element]], na.rm = TRUE)
      total_row <- data.frame(
        element = total_value,      # Use dynamic column name
        Mineral = "total",             # Set the mineral as "total"
        treatment = unique(results_df_element$treatment))
      colnames (total_row)[1] <- element
      
      # Get the treatment from the existing data
      results_df_element <- rbind(results_df_element, total_row)
      available_shapes <- c(17, 18, 15, 19, 3, 4, 7, 20, 9)
      shape_mapping <- setNames(available_shapes[1:length(secondaryminerals)], secondaryminerals)
      shape_mapping["total"] <- 8
    }
    
    secondaryminerals <- c("gibbsite","goethite","calcite","Nesquehonite","CaMontmorillonite","Chlorite","Fe.OH.3.a.","Al.OH.3.a.","SiO2.a.")
    
    
    
    #trick for legend below:
    secondaryminerals <- c("gibbsite","goethite","calcite","Nesquehonite","CaMontmorillonite","Chlorite","Fe.OH.3.a.","Al.OH.3.a.","SiO2.a.")
    mineral_assignments <- rep(secondaryminerals, length.out = nrow(results_df_with_total))
    # Shuffle the mineral_assignments to create random order
    set.seed(123)  # Set seed for reproducibility
    mineral_assignments <- sample(mineral_assignments)
    mineral_assignments[1]<- "total"
    # Add the Minerals column to the dataframe
    results_df_with_total$Mineral <- mineral_assignments
    
    name <- "forlegend"
    available_shapes <- c(17, 18, 15, 19, 3, 4, 7, 20, 9)
    shape_mapping <- setNames(available_shapes[1:length(secondaryminerals)], secondaryminerals)
    shape_mapping["total"] <- 8
    
    assign(name, ggplot(results_df_with_total, aes(x = treatment, y = diff, color = pool)) +
             geom_point(size = 3,position = position_dodge(width = 0.5)) +  # Plot points for 'diff' values
             geom_errorbar(aes(ymin = diff - se, ymax = diff + se), width = 0.2, position = position_dodge(width = 0.5)) +  # Add error bars
             theme_bw() +  # Minimal theme for cleaner look
             theme(axis.title.y = element_text(size= 10), 
                   axis.title.x = element_text(size = 10), 
                   axis.title = element_text(face = "bold"),
                   panel.grid.major = element_blank(), 
                   panel.grid.minor = element_blank(), 
                   panel.background = element_blank(), 
                   axis.line = element_line(colour = "black")) +
             scale_color_manual(values = cols, name = "Pool \n (extracted)") +  # Set colors from cols
             labs(
               x = "Treatment", 
               y = paste("\u0394", "mol", element), 
               color = "Pool"
             )       + 
             geom_point(data = results_df_with_total, 
                        aes(x = treatment, y = diff, shape = Mineral),  # Correct dynamic column access
                        size = 3, fill = NA, col = "black") + 
             scale_shape_manual(values = shape_mapping, name = "Secondary minerals \n (simulated)"))
   
        # Access the plot from the list by index
      # Modify the plot by adding points
      plot_list[[j]] <- plot_list[[j]] + 
        geom_point(data = results_df_element, 
                   aes(x = treatment, y = .data[[element]], shape = Mineral),  # Correct dynamic column access
                   size = 3, fill = NA, col = "black") + 
        scale_shape_manual(values = shape_mapping)
  }
}

  dataframe_2 <- do.call(rbind, empty_list)



library (cowplot)
legend <- get_legend (forlegend)
final_plot <- plot_grid (plot_grid(plotlist = plot_list, labels = "AUTO"), legend, ncol = 2,  rel_widths  = c(1, 0.2))

output_path <- "C:/Users/AVienne/OneDrive - Universiteit Antwerpen/Labexp 3 (start Oct 2022)/PAPER_/Sec_minerals_heckman2.jpeg"
ggsave(filename = output_path, plot = final_plot, width = 10, height = 8, units = "in", dpi = 600)

```

dataframe_1$day <- 130
dataframe_2$day <- 386
dataframe_0 <- dataframe_1
dataframe_0$day <- 0
dataframe_0$diff <- 0
dataframe_0$se <- 0



# Combine the dataframes
combined_data <- rbind(dataframe_0, dataframe_1, dataframe_2)
library(ggplot2)
# for every element 
plotlist <- list ()

for (i in (1: length (elements))) {
element <- elements[i]
combined_datatemp <- combined_data[which(combined_data$element == element & combined_data$pool == "total" ), ]
plotlist [[i]] <- ggplot(combined_datatemp, aes(x = day, y = diff, color = treatment, fill = treatment, shape = treatment, size = treatment)) +
  geom_point() + 
  geom_errorbar(aes(ymin = diff - se, ymax = diff + se), width = 0.02, color = "black", size = 0.3) +
  labs(x = "Time (Days)", y = paste("\u0394", "mol", element)) +
  geom_line(linetype = "dotted", size = 0.5, color = "black") +
  scale_shape_manual(values = c("FB1" = 21, "CB3" = 21, "FB3" = 21)) + 
  scale_size_manual(values = c("FB1" = 2, "CB3" = 4, "FB3" = 2)) + 
  scale_color_manual(values = c("FB1" = "grey", "CB3" = "grey", "FB3" = "grey")) +
    scale_fill_manual(values = c("FB1" = "white", "CB3" = "grey",  "FB3" = "grey"))+
  theme_bw() + theme(axis.title.y = element_text(size= 10), axis.title.x = element_text(size= 10), axis.title = element_text(face = "bold")) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour= "black"))+
  theme(axis.text.x = element_text(size = 10),  # Increase size of x-axis ticks
    axis.text.y = element_text(size = 10)) +  theme(legend.position = "none")
}


forlegend <-  ggplot(combined_datatemp, aes(x = day, y = diff, color = treatment, fill = treatment, shape = treatment, size = treatment)) +
  geom_point() + 
  geom_errorbar(aes(ymin = diff - se, ymax = diff + se), width = 0.02, color = "black", size = 0.3) +
  labs(x = "Time (Days)", y = paste("\u0394", "mol", element)) +
  geom_line(linetype = "dotted", size = 0.5, color = "black") +
  scale_shape_manual(values = c("FB1" = 21, "CB3" = 21, "FB3" = 21)) + 
  scale_size_manual(values = c("FB1" = 2, "CB3" = 4, "FB3" = 2)) + 
  scale_color_manual(values = c("FB1" = "grey", "CB3" = "grey", "FB3" = "grey")) +
    scale_fill_manual(values = c("FB1" = "white", "CB3" = "grey",  "FB3" = "grey"))+
  theme_bw() + theme(axis.title.y = element_text(size= 15), axis.title.x = element_text(size= 20), axis.title = element_text(face = "bold")) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour= "black"))+
  theme(axis.text.x = element_text(size = 6),  # Increase size of x-axis ticks
    axis.text.y = element_text(size = 6))

legend <- get_legend (forlegend)
final_plot <- plot_grid (plot_grid(plotlist = plotlist, labels = "AUTO"), legend, ncol = 2,  rel_widths  = c(1, 0.2))

output_path <- "C:/Users/AVienne/OneDrive - Universiteit Antwerpen/Labexp 3 (start Oct 2022)/PAPER_/Sec_minerals_heckman_timeplot_totals_allelements.jpeg"
ggsave(filename = output_path, plot = final_plot, width = 10, height = 8, units = "in", dpi = 600)

#conclusion: elemental uptake of Ca, Fe and Mg > weathering of those elements. Balance shifts towards formation of unextractable secondary minerals.
#Only extracted Si and Al still increased during the second half of the year.

#stats for timeseries



