Just getting some initializing stuff out of the way and bring in some outside data that was generated in thermoclineLevels.r and solarData.r:

setwd("/Users/Nick/mysisModeling")
library(ggplot2)
depthData = read.csv("data/Depth_Thermocline_Hour.csv")$dist
lightData = read.csv("data/light_day_moon_hour.csv")$x
hours = 1:(365*24)

Do a quick plot to make sure the data came in properly:

lightDataDf = as.data.frame(cbind(hours, lightData))
m = ggplot(lightDataDf, aes(x = hours, y = lightData)) + geom_line() 
m + labs(title = "Light intensity levels by hour", x = "hour", y = "light intensity")

Mylux limit model:

Now we code put in the model for depth at mylux intensity limit.

Distance of light threshold: \(f(I_o) = \frac{1}{k}(ln(I_o) - ln(I_x))\)

lightDepth = function(surfaceLight){
  
  k   = 0.3  #extinction coefficient
  I_x = 0.001 #Mysis light threshold (paper quotes between 10^-2 and 10^-4)
  
  distance = (1/k) * (log(surfaceLight) - log(I_x))
  if (distance < 0){ distance = 0 }
  return(distance)
}

Mylux limit depth

Now we will run the raw data through our light depth function and plot:

isocline = NULL
for (hour in lightData){
  if (hour == 0){
    isocline = c(isocline, 0) #logs dont play nice with 0s
  } else {
    isocline = c(isocline, lightDepth(hour))
  }
}

isoclineDf = as.data.frame(cbind(hours, isocline))
m = ggplot(isoclineDf, aes(x = hours, y = -isocline)) + geom_line()
m + labs(title = "Depth of Mylux Limit", x = "hour", y = "depth from surface")

Mysocline

Combining the isocline with the thermocline to get the mysocline and output to file:

mysocline = NULL
for (i in 1:(24*365)){
  if (depthData[i] > isocline[i]){
    mysocline = c(mysocline, depthData[i])
  } else {
    mysocline = c(mysocline, isocline[i])
  }
}

mysoclineDf = as.data.frame(cbind(hours, mysocline))
m = ggplot(mysoclineDf[1:24,], aes(x = hours, y = -mysocline)) + geom_line()
m + labs(title = "Depth of Mysocline", x = "hour", y = "depth from surface")

#write.csv(mysocline, "data/mysocline_hour.csv", row.names=FALSE)