--- a/insect_phenology_model.R	Sun Nov 12 10:18:59 2017 -0500
+++ b/insect_phenology_model.R	Sun Nov 12 10:19:10 2017 -0500
@@ -25,10 +25,13 @@
 args <- parse_args(parser, positional_arguments=TRUE)
 opt <- args$options
 
-convert_csv_to_rdata=function(temperature_data, data_matrix)
+convert_csv_to_rdata=function(start_doy, end_doy, temperature_data, data_matrix)
 {
-    # Integer day of the year.
-    data_matrix[,1] <- c(1:opt$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
@@ -38,12 +41,12 @@
     namedat
 }
 
-daylength=function(latitude, num_days)
+daylength=function(latitude, start_doy, end_doy)
 {
     # From Forsythe 1995.
     p=0.8333
     dl <- NULL
-    for (i in 1:num_days) {
+    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)))
@@ -70,7 +73,7 @@
     else {
         # Extract daylength data for the number of
         # days specified in the input temperature data.
-        dlprofile <- daylength(latitude, num_days)
+        dlprofile <- daylength(latitude, start_doy, end_doy)
         # Initialize hourly temperature.
         T <- NULL
         # Initialize degree hour vector.
@@ -199,13 +202,19 @@
 }
 
 # Read in the input temperature datafile into a Data Frame object.
-temperature_data <- read.csv(file=opt$input, header=T, sep=",")
+# 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[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(temperature_data, raw_data_matrix)
+temperature_file_path <- convert_csv_to_rdata(start_doy, end_doy, temperature_data, raw_data_matrix)
 latitude <- temperature_data[1, 1]
 
+cat("Start date: ", start_date, "\n")
+cat("End date: ", end_date, "\n")
 cat("Number of days: ", opt$num_days, "\n")
 
 # Initialize matrix for results from all replications.
@@ -223,8 +232,9 @@
     vec.mat <- rep(vec.ini, n)
     # Complete matrix for the population.
     vec.mat <- base::t(matrix(vec.mat, nrow=5))
-    # Complete photoperiod profile in a year, requires daylength function.
-    ph.p <- daylength(latitude, opt$num_days)
+    # 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
@@ -237,7 +247,7 @@
     dd.day <- rep(0, opt$num_days)
 
     # All the days included in the input temperature dataset.
-    for (day in 1:opt$num_days) {
+    for (day in start_doy:end_doy) {
         # Photoperiod in the day.
         photoperiod <- ph.p[day]
         temp.profile <- hourtemp(latitude, day, temperature_file_path, opt$num_days)
@@ -534,11 +544,12 @@
     g2a.rep[,N.rep] <- g2.adult
 }
 
-# Data analysis and visualization
-# default: plot 1 year of result
-# but can be expanded to accommodate multiple years
-n.yr <- 1
-day.all <- c(1:opt$num_days * n.yr)
+# 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)
 
 # mean value for adults
 sa <- apply((S3.rep + S4.rep + S5.rep), 1, mean)