--- a/insect_phenology_model.R	Mon Aug 14 13:44:45 2017 -0400
+++ b/insect_phenology_model.R	Thu Sep 28 10:45:53 2017 -0400
@@ -12,6 +12,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("-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"),
     make_option(c("-q", "--photoperiod"), action="store", dest="photoperiod", type="double", help="Critical photoperiod for diapause induction/termination"),
@@ -28,36 +29,39 @@
 
 data.input=function(loc, year, temperature.dataset)
 {
-    expdata <- matrix(rep(0, 365 * 3), nrow=365)
+    expdata <- matrix(rep(0, opt$num_days * 3), nrow=opt$num_days)
     namedat <- paste(loc,  year, ".Rdat", sep="")
     temp.data <- read.csv(file=temperature.dataset, header=T)
 
-    expdata[,1] <- c(1:365)
+    expdata[,1] <- c(1:opt$num_days)
     # Minimum
-    expdata[,2] <- temp.data[c(1:365), 3]
+    expdata[,2] <- temp.data[c(1:opt$num_days), 3]
     # Maximum
-    expdata[,3] <- temp.data[c(1:365), 2]
+    expdata[,3] <- temp.data[c(1:opt$num_days), 2]
     save(expdata, file=namedat)
     namedat
 }
 
-daylength=function(latitude)
+daylength=function(latitude, num_days)
 {
     # from Forsythe 1995
     p=0.8333
     dl <- NULL
-    for (i in 1:365) {
+    for (i in 1:num_days) {
         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)))
     }
-    dl   # return a vector of daylength in 365 days
+    # return a vector of daylength for the number of
+    # days specifie din the input temperature data
+    dl
 }
 
-hourtemp=function(latitude, date, temperature_file_path)
+hourtemp=function(latitude, date, temperature_file_path, num_days)
 {
     load(temperature_file_path)
-    threshold <- 14.17  # base development threshold for BMSB
+    # base development threshold for BMSB
+    threshold <- 14.17
     dnp <- expdata[date, 2]  # daily minimum
     dxp <- expdata[date, 3]  # daily maximum
     dmean <- 0.5 * (dnp + dxp)
@@ -67,22 +71,28 @@
         dd <- 0
     }
     else {
-        dlprofile <- daylength(latitude)  # extract daylength data for entire year
+        # extract daylength data for the number of
+        # days specified in the input temperature data
+        dlprofile <- daylength(latitude, num_days)
         T <- NULL  # initialize hourly temperature
         dh <- NULL #initialize degree hour vector
-        # date <- 200
-        y <- dlprofile[date]  # calculate daylength in given date
-        z <- 24 - y     # night length
-        a <- 1.86     # lag coefficient
-        b <- 2.20     # night coefficient
-        #tempdata <- read.csv("tempdata.csv") #import raw data set
-        # Should be outside function otherwise its redundant
-        risetime <- 12 - y / 2      # sunrise time
-        settime <- 12 + y / 2       # sunset time
+        # calculate daylength in given date
+        y <- dlprofile[date]
+        # night length
+        z <- 24 - y
+        # lag coefficient
+        a <- 1.86
+        # night coefficient
+        b <- 2.20
+        # sunrise time
+        risetime <- 12 - y / 2
+        # sunset time
+        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) {
-                m <- i - 5  # number of hours after Tmin until sunset
+                # 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) {
                     dh[i] <- 0
@@ -205,8 +215,8 @@
 temperature_file_path <- data.input(opt$location, opt$year, opt$temperature_dataset)
 
 # Initialize matrix for results from all replications
-S0.rep <- S1.rep <- S2.rep <- S3.rep <- S4.rep <- S5.rep <- matrix(rep(0, 365 * opt$replications), ncol = opt$replications)
-newborn.rep <- death.rep <- adult.rep <- pop.rep <- g0.rep <- g1.rep <- g2.rep <- g0a.rep <- g1a.rep <- g2a.rep <- matrix(rep(0, 365 * opt$replications), ncol=opt$replications)
+S0.rep <- S1.rep <- S2.rep <- S3.rep <- S4.rep <- S5.rep <- matrix(rep(0, opt$num_days * opt$replications), ncol = opt$replications)
+newborn.rep <- death.rep <- adult.rep <- pop.rep <- g0.rep <- g1.rep <- g2.rep <- g0a.rep <- g1a.rep <- g2a.rep <- matrix(rep(0, opt$num_days * opt$replications), ncol=opt$replications)
 
 # loop through replications
 for (N.rep in 1:opt$replications) {
@@ -220,27 +230,27 @@
     # complete matrix for the population
     vec.mat <- t(matrix(vec.mat, nrow=5))
     # complete photoperiod profile in a year, requires daylength function
-    ph.p <- daylength(opt$latitude)
+    ph.p <- daylength(opt$latitude, opt$num_days)
 
     # time series of population size
     tot.pop <- NULL
     # gen.0 pop size
-    gen0.pop <- rep(0, 365)
-    gen1.pop <- rep(0, 365)
-    gen2.pop <- rep(0, 365)
-    S0 <- S1 <- S2 <- S3 <- S4 <- S5 <- rep(0, 365)
-    g0.adult <- g1.adult <- g2.adult <- rep(0, 365)
-    N.newborn <- N.death <- N.adult <- rep(0, 365)
-    dd.day <- rep(0, 365)
+    gen0.pop <- rep(0, opt$num_days)
+    gen1.pop <- rep(0, opt$num_days)
+    gen2.pop <- rep(0, opt$num_days)
+    S0 <- S1 <- S2 <- S3 <- S4 <- S5 <- rep(0, opt$num_days)
+    g0.adult <- g1.adult <- g2.adult <- rep(0, opt$num_days)
+    N.newborn <- N.death <- N.adult <- rep(0, opt$num_days)
+    dd.day <- rep(0, opt$num_days)
 
     # start tick
     ptm <- proc.time()
 
     # all the days
-    for (day in 1:365) {
+    for (day in 1:opt$num_days) {
         # photoperiod in the day
         photoperiod <- ph.p[day]
-        temp.profile <- hourtemp(opt$latitude, day, temperature_file_path)
+        temp.profile <- hourtemp(opt$latitude, day, temperature_file_path, opt$num_days)
         mean.temp <- temp.profile[1]
         dd.temp <- temp.profile[2]
         dd.day[day] <- dd.temp
@@ -514,8 +524,7 @@
         N.newborn[day] <- n.newborn
         N.death[day] <- n.death
         N.adult[day] <- n.adult
-        #print(c(N.rep, day, n, n.adult))
-    }   # end of 365 days
+    }   # end of days specified in the input temperature data
 
     dd.cum <- cumsum(dd.day)
     # collect all the outputs
@@ -537,15 +546,11 @@
     g2a.rep[,N.rep] <- g2.adult
 }
 
-# save(dd.day, dd.cum, S0.rep, S1.rep, S2.rep, S3.rep, S4.rep, S5.rep, newborn.rep, death.rep, adult.rep, pop.rep, g0.rep, g1.rep, g2.rep, g0a.rep, g1a.rep, g2a.rep, file=opt$output)
-# maybe do not need to export this bit, but for now just leave it as-is
-# do we need to export this Rdat file? 
-
 # 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:365 * n.yr)
+day.all <- c(1:opt$num_days * n.yr)
 
 # mean value for adults
 sa <- apply((S3.rep + S4.rep + S5.rep), 1, mean)