--- a/insect_phenology_model.R	Sun Nov 12 10:38:18 2017 -0500
+++ b/insect_phenology_model.R	Mon Nov 13 10:42:31 2017 -0500
@@ -25,66 +25,57 @@
 args <- parse_args(parser, positional_arguments=TRUE)
 opt <- args$options
 
-convert_csv_to_rdata=function(start_doy, end_doy, temperature_data, data_matrix)
+get_daylight_length = function(latitude, temperature_data, num_days)
 {
-    # Make sure the first column starts with the start date in
-    # the raw csv data converted to the integer day of the year
-    # and continues through the end date in the raw csv data
-    # converted to the integer day of the year.
-    data_matrix[,1] <- c(start_doy:end_doy)
-    # Minimum
-    data_matrix[,2] <- temperature_data[c(1:opt$num_days), 5]
-    # Maximum
-    data_matrix[,3] <- temperature_data[c(1:opt$num_days), 6]
-    namedat <- "tempdata.Rdat"
-    save(data_matrix, file=namedat)
-    namedat
+    # Return a vector of daylight length (photoperido profile) for
+    # the number of days specified in the input temperature data
+    # (from Forsythe 1995).
+    p = 0.8333
+    daylight_length_vector <- NULL
+    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[c(i: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
 }
 
-daylength=function(latitude, start_doy, end_doy)
+get_temperature_at_hour = function(latitude, temperature_data, daylight_length_vector, row, num_days)
 {
-    # From Forsythe 1995.
-    p=0.8333
-    dl <- NULL
-    for (i in start_doy:end_doy) {
-        theta <- 0.2163108 + 2 * atan(0.9671396 * tan(0.00860 * (i - 186)))
-        phi <- asin(0.39795 * cos(theta))
-        dl[i] <- 24 - 24 / pi * acos((sin(p * pi / 180) + sin(latitude * pi / 180) * sin(phi)) / (cos(latitude * pi / 180) * cos(phi)))
-    }
-    # Return a vector of daylength for the number of
-    # days specified in the input temperature data.
-    dl
-}
-
-hourtemp=function(latitude, date, temperature_file_path, num_days)
-{
-    load(temperature_file_path)
     # Base development threshold for Brown Marmolated Stink Bug
     # insect phenology model.
+    # TODO: Pass insect on the command line to accomodate more
+    # the just the Brown Marmolated Stink Bub.
     threshold <- 14.17
-    dnp <- data_matrix[date, 2]  # daily minimum
-    dxp <- data_matrix[date, 3]  # daily maximum
+
+    # Input temperature currently has the following columns.
+    # # LATITUDE, LONGITUDE, DATE, DOY, TMIN, TMAX
+    # Minimum temperature for current row.
+    dnp <- temperature_data[row, 5]
+    # Maximum temperature for current row.
+    dxp <- temperature_data[row, 6]
+    # Mean temperature for current row.
     dmean <- 0.5 * (dnp + dxp)
     dd <- 0  # initialize degree day accumulation
-
-    if (dxp<threshold) {
+    if (dxp < threshold) {
         dd <- 0
     }
     else {
-        # Extract daylength data for the number of
-        # days specified in the input temperature data.
-        dlprofile <- daylength(latitude, start_doy, end_doy)
         # Initialize hourly temperature.
         T <- NULL
         # Initialize degree hour vector.
         dh <- NULL
-        # Calculate daylength in given date.
-        y <- dlprofile[date]
-        # Night length.
+        # Daylight length for current row.
+        y <- daylight_length_vector[row]
+        # Darkness length.
         z <- 24 - y
         # Lag coefficient.
         a <- 1.86
-        # Night coefficient.
+        # Darkness coefficient.
         b <- 2.20
         # Sunrise time.
         risetime <- 12 - y / 2
@@ -92,11 +83,11 @@
         settime <- 12 + y / 2
         ts <- (dxp - dnp) * sin(pi * (settime - 5) / (y + 2 * a)) + dnp
         for (i in 1:24) {
-            if (i > risetime && i<settime) {
+            if (i > risetime && i < settime) {
                 # Number of hours after Tmin until sunset.
                 m <- i - 5
-                T[i]=(dxp - dnp) * sin(pi * m / (y + 2 * a)) + dnp
-                if (T[i]<8.4) {
+                T[i] = (dxp - dnp) * sin(pi * m / (y + 2 * a)) + dnp
+                if (T[i] < 8.4) {
                     dh[i] <- 0
                 }
                 else {
@@ -105,8 +96,8 @@
             }
             else if (i > settime) { 
                 n <- i - settime
-                T[i]=dnp + (ts - dnp) * exp( - b * n / z)
-                if (T[i]<8.4) {
+                T[i] = dnp + (ts - dnp) * exp( - b * n / z)
+                if (T[i] < 8.4) {
                     dh[i] <- 0
                 }
                 else {
@@ -116,7 +107,7 @@
             else {
                 n <- i + 24 - settime
                 T[i]=dnp + (ts - dnp) * exp( - b * n / z)
-                if (T[i]<8.4) {
+                if (T[i] < 8.4) {
                     dh[i] <- 0
                 }
                 else {
@@ -206,12 +197,9 @@
 # LATITUDE, LONGITUDE, DATE, DOY, TMIN, TMAX
 temperature_data <- read.csv(file=opt$input, header=T, strip.white=TRUE, sep=",")
 start_date <- temperature_data[c(1:1), 3]
-start_doy <- temperature_data[c(1:1), 4]
 end_date <- temperature_data[c(opt$num_days:opt$num_days), 3]
-end_doy <- temperature_data[c(opt$num_days:opt$num_days), 4]
-raw_data_matrix <- matrix(rep(0, opt$num_days * 6), nrow=opt$num_days)
-temperature_file_path <- convert_csv_to_rdata(start_doy, end_doy, temperature_data, raw_data_matrix)
 latitude <- temperature_data[1, 1]
+daylight_length_vector <- get_daylight_length(latitude, temperature_data, opt$num_days)
 
 cat("Start date: ", start_date, "\n")
 cat("End date: ", end_date, "\n")
@@ -232,10 +220,6 @@
     vec.mat <- rep(vec.ini, n)
     # Complete matrix for the population.
     vec.mat <- base::t(matrix(vec.mat, nrow=5))
-    # Complete photoperiod profile for the total
-    # number of days in the input data.
-    ph.p <- daylength(latitude, start_doy, end_doy)
-
     # Time series of population size.
     tot.pop <- NULL
     gen0.pop <- rep(0, opt$num_days)
@@ -247,10 +231,12 @@
     dd.day <- rep(0, opt$num_days)
 
     # All the days included in the input temperature dataset.
-    for (day in start_doy:end_doy) {
+    for (row in 1:opt$num_days) {
+        # Get the integer day of the year for the current row.
+        doy <- temperature_data[c(row:row), 4]
         # Photoperiod in the day.
-        photoperiod <- ph.p[day]
-        temp.profile <- hourtemp(latitude, day, temperature_file_path, opt$num_days)
+        photoperiod <- daylight_length_vector[row]
+        temp.profile <- get_temperature_at_hour(latitude, temperature_data, daylight_length_vector, row, opt$num_days)
         mean.temp <- temp.profile[1]
         dd.temp <- temp.profile[2]
         dd.day[day] <- dd.temp
@@ -545,11 +531,10 @@
 }
 
 # Data analysis and visualization can currently
-# plot only within a single calendar year.  TODO:
-# enhance this to accomodate multiple calendar years.
-#n.yr <- 1
-#day.all <- c(1:opt$num_days * n.yr)
-day.all <- c(start_doy:end_doy)
+# plot only within a single calendar year.
+# TODO: enhance this to accomodate multiple calendar years.
+n.yr <- 1
+day.all <- c(1:opt$num_days * n.yr)
 
 # mean value for adults
 sa <- apply((S3.rep + S4.rep + S5.rep), 1, mean)