Random Walk Particle Transport in Porous Heterogeneous Media

library(readr)
library(purrr)
library(dplyr)

dat <- read_csv("part_evol.csv")
part_df <- read.csv("part_df_tort.csv")

library(cowplot)
library(ggplot2)
part_df %>% 
  ggplot() +
  geom_histogram(aes(x = tort), fill = "grey", color = "black", binwidth = 2, center = 1) +
  #geom_freqpoly(aes(x = tort), binwidth = 3, color = "red") +
  theme_bw() +
  theme_cowplot() -> tort_hist
#tort_hist %>% ggsave(file = "tort_hist.png")

# get this data in a table
b <- ggplot_build(tort_hist)
b$data %>% flatten_df() -> b_dat
b_dat %>% select(xmin, xmax, density, ymax) %>% rename(frequency = ymax) -> tort_table

tort_table %>% 
  mutate(cum_frequency = cumsum(density)*2) %>% 
  ggplot(aes(x=xmax, y=cum_frequency)) +  
  geom_step() + xlab("Tortuosity") + ylab("Cumulative Density") -> tort_cum_density

#tort_cum_density %>% ggsave(file = "tort_density.png")

library(rowr)
library(plotly)

plot_3d_line <- function(df, dat, N){ # df = tortousity summary of subset dat, dat = all data, N = max lines to plot
  temp <- dat %>% filter(id %in% df$id) # get all particles in df from dat
  num_part <- temp$id %>% unique() %>% length() # number of times to repeat column names
  temp2 <- split(temp, temp$id)  # split by particle id
  
  # combine into df
  part_df <- do.call(cbind.fill, c(temp2, fill = NA)) # slow.
  
  # add new column names for iterating traces in a loop
  id <- df$id %>% unique() # grab particle ids
  id <- rep(id, each = 6) 
  prefix <- rep(c("t_","x_","y_","z_","sp_","id_"), times = num_part)
  nam <- paste0(prefix, id)
  
  # assign colnames
  colnames(part_df) <- nam

  # add time
  part_df$t <- temp$t %>% unique()
  
  
  # Plot
  temp %>% 
    group_by(id) %>% 
    summarise(tsteps = n()) %>% 
    arrange(desc(tsteps)) %>% 
    top_n(1000, wt = tsteps) -> ranked_par

    p <- plot_ly(type = 'scatter3d', mode = 'lines')
      for(i in sample(ranked_par$id, N)){
        x = paste0("x_",i)
        y = paste0("y_",i)
        z = paste0("z_",i)
        t = paste0("t_",i)
        p <- add_trace(p, 
                       x = part_df[[x]], 
                       y = part_df[[y]], 
                       z = part_df[[z]],
                       #frame = part_df[[t]],
                       type = 'scatter3d', 
                       mode = 'lines',
                       line = list(width = 1))
        }
    p <- p %>% layout(showlegend = FALSE,
                      scene = list(xaxis = list(title = 'x',
                                   range = c(min(dat$x), max(dat$x))),
                      yaxis = list(title = 'y',
                                   range = c(min(dat$y), max(dat$y))),
                      zaxis = list(title = 'z ',
                                   range = c(min(dat$z), max(dat$z)))
                      ))
    return(p)
}


# check out bins
temp <- part_df %>% 
        filter(tort >= 1 & tort <=2) %>% 
        sample_n(100, replace = FALSE)

plot_3d_line(temp, dat, 100)