--- a/insect_phenology_model.R	Mon Dec 04 10:57:59 2017 -0500
+++ b/insect_phenology_model.R	Mon Dec 04 11:21:25 2017 -0500
@@ -201,35 +201,35 @@
 
 cat("Number of days: ", opt$num_days, "\n")
 
-# Initialize matrices for results from all replications.
-S0.rep <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
-S1.rep <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
-S2.rep <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
-S3.rep <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
-S4.rep <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
-S5.rep <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
-newborn.rep <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
-death.rep <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
-adult.rep <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
-pop.rep <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
-g0.rep <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
-g1.rep <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
-g2.rep <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
-g0a.rep <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
-g1a.rep <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
-g2a.rep <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
+# Initialize matrices.
+S0.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
+S1.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
+S2.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
+S3.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
+S4.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
+S5.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
+newborn.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
+death.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
+adult.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
+pop.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
+g0.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
+g1.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
+g2.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
+g0a.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
+g1a.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
+g2a.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
 
 # Loop through replications.
-for (N.rep in 1:opt$replications) {
+for (N.replications in 1:opt$replications) {
     # During each replication start with 1000 individuals.
     # TODO: user definable as well?
     num_insects <- 1000
     # Generation, Stage, degree-days, T, Diapause.
-    vec.ini <- c(0, 3, 0, 0, 0)
+    vector.ini <- c(0, 3, 0, 0, 0)
     # Overwintering, previttelogenic, degree-days=0, T=0, no-diapause.
-    vec.mat <- rep(vec.ini, num_insects)
+    vector.matrix <- rep(vector.ini, num_insects)
     # Complete matrix for the population.
-    vec.mat <- base::t(matrix(vec.mat, nrow=5))
+    vector.matrix <- base::t(matrix(vector.matrix, nrow=5))
     # Time series of population size.
     tot.pop <- NULL
     gen0.pop <- rep(0, opt$num_days)
@@ -250,13 +250,13 @@
         degree_days.temp <- temp.profile[2]
         degree_days.day[row] <- degree_days.temp
         # Trash bin for death.
-        death.vec <- NULL
+        death.vector <- NULL
         # Newborn.
-        birth.vec <- NULL
+        birth.vector <- NULL
         # All individuals.
         for (i in 1:num_insects) {
             # Find individual record.
-            vec.ind <- vec.mat[i,]
+            vector.ind <- vector.matrix[i,]
             # First of all, still alive?
             # Adjustment for late season mortality rate.
             if (latitude < 40.0) {
@@ -267,14 +267,14 @@
                 post.mort <- 2
                 day.kill <- 250
             }
-            if (vec.ind[2] == 0) {
+            if (vector.ind[2] == 0) {
                 # Egg.
                 death.prob = opt$egg_mort * mortality.egg(mean.temp)
             }
-            else if (vec.ind[2] == 1 | vec.ind[2] == 2) {
+            else if (vector.ind[2] == 1 | vector.ind[2] == 2) {
                 death.prob = opt$nymph_mort * mortality.nymph(mean.temp)
             }
-            else if (vec.ind[2] == 3 | vec.ind[2] == 4 | vec.ind[2] == 5) {
+            else if (vector.ind[2] == 3 | vector.ind[2] == 4 | vector.ind[2] == 5) {
                 # For adult.
                 if (doy < day.kill) {
                     death.prob = opt$adult_mort * mortality.adult(mean.temp)
@@ -287,48 +287,48 @@
             # (or dependent on temperature and life stage?)
             u.d <- runif(1)
             if (u.d < death.prob) {
-                death.vec <- c(death.vec, i)
+                death.vector <- c(death.vector, i)
             }
             else {
                 # Aggregrate index of dead bug.
                 # Event 1 end of diapause.
-                if (vec.ind[1] == 0 && vec.ind[2] == 3) {
+                if (vector.ind[1] == 0 && vector.ind[2] == 3) {
                     # Overwintering adult (previttelogenic).
-                    if (photoperiod > opt$photoperiod && vec.ind[3] > 68 && doy < 180) {
+                    if (photoperiod > opt$photoperiod && vector.ind[3] > 68 && doy < 180) {
                         # Add 68C to become fully reproductively matured.
                         # Transfer to vittelogenic.
-                        vec.ind <- c(0, 4, 0, 0, 0)
-                        vec.mat[i,] <- vec.ind
+                        vector.ind <- c(0, 4, 0, 0, 0)
+                        vector.matrix[i,] <- vector.ind
                     }
                     else {
                         # Add to degree_days.
-                        vec.ind[3] <- vec.ind[3] + degree_days.temp
+                        vector.ind[3] <- vector.ind[3] + degree_days.temp
                         # Add 1 day in current stage.
-                        vec.ind[4] <- vec.ind[4] + 1
-                        vec.mat[i,] <- vec.ind
+                        vector.ind[4] <- vector.ind[4] + 1
+                        vector.matrix[i,] <- vector.ind
                     }
                 }
-                if (vec.ind[1] != 0 && vec.ind[2] == 3) {
+                if (vector.ind[1] != 0 && vector.ind[2] == 3) {
                     # Not overwintering adult (previttelogenic).
-                    current.gen <- vec.ind[1]
-                    if (vec.ind[3] > 68) {
+                    current.gen <- vector.ind[1]
+                    if (vector.ind[3] > 68) {
                         # Add 68C to become fully reproductively matured.
                         # Transfer to vittelogenic.
-                        vec.ind <- c(current.gen, 4, 0, 0, 0)
-                        vec.mat[i,] <- vec.ind
+                        vector.ind <- c(current.gen, 4, 0, 0, 0)
+                        vector.matrix[i,] <- vector.ind
                     }
                     else {
                         # Add to degree_days.
-                        vec.ind[3] <- vec.ind[3] + degree_days.temp
+                        vector.ind[3] <- vector.ind[3] + degree_days.temp
                         # Add 1 day in current stage.
-                        vec.ind[4] <- vec.ind[4] + 1
-                        vec.mat[i,] <- vec.ind
+                        vector.ind[4] <- vector.ind[4] + 1
+                        vector.matrix[i,] <- vector.ind
                     }
                 }
                 # Event 2 oviposition -- where population dynamics comes from.
-                if (vec.ind[2] == 4 && vec.ind[1] == 0 && mean.temp > 10) {
+                if (vector.ind[2] == 4 && vector.ind[1] == 0 && mean.temp > 10) {
                     # Vittelogenic stage, overwintering generation.
-                    if (vec.ind[4] == 0) {
+                    if (vector.ind[4] == 0) {
                         # Just turned in vittelogenic stage.
                         num_insects.birth = round(runif(1, 2 + opt$min_clutch_size, 8 + opt$max_clutch_size))
                     }
@@ -341,26 +341,26 @@
                         }
                     }
                     # Add to degree_days.
-                    vec.ind[3] <- vec.ind[3] + degree_days.temp
+                    vector.ind[3] <- vector.ind[3] + degree_days.temp
                     # Add 1 day in current stage.
-                    vec.ind[4] <- vec.ind[4] + 1
-                    vec.mat[i,] <- vec.ind
+                    vector.ind[4] <- vector.ind[4] + 1
+                    vector.matrix[i,] <- vector.ind
                     if (num_insects.birth > 0) {
                         # Add new birth -- might be in different generations.
-                        new.gen <- vec.ind[1] + 1
+                        new.gen <- vector.ind[1] + 1
                         # Egg profile.
                         new.ind <- c(new.gen, 0, 0, 0, 0)
-                        new.vec <- rep(new.ind, num_insects.birth)
+                        new.vector <- rep(new.ind, num_insects.birth)
                         # Update batch of egg profile.
-                        new.vec <- t(matrix(new.vec, nrow=5))
+                        new.vector <- t(matrix(new.vector, nrow=5))
                         # Group with total eggs laid in that day.
-                        birth.vec <- rbind(birth.vec, new.vec)
+                        birth.vector <- rbind(birth.vector, new.vector)
                     }
                 }
                 # Event 2 oviposition -- for generation 1.
-                if (vec.ind[2] == 4 && vec.ind[1] == 1 && mean.temp > 12.5 && doy < 222) {
+                if (vector.ind[2] == 4 && vector.ind[1] == 1 && mean.temp > 12.5 && doy < 222) {
                     # Vittelogenic stage, 1st generation
-                    if (vec.ind[4] == 0) {
+                    if (vector.ind[4] == 0) {
                         # Just turned in vittelogenic stage.
                         num_insects.birth=round(runif(1, 2 + opt$min_clutch_size, 8 + opt$max_clutch_size))
                     }
@@ -373,125 +373,125 @@
                         }
                     }
                     # Add to degree_days.
-                    vec.ind[3] <- vec.ind[3] + degree_days.temp
+                    vector.ind[3] <- vector.ind[3] + degree_days.temp
                     # Add 1 day in current stage.
-                    vec.ind[4] <- vec.ind[4] + 1
-                    vec.mat[i,] <- vec.ind
+                    vector.ind[4] <- vector.ind[4] + 1
+                    vector.matrix[i,] <- vector.ind
                     if (num_insects.birth > 0) {
                         # Add new birth -- might be in different generations.
-                        new.gen <- vec.ind[1] + 1
+                        new.gen <- vector.ind[1] + 1
                         # Egg profile.
                         new.ind <- c(new.gen, 0, 0, 0, 0)
-                        new.vec <- rep(new.ind, num_insects.birth)
+                        new.vector <- rep(new.ind, num_insects.birth)
                         # Update batch of egg profile.
-                        new.vec <- t(matrix(new.vec, nrow=5))
+                        new.vector <- t(matrix(new.vector, nrow=5))
                         # Group with total eggs laid in that day.
-                        birth.vec <- rbind(birth.vec, new.vec)
+                        birth.vector <- rbind(birth.vector, new.vector)
                     }
                 }
                 # Event 3 development (with diapause determination).
-                # Event 3.1 egg development to young nymph (vec.ind[2]=0 -> egg).
-                if (vec.ind[2] == 0) {
+                # Event 3.1 egg development to young nymph (vector.ind[2]=0 -> egg).
+                if (vector.ind[2] == 0) {
                     # Egg stage.
                     # Add to degree_days.
-                    vec.ind[3] <- vec.ind[3] + degree_days.temp
-                    if (vec.ind[3] >= (68 + opt$young_nymph_accum)) {
+                    vector.ind[3] <- vector.ind[3] + degree_days.temp
+                    if (vector.ind[3] >= (68 + opt$young_nymph_accum)) {
                         # From egg to young nymph, degree-days requirement met.
-                        current.gen <- vec.ind[1]
+                        current.gen <- vector.ind[1]
                         # Transfer to young nymph stage.
-                        vec.ind <- c(current.gen, 1, 0, 0, 0)
+                        vector.ind <- c(current.gen, 1, 0, 0, 0)
                     }
                     else {
                         # Add 1 day in current stage.
-                        vec.ind[4] <- vec.ind[4] + 1
+                        vector.ind[4] <- vector.ind[4] + 1
                     }
-                    vec.mat[i,] <- vec.ind
+                    vector.matrix[i,] <- vector.ind
                 }
-                # Event 3.2 young nymph to old nymph (vec.ind[2]=1 -> young nymph: determines diapause).
-                if (vec.ind[2] == 1) {
+                # Event 3.2 young nymph to old nymph (vector.ind[2]=1 -> young nymph: determines diapause).
+                if (vector.ind[2] == 1) {
                     # Young nymph stage.
                     # Add to degree_days.
-                    vec.ind[3] <- vec.ind[3] + degree_days.temp
-                    if (vec.ind[3] >= (250 + opt$old_nymph_accum)) {
+                    vector.ind[3] <- vector.ind[3] + degree_days.temp
+                    if (vector.ind[3] >= (250 + opt$old_nymph_accum)) {
                         # From young to old nymph, degree_days requirement met.
-                        current.gen <- vec.ind[1]
+                        current.gen <- vector.ind[1]
                         # Transfer to old nym stage.
-                        vec.ind <- c(current.gen, 2, 0, 0, 0)
+                        vector.ind <- c(current.gen, 2, 0, 0, 0)
                         if (photoperiod < opt$photoperiod && doy > 180) {
-                            vec.ind[5] <- 1
+                            vector.ind[5] <- 1
                         } # Prepare for diapausing.
                     }
                     else {
                         # Add 1 day in current stage.
-                        vec.ind[4] <- vec.ind[4] + 1
+                        vector.ind[4] <- vector.ind[4] + 1
                     }
-                    vec.mat[i,] <- vec.ind
+                    vector.matrix[i,] <- vector.ind
                 }
                 # Event 3.3 old nymph to adult: previttelogenic or diapausing?
-                if (vec.ind[2] == 2) {
+                if (vector.ind[2] == 2) {
                     # Old nymph stage.
                     # Add to degree_days.
-                    vec.ind[3] <- vec.ind[3] + degree_days.temp
-                    if (vec.ind[3] >= (200 + opt$adult_accum)) {
+                    vector.ind[3] <- vector.ind[3] + degree_days.temp
+                    if (vector.ind[3] >= (200 + opt$adult_accum)) {
                         # From old to adult, degree_days requirement met.
-                        current.gen <- vec.ind[1]
-                        if (vec.ind[5] == 0) {
+                        current.gen <- vector.ind[1]
+                        if (vector.ind[5] == 0) {
                             # Non-diapausing adult -- previttelogenic.
-                            vec.ind <- c(current.gen, 3, 0, 0, 0)
+                            vector.ind <- c(current.gen, 3, 0, 0, 0)
                         }
                         else {
                             # Diapausing.
-                            vec.ind <- c(current.gen, 5, 0, 0, 1)
+                            vector.ind <- c(current.gen, 5, 0, 0, 1)
                         }
                     }
                     else {
                         # Add 1 day in current stage.
-                        vec.ind[4] <- vec.ind[4] + 1
+                        vector.ind[4] <- vector.ind[4] + 1
                     }
-                    vec.mat[i,] <- vec.ind
+                    vector.matrix[i,] <- vector.ind
                 }
                 # Event 4 growing of diapausing adult (unimportant, but still necessary).
-                if (vec.ind[2] == 5) {
-                    vec.ind[3] <- vec.ind[3] + degree_days.temp
-                    vec.ind[4] <- vec.ind[4] + 1
-                    vec.mat[i,] <- vec.ind
+                if (vector.ind[2] == 5) {
+                    vector.ind[3] <- vector.ind[3] + degree_days.temp
+                    vector.ind[4] <- vector.ind[4] + 1
+                    vector.matrix[i,] <- vector.ind
                 }
             } # Else if it is still alive.
         } # End of the individual bug loop.
         # Find how many died.
-        num_insects.death <- length(death.vec)
+        num_insects.death <- length(death.vector)
         if (num_insects.death > 0) {
-            vec.mat <- vec.mat[-death.vec, ]
+            vector.matrix <- vector.matrix[-death.vector, ]
         }
         # Remove record of dead.
         # Find how many new born.
-        num_insects.newborn <- length(birth.vec[,1])
-        vec.mat <- rbind(vec.mat, birth.vec)
+        num_insects.newborn <- length(birth.vector[,1])
+        vector.matrix <- rbind(vector.matrix, birth.vector)
         # Update population size for the next day.
         num_insects <- num_insects - num_insects.death + num_insects.newborn
 
         # Aggregate results by day.
         tot.pop <- c(tot.pop, num_insects)
         # Egg.
-        sum_eggs <- sum(vec.mat[,2] == 0)
+        sum_eggs <- sum(vector.matrix[,2] == 0)
         # Young nymph.
-        sum_young_nymphs <- sum(vec.mat[,2] == 1)
+        sum_young_nymphs <- sum(vector.matrix[,2] == 1)
         # Old nymph.
-        sum_old_nymphs <- sum(vec.mat[,2] == 2)
+        sum_old_nymphs <- sum(vector.matrix[,2] == 2)
         # Previtellogenic.
-        sum_previtellogenic <- sum(vec.mat[,2] == 3)
+        sum_previtellogenic <- sum(vector.matrix[,2] == 3)
         # Vitellogenic.
-        sum_vitellogenic <- sum(vec.mat[,2] == 4)
+        sum_vitellogenic <- sum(vector.matrix[,2] == 4)
         # Diapausing.
-        sum_diapausing <- sum(vec.mat[,2] == 5)
+        sum_diapausing <- sum(vector.matrix[,2] == 5)
         # Overwintering adult.
-        gen0 <- sum(vec.mat[,1] == 0)
+        gen0 <- sum(vector.matrix[,1] == 0)
         # First generation.
-        gen1 <- sum(vec.mat[,1] == 1)
+        gen1 <- sum(vector.matrix[,1] == 1)
         # Second generation.
-        gen2 <- sum(vec.mat[,1] == 2)
+        gen2 <- sum(vector.matrix[,1] == 2)
         # Sum of all adults.
-        num_insects.adult <- sum(vec.mat[,2] == 3) + sum(vec.mat[,2] == 4) + sum(vec.mat[,2] == 5)
+        num_insects.adult <- sum(vector.matrix[,2] == 3) + sum(vector.matrix[,2] == 4) + sum(vector.matrix[,2] == 5)
 
         # Generation 0 pop size.
         gen0.pop[row] <- gen0
@@ -505,9 +505,9 @@
         S4[row] <- sum_vitellogenic
         S5[row] <- sum_diapausing
 
-        g0.adult[row] <- sum(vec.mat[,1] == 0)
-        g1.adult[row] <- sum((vec.mat[,1] == 1 & vec.mat[,2] == 3) | (vec.mat[,1] == 1 & vec.mat[,2] == 4) | (vec.mat[,1] == 1 & vec.mat[,2] == 5))
-        g2.adult[row] <- sum((vec.mat[,1]== 2 & vec.mat[,2] == 3) | (vec.mat[,1] == 2 & vec.mat[,2] == 4) | (vec.mat[,1] == 2 & vec.mat[,2] == 5))
+        g0.adult[row] <- sum(vector.matrix[,1] == 0)
+        g1.adult[row] <- sum((vector.matrix[,1] == 1 & vector.matrix[,2] == 3) | (vector.matrix[,1] == 1 & vector.matrix[,2] == 4) | (vector.matrix[,1] == 1 & vector.matrix[,2] == 5))
+        g2.adult[row] <- sum((vector.matrix[,1]== 2 & vector.matrix[,2] == 3) | (vector.matrix[,1] == 2 & vector.matrix[,2] == 4) | (vector.matrix[,1] == 2 & vector.matrix[,2] == 5))
 
         N.newborn[row] <- num_insects.newborn
         N.death[row] <- num_insects.death
@@ -517,61 +517,61 @@
     degree_days.cum <- cumsum(degree_days.day)
 
     # Collect all the outputs.
-    S0.rep[,N.rep] <- S0
-    S1.rep[,N.rep] <- S1
-    S2.rep[,N.rep] <- S2
-    S3.rep[,N.rep] <- S3
-    S4.rep[,N.rep] <- S4
-    S5.rep[,N.rep] <- S5
-    newborn.rep[,N.rep] <- N.newborn
-    death.rep[,N.rep] <- N.death
-    adult.rep[,N.rep] <- N.adult
-    pop.rep[,N.rep] <- tot.pop
-    g0.rep[,N.rep] <- gen0.pop
-    g1.rep[,N.rep] <- gen1.pop
-    g2.rep[,N.rep] <- gen2.pop
-    g0a.rep[,N.rep] <- g0.adult
-    g1a.rep[,N.rep] <- g1.adult
-    g2a.rep[,N.rep] <- g2.adult
+    S0.replications[,N.replications] <- S0
+    S1.replications[,N.replications] <- S1
+    S2.replications[,N.replications] <- S2
+    S3.replications[,N.replications] <- S3
+    S4.replications[,N.replications] <- S4
+    S5.replications[,N.replications] <- S5
+    newborn.replications[,N.replications] <- N.newborn
+    death.replications[,N.replications] <- N.death
+    adult.replications[,N.replications] <- N.adult
+    pop.replications[,N.replications] <- tot.pop
+    g0.replications[,N.replications] <- gen0.pop
+    g1.replications[,N.replications] <- gen1.pop
+    g2.replications[,N.replications] <- gen2.pop
+    g0a.replications[,N.replications] <- g0.adult
+    g1a.replications[,N.replications] <- g1.adult
+    g2a.replications[,N.replications] <- g2.adult
 }
 
 # Mean value for adults.
-mean_value_adult <- apply((S3.rep + S4.rep + S5.rep), 1, mean)
+adult_mean <- apply((S3.replications + S4.replications + S5.replications), 1, mean)
 # Mean value for nymphs.
-mean_value_nymph <- apply((S1.rep + S2.rep), 1, mean)
+nymph_mean <- apply((S1.replications + S2.replications), 1, mean)
 # Mean value for eggs.
-mean_value_egg <- apply(S0.rep, 1, mean)
+egg_mean <- apply(S0.replications, 1, mean)
 # Mean value for P.
-g0 <- apply(g0.rep, 1, mean)
+g0 <- apply(g0.replications, 1, mean)
 # Mean value for F1.
-g1 <- apply(g1.rep, 1, mean)
+g1 <- apply(g1.replications, 1, mean)
 # Mean value for F2.
-g2 <- apply(g2.rep, 1, mean)
+g2 <- apply(g2.replications, 1, mean)
 # Mean value for P adult.
-g0a <- apply(g0a.rep, 1, mean)
+g0a <- apply(g0a.replications, 1, mean)
 # Mean value for F1 adult.
-g1a <- apply(g1a.rep, 1, mean)
+g1a <- apply(g1a.replications, 1, mean)
 # Mean value for F2 adult.
-g2a <- apply(g2a.rep, 1, mean)
+g2a <- apply(g2a.replications, 1, mean)
 
 # Standard error for adults.
-mean_value_adult.std_error <- apply((S3.rep + S4.rep + S5.rep), 1, sd) / sqrt(opt$replications)
+adult_mean.std_error <- apply((S3.replications + S4.replications + S5.replications), 1, sd) / sqrt(opt$replications)
 # Standard error for nymphs.
-mean_value_nymph.std_error <- apply((S1.rep + S2.rep) / sqrt(opt$replications), 1, sd)
+nymph_mean.std_error <- apply((S1.replications + S2.replications) / sqrt(opt$replications), 1, sd)
 # Standard error for eggs.
-mean_value_egg.std_error <- apply(S0.rep, 1, sd) / sqrt(opt$replications)
+egg_mean.std_error <- apply(S0.replications, 1, sd) / sqrt(opt$replications)
 # Standard error value for P.
-g0.std_error <- apply(g0.rep, 1, sd) / sqrt(opt$replications)
+g0.std_error <- apply(g0.replications, 1, sd) / sqrt(opt$replications)
 # Standard error for F1.
-g1.std_error <- apply(g1.rep, 1, sd) / sqrt(opt$replications)
+g1.std_error <- apply(g1.replications, 1, sd) / sqrt(opt$replications)
 # Standard error for F2.
-g2.std_error <- apply(g2.rep, 1, sd) / sqrt(opt$replications)
+g2.std_error <- apply(g2.replications, 1, sd) / sqrt(opt$replications)
 # Standard error for P adult.
-g0a.std_error <- apply(g0a.rep, 1, sd) / sqrt(opt$replications)
+g0a.std_error <- apply(g0a.replications, 1, sd) / sqrt(opt$replications)
 # Standard error for F1 adult.
-g1a.std_error <- apply(g1a.rep, 1, sd) / sqrt(opt$replications)
+g1a.std_error <- apply(g1a.replications, 1, sd) / sqrt(opt$replications)
 # Standard error for F2 adult.
-g2a.std_error <- apply(g2a.rep, 1, sd) / sqrt(opt$replications)
+g2a.std_error <- apply(g2a.replications, 1, sd) / sqrt(opt$replications)
 
 dev.new(width=20, height=30)
 
@@ -588,26 +588,25 @@
 
 # Subfigure 1: population size by life stage
 title <- paste(opt$insect, ": Total pop. by life stage :", opt$location, ": Lat:", latitude, ":", start_date, "-", end_date, sep=" ")
-plot(day.all, mean_value_adult, main=title, type="l", ylim=c(0, max(mean_value_egg + mean_value_egg.std_error, mean_value_nymph + mean_value_nymph.std_error, mean_value_adult + mean_value_adult.std_error)), axes=F, lwd=2, xlab="", ylab="", cex=3, cex.lab=3, cex.axis=3, cex.main=3)
+plot(day.all, adult_mean, main=title, type="l", ylim=c(0, max(egg_mean + egg_mean.std_error, nymph_mean + nymph_mean.std_error, adult_mean + adult_mean.std_error)), axes=F, lwd=2, xlab="", ylab="", cex=3, cex.lab=3, cex.axis=3, cex.main=3)
 # Young and old nymphs.
-lines(day.all, mean_value_nymph, lwd=2, lty=1, col=2)
+lines(day.all, nymph_mean, lwd=2, lty=1, col=2)
 # Eggs
-lines(day.all, mean_value_egg, lwd=2, lty=1, col=4)
+lines(day.all, egg_mean, lwd=2, lty=1, col=4)
 axis(1, at=c(1:12) * 30 - 15, cex.axis=3, labels=c("Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"))
 axis(2, cex.axis=3)
-leg.text <- c("Egg", "Nymph", "Adult")
-legend("topleft", leg.text, lty=c(1, 1, 1), col=c(4, 2, 1), cex=3)
+legend("topleft", c("Egg", "Nymph", "Adult"), lty=c(1, 1, 1), col=c(4, 2, 1), cex=3)
 if (opt$std_error_plot == 1) {
     # Add Standard error lines to plot
     # Standard error for adults
-    lines (day.all, mean_value_adult+mean_value_adult.std_error, lty=2)
-    lines (day.all, mean_value_adult-mean_value_adult.std_error, lty=2)
+    lines (day.all, adult_mean+adult_mean.std_error, lty=2)
+    lines (day.all, adult_mean-adult_mean.std_error, lty=2)
     # Standard error for nymphs
-    lines (day.all, mean_value_nymph+mean_value_nymph.std_error, col=2, lty=2)
-    lines (day.all, mean_value_nymph-mean_value_nymph.std_error, col=2, lty=2)
+    lines (day.all, nymph_mean+nymph_mean.std_error, col=2, lty=2)
+    lines (day.all, nymph_mean-nymph_mean.std_error, col=2, lty=2)
     # Standard error for eggs
-    lines (day.all, mean_value_egg+mean_value_egg.std_error, col=4, lty=2)
-    lines (day.all, mean_value_egg-mean_value_egg.std_error, col=4, lty=2)
+    lines (day.all, egg_mean+egg_mean.std_error, col=4, lty=2)
+    lines (day.all, egg_mean-egg_mean.std_error, col=4, lty=2)
 }
 
 # Subfigure 2: population size by generation
@@ -617,8 +616,7 @@
 lines(day.all, g2, lwd = 2, lty = 1, col=4)
 axis(1, at=c(1:12) * 30 - 15, cex.axis=3, labels = c("Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"))
 axis(2, cex.axis=3)
-leg.text <- c("P", "F1", "F2")
-legend("topleft", leg.text, lty=c(1, 1, 1), col=c(1, 2, 4), cex=3)
+legend("topleft", c("P", "F1", "F2"), lty=c(1, 1, 1), col=c(1, 2, 4), cex=3)
 if (opt$std_error_plot == 1) {
     # Add Standard error lines to plot
     # Standard error for adults
@@ -639,8 +637,7 @@
 lines(day.all, g2a, lwd = 2, lty = 1, col=4)
 axis(1, at=c(1:12) * 30 - 15, cex.axis=3, labels = c("Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"))
 axis(2, cex.axis=3)
-leg.text <- c("P", "F1", "F2")
-legend("topleft", leg.text, lty=c(1, 1, 1), col=c(1, 2, 4), cex=3)
+legend("topleft", c("P", "F1", "F2"), lty=c(1, 1, 1), col=c(1, 2, 4), cex=3)
 if (opt$std_error_plot == 1) {
     # Add Standard error lines to plot
     # Standard error for adults