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")
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)
}
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")
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)