--- a/bmsb.R	Fri Dec 16 10:16:11 2016 -0500
+++ b/bmsb.R	Sun Dec 18 13:28:04 2016 -0500
@@ -26,20 +26,22 @@
 args <- parse_args(parser, positional_arguments=TRUE)
 opt <- args$options
 
-data.input=function(loc, start.yr, temperature.dataset)
+data.input=function(loc, year, temperature.dataset)
 {
     expdata <- matrix(rep(0, 365 * 3), nrow=365)
-    # replace 2004 with start. yr
-    yr <- start.yr
-    namedat <- paste(loc,  yr, ".Rdat", sep="")
+    namedat <- paste(loc,  year, ".Rdat", sep="")
     temp.data <- read.csv(file=temperature.dataset, header=T)
 
     expdata[,1] <- c(1:365)
+    # Minimum
+    expdata[,2] <- temp.data[c(1:365), 3]
+    # Maximum
+    expdata[,3] <- temp.data[c(1:365), 2]
     save(expdata, file=namedat)
     namedat
 }
 
-daylength=function(L)
+daylength=function(latitude)
 {
     # from Forsythe 1995
     p=0.8333
@@ -48,12 +50,12 @@
     {
         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(L * pi / 180) * sin(phi)) / (cos(L * pi / 180) * cos(phi)))
+        dl[i] <- 24 - 24 / pi * acos((sin(p * pi / 180) + sin(latitude * pi / 180) * sin(phi)) / (cos(latitude * pi / 180) * cos(phi)))
     }
     dl   # return a vector of daylength in 365 days
 }
 
-hourtemp=function(L, date, temperature_file_path)
+hourtemp=function(latitude, date, temperature_file_path)
 {
     load(temperature_file_path)
     threshold <- 14.17  # base development threshold for BMSB
@@ -68,7 +70,7 @@
     }
     else
     {
-        dlprofile <- daylength(L)  # extract daylength data for entire year
+        dlprofile <- daylength(latitude)  # extract daylength data for entire year
         T <- NULL  # initialize hourly temperature
         dh <- NULL #initialize degree hour vector
         # date <- 200
@@ -365,7 +367,7 @@
                     if (vec.ind[4] == 0)
                     {
                         # just turned in vittelogenic stage
-                        n.birth=round(runif(1, 2 + min.ovi.adj, 8 + max.ovi.adj))
+                        n.birth=round(runif(1, 2 + opt$min_clutch_size, 8 + opt$max_clutch_size))
                     }
                     else
                     {
@@ -404,7 +406,7 @@
                     if (vec.ind[4] == 0)
                     {
                         # just turned in vittelogenic stage
-                        n.birth=round(runif(1, 2 + min.ovi.adj, 8 + max.ovi.adj))
+                        n.birth=round(runif(1, 2 + opt$min_clutch_size, 8 + opt$max_clutch_size))
                     }
                     else
                     {
@@ -443,7 +445,7 @@
                     # egg stage
                     # add to DD
                     vec.ind[3] <- vec.ind[3] + dd.temp
-                    if (vec.ind[3] >= (68 + dd.adj1))
+                    if (vec.ind[3] >= (68 + opt$young_nymph_accum))
                     {
                         # from egg to young nymph, DD requirement met
                         current.gen <- vec.ind[1]
@@ -464,7 +466,7 @@
                     # young nymph stage
                     # add to DD
                     vec.ind[3] <- vec.ind[3] + dd.temp
-                    if (vec.ind[3] >= (250 +dd.adj2))
+                    if (vec.ind[3] >= (250 + opt$old_nymph_accum))
                     {
                         # from young to old nymph, DD requirement met
                         current.gen <- vec.ind[1]
@@ -489,7 +491,7 @@
                     # old nymph stage
                     # add to DD
                     vec.ind[3] <- vec.ind[3] + dd.temp
-                    if (vec.ind[3] >= (200 + dd.adj3))
+                    if (vec.ind[3] >= (200 + opt$adult_accum))
                     {
                         # from old to adult, DD requirement met
                         current.gen <- vec.ind[1]