--- a/insect_phenology_model.R	Tue Nov 14 13:23:36 2017 -0500
+++ b/insect_phenology_model.R	Tue Nov 14 13:54:27 2017 -0500
@@ -26,7 +26,19 @@
 args <- parse_args(parser, positional_arguments=TRUE)
 opt <- args$options
 
-add_daylight_length = function(latitude, temperature_data_matrix, num_days) {
+parse_input_data = function(input_file, num_rows, num_columns) {
+    # Read in the input temperature datafile into a data frame.
+    if (num_rows == 6) {
+        # The input data has the following 6 columns:
+        # LATITUDE, LONGITUDE, DATE, DOY, TMIN, TMAX
+        temperature_data_frame <- read.csv(file=input_file, header=T, strip.white=TRUE, sep=",")
+        # Add a column containing the daylight length for each day.
+        temperature_data_frame <- add_daylight_length(latitude, temperature_data_frame, opt$num_days)
+        temperature_data_frame
+    }
+}
+
+add_daylight_length = function(latitude, temperature_data_frame, 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,19 +47,19 @@
     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_matrix[i, 4]
+        doy <- temperature_data_frame[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))))
     }
-    # Append daylight_length_vector as a new column to temperature_data_matrix.
-    temperature_data_matrix[, 7] <- daylight_length_vector
-    # Return the altered temperature_data_matrix.
-    temperature_data_matrix
+    # Append daylight_length_vector as a new column to temperature_data_frame.
+    temperature_data_frame[, 7] <- daylight_length_vector
+    # Return the altered temperature_data_frame.
+    temperature_data_frame
 }
 
-get_temperature_at_hour = function(latitude, temperature_data_matrix, row, num_days) {
+get_temperature_at_hour = function(latitude, temperature_data_frame, row, num_days) {
     # Base development threshold for Brown Marmolated Stink Bug
     # insect phenology model.
     # TODO: Pass insect on the command line to accomodate more
@@ -57,9 +69,9 @@
     # Input temperature currently has the following columns.
     # # LATITUDE, LONGITUDE, DATE, DOY, TMIN, TMAX
     # Minimum temperature for current row.
-    dnp <- temperature_data_matrix[row, 5]
+    dnp <- temperature_data_frame[row, 5]
     # Maximum temperature for current row.
-    dxp <- temperature_data_matrix[row, 6]
+    dxp <- temperature_data_frame[row, 6]
     # Mean temperature for current row.
     dmean <- 0.5 * (dnp + dxp)
     # Initialize degree day accumulation
@@ -73,7 +85,7 @@
         # Initialize degree hour vector.
         dh <- NULL
         # Daylight length for current row.
-        y <- temperature_data_matrix[row, 7]
+        y <- temperature_data_frame[row, 7]
         # Darkness length.
         z <- 24 - y
         # Lag coefficient.
@@ -188,17 +200,10 @@
     return
 }
 
-# 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_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]
-temperature_data_matrix <- add_daylight_length(latitude, temperature_data_matrix, opt$num_days)
+temperature_data_frame <- parse_input_data(opt$input, opt$num_rows, opt$num_columns)
+latitude <- temperature_data_frame[1, 1]
 
 cat("Number of days: ", opt$num_days, "\n")
-print(temperature_data_matrix)
 
 # Initialize matrix for results from all replications.
 S0.rep <- S1.rep <- S2.rep <- S3.rep <- S4.rep <- S5.rep <- matrix(rep(0, opt$num_days * opt$replications), ncol = opt$replications)
@@ -228,10 +233,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_matrix[row, 4]
+        doy <- temperature_data_frame[row, 4]
         # Photoperiod in the day.
-        photoperiod <- temperature_data_matrix[row, 7]
-        temp.profile <- get_temperature_at_hour(latitude, temperature_data_matrix, row, opt$num_days)
+        photoperiod <- temperature_data_frame[row, 7]
+        temp.profile <- get_temperature_at_hour(latitude, temperature_data_frame, row, opt$num_days)
         mean.temp <- temp.profile[1]
         dd.temp <- temp.profile[2]
         dd.day[row] <- dd.temp
@@ -532,6 +537,9 @@
 # Data analysis and visualization can currently
 # plot only within a single calendar year.
 # TODO: enhance this to accomodate multiple calendar years.
+start_date <- temperature_data_frame[1, 3]
+end_date <- temperature_data_frame[opt$num_days, 3]
+
 n.yr <- 1
 day.all <- c(1:opt$num_days * n.yr)