--- a/insect_phenology_model.R	Tue Nov 14 11:35:02 2017 -0500
+++ b/insect_phenology_model.R	Tue Nov 14 11:35:11 2017 -0500
@@ -11,6 +11,7 @@
     make_option(c("-i", "--max_clutch_size"), action="store", dest="max_clutch_size", type="integer", help="Adjustment of maximum clutch size"),
     make_option(c("-j", "--nymph_mort"), action="store", dest="nymph_mort", type="integer", help="Adjustment rate for nymph mortality"),
     make_option(c("-k", "--old_nymph_accum"), action="store", dest="old_nymph_accum", type="integer", help="Adjustment of DD accumulation (young nymph->old nymph)"),
+    make_option(c("-m", "--num_columns"), action="store", dest="num_columns", type="integer", help="Total number of columns in the temperature dataset"),
     make_option(c("-n", "--num_days"), action="store", dest="num_days", type="integer", help="Total number of days in the temperature dataset"),
     make_option(c("-o", "--output"), action="store", dest="output", help="Output dataset"),
     make_option(c("-p", "--oviposition"), action="store", dest="oviposition", type="integer", help="Adjustment for oviposition rate"),
@@ -25,8 +26,7 @@
 args <- parse_args(parser, positional_arguments=TRUE)
 opt <- args$options
 
-get_daylight_length = function(latitude, temperature_data, num_days)
-{
+add_daylight_length = function(latitude, temperature_data_matrix, num_days) {
     # Return a vector of daylight length (photoperido profile) for
     # the number of days specified in the input temperature data
     # (from Forsythe 1995).
@@ -35,17 +35,17 @@
     for (i in 1:num_days) {
         # Get the day of the year from the current row
         # of the temperature data for computation.
-        doy <- temperature_data[i, 4]
+        doy <- temperature_data_matrix[i, 4]
         theta <- 0.2163108 + 2 * atan(0.9671396 * tan(0.00860 * (doy - 186)))
         phi <- asin(0.39795 * cos(theta))
         # Compute the length of daylight for the day of the year.
         daylight_length_vector[i] <- 24 - (24 / pi * acos((sin(p * pi / 180) + sin(latitude * pi / 180) * sin(phi)) / (cos(latitude * pi / 180) * cos(phi))))
     }
-    daylight_length_vector
+    # Append daylight_length_vector as a new column to temperature_data_matrix.
+    temperature_data_matrix <- cbind(temperature_data_matrix, daylight_length_vector)
 }
 
-get_temperature_at_hour = function(latitude, temperature_data, daylight_length_vector, row, num_days)
-{
+get_temperature_at_hour = function(latitude, temperature_data_matrix, row, num_days) {
     # Base development threshold for Brown Marmolated Stink Bug
     # insect phenology model.
     # TODO: Pass insect on the command line to accomodate more
@@ -55,9 +55,9 @@
     # Input temperature currently has the following columns.
     # # LATITUDE, LONGITUDE, DATE, DOY, TMIN, TMAX
     # Minimum temperature for current row.
-    dnp <- temperature_data[row, 5]
+    dnp <- temperature_data_matrix[row, 5]
     # Maximum temperature for current row.
-    dxp <- temperature_data[row, 6]
+    dxp <- temperature_data_matrix[row, 6]
     # Mean temperature for current row.
     dmean <- 0.5 * (dnp + dxp)
     # Initialize degree day accumulation
@@ -71,7 +71,7 @@
         # Initialize degree hour vector.
         dh <- NULL
         # Daylight length for current row.
-        y <- daylight_length_vector[row]
+        y <- temperature_data_matrix[row, 7]
         # Darkness length.
         z <- 24 - y
         # Lag coefficient.
@@ -122,15 +122,13 @@
     return
 }
 
-dev.egg = function(temperature)
-{
+dev.egg = function(temperature) {
     dev.rate= -0.9843 * temperature + 33.438
     return = dev.rate
     return
 }
 
-dev.young = function(temperature)
-{
+dev.young = function(temperature) {
     n12 <- -0.3728 * temperature + 14.68
     n23 <- -0.6119 * temperature + 25.249
     dev.rate = mean(n12 + n23)
@@ -138,8 +136,7 @@
     return
 }
 
-dev.old = function(temperature)
-{
+dev.old = function(temperature) {
     n34 <- -0.6119 * temperature + 17.602
     n45 <- -0.4408 * temperature + 19.036
     dev.rate = mean(n34 + n45)
@@ -147,15 +144,13 @@
     return
 }
 
-dev.emerg = function(temperature)
-{
+dev.emerg = function(temperature) {
     emerg.rate <- -0.5332 * temperature + 24.147
     return = emerg.rate
     return
 }
 
-mortality.egg = function(temperature)
-{
+mortality.egg = function(temperature) {
     if (temperature < 12.7) {
         mort.prob = 0.8
     }
@@ -169,8 +164,7 @@
     return
 }
 
-mortality.nymph = function(temperature)
-{
+mortality.nymph = function(temperature) {
     if (temperature < 12.7) {
         mort.prob = 0.03
     }
@@ -181,8 +175,7 @@
     return
 }
 
-mortality.adult = function(temperature)
-{
+mortality.adult = function(temperature) {
     if (temperature < 12.7) {
         mort.prob = 0.002
     }
@@ -196,11 +189,11 @@
 # Read in the input temperature datafile into a Data Frame object.
 # The input data currently must have 6 columns:
 # LATITUDE, LONGITUDE, DATE, DOY, TMIN, TMAX
-temperature_data <- read.csv(file=opt$input, header=T, strip.white=TRUE, sep=",")
-start_date <- temperature_data[1, 3]
-end_date <- temperature_data[opt$num_days, 3]
-latitude <- temperature_data[1, 1]
-daylight_length_vector <- get_daylight_length(latitude, temperature_data, opt$num_days)
+temperature_data_matrix <- read.csv(file=opt$input, header=T, strip.white=TRUE, sep=",")
+start_date <- temperature_data_matrix[1, 3]
+end_date <- temperature_data_matrix[opt$num_days, 3]
+latitude <- temperature_data_matrix[1, 1]
+add_daylight_length(latitude, temperature_data_matrix, opt$num_days)
 
 cat("Number of days: ", opt$num_days, "\n")
 
@@ -232,10 +225,10 @@
     # All the days included in the input temperature dataset.
     for (row in 1:opt$num_days) {
         # Get the integer day of the year for the current row.
-        doy <- temperature_data[row, 4]
+        doy <- temperature_data_matrix[row, 4]
         # Photoperiod in the day.
-        photoperiod <- daylight_length_vector[row]
-        temp.profile <- get_temperature_at_hour(latitude, temperature_data, daylight_length_vector, row, opt$num_days)
+        photoperiod <- temperature_data_matrix[row, 7]
+        temp.profile <- get_temperature_at_hour(latitude, temperature_data_matrix, row, opt$num_days)
         mean.temp <- temp.profile[1]
         dd.temp <- temp.profile[2]
         dd.day[row] <- dd.temp