5
+ − 1 #!/usr/bin/env Rscript
+ − 2
+ − 3 suppressPackageStartupMessages(library("optparse"))
+ − 4
+ − 5 option_list <- list(
6
+ − 6 make_option(c("--adult_mortality"), action="store", dest="adult_mortality", type="integer", help="Adjustment rate for adult mortality"),
+ − 7 make_option(c("--adult_accumulation"), action="store", dest="adult_accumulation", type="integer", help="Adjustment of degree-days accumulation (old nymph->adult)"),
+ − 8 make_option(c("--egg_mortality"), action="store", dest="egg_mortality", type="integer", help="Adjustment rate for egg mortality"),
38
+ − 9 make_option(c("--input_norm"), action="store", dest="input_norm", help="30 year normals temperature data for selected station"),
43
+ − 10 make_option(c("--input_ytd"), action="store", dest="input_ytd", default=NULL, help="Year-to-date temperature data for selected location"),
6
+ − 11 make_option(c("--insect"), action="store", dest="insect", help="Insect name"),
+ − 12 make_option(c("--insects_per_replication"), action="store", dest="insects_per_replication", type="integer", help="Number of insects with which to start each replication"),
10
+ − 13 make_option(c("--life_stages"), action="store", dest="life_stages", help="Selected life stages for plotting"),
+ − 14 make_option(c("--life_stages_adult"), action="store", dest="life_stages_adult", default=NULL, help="Adult life stages for plotting"),
16
+ − 15 make_option(c("--life_stages_nymph"), action="store", dest="life_stages_nymph", default=NULL, help="Nymph life stages for plotting"),
6
+ − 16 make_option(c("--location"), action="store", dest="location", help="Selected location"),
+ − 17 make_option(c("--min_clutch_size"), action="store", dest="min_clutch_size", type="integer", help="Adjustment of minimum clutch size"),
+ − 18 make_option(c("--max_clutch_size"), action="store", dest="max_clutch_size", type="integer", help="Adjustment of maximum clutch size"),
43
+ − 19 make_option(c("--num_days_ytd"), action="store", dest="num_days_ytd", default=NULL, type="integer", help="Total number of days in the year-to-date temperature dataset"),
6
+ − 20 make_option(c("--nymph_mortality"), action="store", dest="nymph_mortality", type="integer", help="Adjustment rate for nymph mortality"),
+ − 21 make_option(c("--old_nymph_accumulation"), action="store", dest="old_nymph_accumulation", type="integer", help="Adjustment of degree-days accumulation (young nymph->old nymph)"),
+ − 22 make_option(c("--oviposition"), action="store", dest="oviposition", type="integer", help="Adjustment for oviposition rate"),
+ − 23 make_option(c("--photoperiod"), action="store", dest="photoperiod", type="double", help="Critical photoperiod for diapause induction/termination"),
10
+ − 24 make_option(c("--plot_generations_separately"), action="store", dest="plot_generations_separately", help="Plot Plot P, F1 and F2 as separate lines or pool across them"),
+ − 25 make_option(c("--plot_std_error"), action="store", dest="plot_std_error", help="Plot Standard error"),
27
+ − 26 make_option(c("--replications"), action="store", dest="replications", type="integer", help="Number of replications"),
6
+ − 27 make_option(c("--young_nymph_accumulation"), action="store", dest="young_nymph_accumulation", type="integer", help="Adjustment of degree-days accumulation (egg->young nymph)")
5
+ − 28 )
+ − 29
8
+ − 30 parser <- OptionParser(usage="%prog [options] file", option_list=option_list);
+ − 31 args <- parse_args(parser, positional_arguments=TRUE);
+ − 32 opt <- args$options;
5
+ − 33
27
+ − 34 add_daylight_length = function(temperature_data_frame, num_rows) {
5
+ − 35 # Return a vector of daylight length (photoperido profile) for
38
+ − 36 # the number of days specified in the input_ytd temperature data
5
+ − 37 # (from Forsythe 1995).
8
+ − 38 p = 0.8333;
+ − 39 latitude = temperature_data_frame$LATITUDE[1];
+ − 40 daylight_length_vector = NULL;
5
+ − 41 for (i in 1:num_rows) {
+ − 42 # Get the day of the year from the current row
+ − 43 # of the temperature data for computation.
8
+ − 44 doy = temperature_data_frame$DOY[i];
+ − 45 theta = 0.2163108 + 2 * atan(0.9671396 * tan(0.00860 * (doy - 186)));
+ − 46 phi = asin(0.39795 * cos(theta));
5
+ − 47 # Compute the length of daylight for the day of the year.
8
+ − 48 darkness_length = 24 / pi * acos((sin(p * pi / 180) + sin(latitude * pi / 180) * sin(phi)) / (cos(latitude * pi / 180) * cos(phi)));
+ − 49 daylight_length_vector[i] = 24 - darkness_length;
5
+ − 50 }
+ − 51 # Append daylight_length_vector as a new column to temperature_data_frame.
27
+ − 52 temperature_data_frame = append_vector(temperature_data_frame, daylight_length_vector, "DAYLEN");
8
+ − 53 return(temperature_data_frame);
5
+ − 54 }
+ − 55
27
+ − 56 append_vector = function(data_frame, vec, new_column_name) {
+ − 57 num_columns = dim(data_frame)[2];
+ − 58 current_column_names = colnames(data_frame);
+ − 59 # Append vector vec as a new column to data_frame.
+ − 60 data_frame[,num_columns+1] = vec;
+ − 61 # Reset the column names with the additional column for later access.
+ − 62 colnames(data_frame) = append(current_column_names, new_column_name);
+ − 63 return(data_frame);
+ − 64 }
+ − 65
19
+ − 66 get_file_path = function(life_stage, base_name, life_stage_nymph=NULL, life_stage_adult=NULL) {
+ − 67 if (!is.null(life_stage_nymph)) {
+ − 68 lsi = get_life_stage_index(life_stage, life_stage_nymph=life_stage_nymph);
+ − 69 file_name = paste(lsi, tolower(life_stage_nymph), base_name, sep="_");
+ − 70 } else if (!is.null(life_stage_adult)) {
+ − 71 lsi = get_life_stage_index(life_stage, life_stage_adult=life_stage_adult);
+ − 72 file_name = paste(lsi, tolower(life_stage_adult), base_name, sep="_");
+ − 73 } else {
+ − 74 lsi = get_life_stage_index(life_stage);
+ − 75 file_name = paste(lsi, base_name, sep="_");
+ − 76 }
34
+ − 77 file_path = paste("output_plots_dir", file_name, sep="/");
19
+ − 78 return(file_path);
+ − 79 }
+ − 80
18
+ − 81 get_life_stage_index = function(life_stage, life_stage_nymph=NULL, life_stage_adult=NULL) {
+ − 82 # Name collection elements so that they
+ − 83 # are displayed in logical order.
+ − 84 if (life_stage=="Egg") {
+ − 85 lsi = "01";
+ − 86 } else if (life_stage=="Nymph") {
+ − 87 if (life_stage_nymph=="Young") {
+ − 88 lsi = "02";
+ − 89 } else if (life_stage_nymph=="Old") {
+ − 90 lsi = "03";
+ − 91 } else if (life_stage_nymph=="Total") {
+ − 92 lsi="04";
+ − 93 }
+ − 94 } else if (life_stage=="Adult") {
+ − 95 if (life_stage_adult=="Pre-vittelogenic") {
+ − 96 lsi = "05";
+ − 97 } else if (life_stage_adult=="Vittelogenic") {
+ − 98 lsi = "06";
+ − 99 } else if (life_stage_adult=="Diapausing") {
+ − 100 lsi = "07";
+ − 101 } else if (life_stage_adult=="Total") {
+ − 102 lsi = "08";
+ − 103 }
+ − 104 } else if (life_stage=="Total") {
+ − 105 lsi = "09";
+ − 106 }
+ − 107 return(lsi);
+ − 108 }
+ − 109
20
+ − 110 get_mean_and_std_error = function(p_replications, f1_replications, f2_replications) {
+ − 111 # P mean.
+ − 112 p_m = apply(p_replications, 1, mean);
+ − 113 # P standard error.
+ − 114 p_se = apply(p_replications, 1, sd) / sqrt(opt$replications);
+ − 115 # F1 mean.
+ − 116 f1_m = apply(f1_replications, 1, mean);
+ − 117 # F1 standard error.
+ − 118 f1_se = apply(f1_replications, 1, sd) / sqrt(opt$replications);
+ − 119 # F2 mean.
+ − 120 f2_m = apply(f2_replications, 1, mean);
+ − 121 # F2 standard error.
+ − 122 f2_se = apply(f2_replications, 1, sd) / sqrt(opt$replications);
+ − 123 return(list(p_m, p_se, f1_m, f1_se, f2_m, f2_se))
+ − 124 }
+ − 125
39
+ − 126 get_next_normals_row = function(norm_data_frame, year, is_leap_year, index) {
+ − 127 # Return the next 30 year normals row formatted
+ − 128 # appropriately for the year-to-date data.
+ − 129 latitude = norm_data_frame[index,"LATITUDE"][1];
+ − 130 longitude = norm_data_frame[index,"LONGITUDE"][1];
+ − 131 # Format the date.
+ − 132 mmdd = norm_data_frame[index,"MMDD"][1];
+ − 133 date_str = paste(year, mmdd, sep="-");
+ − 134 doy = norm_data_frame[index,"DOY"][1];
+ − 135 if (!is_leap_year) {
+ − 136 # Since all normals data includes Feb 29, we have to
+ − 137 # subtract 1 from DOY if we're not in a leap year since
+ − 138 # we removed the Feb 29 row from the data frame above.
+ − 139 doy = as.integer(doy) - 1;
+ − 140 }
+ − 141 tmin = norm_data_frame[index,"TMIN"][1];
+ − 142 tmax = norm_data_frame[index,"TMAX"][1];
+ − 143 return(list(latitude, longitude, date_str, doy, tmin, tmax));
+ − 144 }
+ − 145
5
+ − 146 get_temperature_at_hour = function(latitude, temperature_data_frame, row, num_days) {
8
+ − 147 # Base development threshold for Brown Marmorated Stink Bug
5
+ − 148 # insect phenology model.
8
+ − 149 threshold = 14.17;
5
+ − 150 # Minimum temperature for current row.
8
+ − 151 curr_min_temp = temperature_data_frame$TMIN[row];
5
+ − 152 # Maximum temperature for current row.
8
+ − 153 curr_max_temp = temperature_data_frame$TMAX[row];
5
+ − 154 # Mean temperature for current row.
8
+ − 155 curr_mean_temp = 0.5 * (curr_min_temp + curr_max_temp);
5
+ − 156 # Initialize degree day accumulation
8
+ − 157 averages = 0;
6
+ − 158 if (curr_max_temp < threshold) {
8
+ − 159 averages = 0;
5
+ − 160 }
+ − 161 else {
+ − 162 # Initialize hourly temperature.
8
+ − 163 T = NULL;
5
+ − 164 # Initialize degree hour vector.
8
+ − 165 dh = NULL;
5
+ − 166 # Daylight length for current row.
8
+ − 167 y = temperature_data_frame$DAYLEN[row];
5
+ − 168 # Darkness length.
8
+ − 169 z = 24 - y;
5
+ − 170 # Lag coefficient.
8
+ − 171 a = 1.86;
5
+ − 172 # Darkness coefficient.
8
+ − 173 b = 2.20;
5
+ − 174 # Sunrise time.
8
+ − 175 risetime = 12 - y / 2;
5
+ − 176 # Sunset time.
8
+ − 177 settime = 12 + y / 2;
+ − 178 ts = (curr_max_temp - curr_min_temp) * sin(pi * (settime - 5) / (y + 2 * a)) + curr_min_temp;
5
+ − 179 for (i in 1:24) {
+ − 180 if (i > risetime && i < settime) {
+ − 181 # Number of hours after Tmin until sunset.
8
+ − 182 m = i - 5;
+ − 183 T[i] = (curr_max_temp - curr_min_temp) * sin(pi * m / (y + 2 * a)) + curr_min_temp;
5
+ − 184 if (T[i] < 8.4) {
8
+ − 185 dh[i] = 0;
5
+ − 186 }
+ − 187 else {
8
+ − 188 dh[i] = T[i] - 8.4;
5
+ − 189 }
+ − 190 }
6
+ − 191 else if (i > settime) {
8
+ − 192 n = i - settime;
+ − 193 T[i] = curr_min_temp + (ts - curr_min_temp) * exp( - b * n / z);
5
+ − 194 if (T[i] < 8.4) {
8
+ − 195 dh[i] = 0;
5
+ − 196 }
+ − 197 else {
8
+ − 198 dh[i] = T[i] - 8.4;
5
+ − 199 }
+ − 200 }
+ − 201 else {
8
+ − 202 n = i + 24 - settime;
+ − 203 T[i] = curr_min_temp + (ts - curr_min_temp) * exp( - b * n / z);
5
+ − 204 if (T[i] < 8.4) {
8
+ − 205 dh[i] = 0;
5
+ − 206 }
+ − 207 else {
8
+ − 208 dh[i] = T[i] - 8.4;
5
+ − 209 }
+ − 210 }
+ − 211 }
8
+ − 212 averages = sum(dh) / 24;
5
+ − 213 }
6
+ − 214 return(c(curr_mean_temp, averages))
5
+ − 215 }
+ − 216
35
+ − 217 get_tick_index = function(index, last_tick, ticks, month_labels) {
+ − 218 # The R code tries hard not to draw overlapping tick labels, and so
+ − 219 # will omit labels where they would abut or overlap previously drawn
+ − 220 # labels. This can result in, for example, every other tick being
+ − 221 # labelled. We'll keep track of the last tick to make sure all of
+ − 222 # the month labels are displayed, and missing ticks are restricted
+ − 223 # to Sundays which have no labels anyway.
+ − 224 if (last_tick==0) {
+ − 225 return(length(ticks)+1);
+ − 226 }
+ − 227 last_saved_tick = ticks[[length(ticks)]];
40
+ − 228 if (index-last_saved_tick<3) {
35
+ − 229 last_saved_month = month_labels[[length(month_labels)]];
+ − 230 if (last_saved_month=="") {
+ − 231 # We're safe overwriting a tick
+ − 232 # with no label (i.e., a Sunday tick).
+ − 233 return(length(ticks));
+ − 234 } else {
+ − 235 # Don't eliminate a Month label.
+ − 236 return(NULL);
+ − 237 }
+ − 238 }
+ − 239 return(length(ticks)+1);
+ − 240 }
+ − 241
38
+ − 242 get_total_days = function(is_leap_year) {
+ − 243 # Get the total number of days in the current year.
+ − 244 if (is_leap_year) {
39
+ − 245 return(366);
38
+ − 246 } else {
39
+ − 247 return(365);
38
+ − 248 }
+ − 249 }
+ − 250
39
+ − 251 get_x_axis_ticks_and_labels = function(temperature_data_frame, num_rows, start_doy_ytd, end_doy_ytd) {
35
+ − 252 # Keep track of the years to see if spanning years.
+ − 253 month_labels = list();
+ − 254 ticks = list();
+ − 255 current_month_label = NULL;
+ − 256 last_tick = 0;
+ − 257 for (i in 1:num_rows) {
39
+ − 258 if (start_doy_ytd > 1 & i==start_doy_ytd-1) {
+ − 259 # Add a tick for the end of the 30 year normnals data
+ − 260 # that was prepended to the year-to-date data.
38
+ − 261 tick_index = get_tick_index(i, last_tick, ticks, month_labels)
+ − 262 ticks[tick_index] = i;
39
+ − 263 month_labels[tick_index] = "End prepended 30 year normals";
38
+ − 264 last_tick = i;
43
+ − 265 } else if (end_doy_ytd > 0 & i==end_doy_ytd+1) {
39
+ − 266 # Add a tick for the start of the 30 year normnals data
+ − 267 # that was appended to the year-to-date data.
+ − 268 tick_index = get_tick_index(i, last_tick, ticks, month_labels)
38
+ − 269 ticks[tick_index] = i;
39
+ − 270 month_labels[tick_index] = "Start appended 30 year normals";
38
+ − 271 last_tick = i;
39
+ − 272 } else if (i==num_rows) {
38
+ − 273 # Add a tick for the last day of the year.
+ − 274 tick_index = get_tick_index(i, last_tick, ticks, month_labels)
+ − 275 ticks[tick_index] = i;
+ − 276 month_labels[tick_index] = "";
+ − 277 last_tick = i;
39
+ − 278 } else {
+ − 279 # Get the year and month from the date which
+ − 280 # has the format YYYY-MM-DD.
+ − 281 date = format(temperature_data_frame$DATE[i]);
+ − 282 # Get the month label.
+ − 283 items = strsplit(date, "-")[[1]];
+ − 284 month = items[2];
+ − 285 month_label = month.abb[as.integer(month)];
+ − 286 if (!identical(current_month_label, month_label)) {
+ − 287 # Add an x-axis tick for the month.
+ − 288 tick_index = get_tick_index(i, last_tick, ticks, month_labels)
+ − 289 ticks[tick_index] = i;
+ − 290 month_labels[tick_index] = month_label;
+ − 291 current_month_label = month_label;
+ − 292 last_tick = i;
+ − 293 }
+ − 294 tick_index = get_tick_index(i, last_tick, ticks, month_labels)
+ − 295 if (!is.null(tick_index)) {
+ − 296 # Get the day.
+ − 297 day = weekdays(as.Date(date));
+ − 298 if (day=="Sunday") {
+ − 299 # Add an x-axis tick if we're on a Sunday.
+ − 300 ticks[tick_index] = i;
+ − 301 # Add a blank month label so it is not displayed.
+ − 302 month_labels[tick_index] = "";
+ − 303 last_tick = i;
+ − 304 }
+ − 305 }
38
+ − 306 }
35
+ − 307 }
+ − 308 return(list(ticks, month_labels));
+ − 309 }
+ − 310
38
+ − 311 is_leap_year = function(date_str) {
+ − 312 # Extract the year from the date_str.
+ − 313 date = format(date_str);
+ − 314 items = strsplit(date, "-")[[1]];
+ − 315 year = as.integer(items[1]);
+ − 316 if (((year %% 4 == 0) & (year %% 100 != 0)) | (year %% 400 == 0)) {
39
+ − 317 return(TRUE);
38
+ − 318 } else {
39
+ − 319 return(FALSE);
38
+ − 320 }
+ − 321 }
+ − 322
6
+ − 323 mortality.adult = function(temperature) {
+ − 324 if (temperature < 12.7) {
8
+ − 325 mortality.probability = 0.002;
6
+ − 326 }
+ − 327 else {
8
+ − 328 mortality.probability = temperature * 0.0005 + 0.02;
6
+ − 329 }
+ − 330 return(mortality.probability)
5
+ − 331 }
+ − 332
+ − 333 mortality.egg = function(temperature) {
+ − 334 if (temperature < 12.7) {
8
+ − 335 mortality.probability = 0.8;
5
+ − 336 }
+ − 337 else {
8
+ − 338 mortality.probability = 0.8 - temperature / 40.0;
6
+ − 339 if (mortality.probability < 0) {
8
+ − 340 mortality.probability = 0.01;
5
+ − 341 }
+ − 342 }
6
+ − 343 return(mortality.probability)
5
+ − 344 }
+ − 345
+ − 346 mortality.nymph = function(temperature) {
+ − 347 if (temperature < 12.7) {
8
+ − 348 mortality.probability = 0.03;
5
+ − 349 }
+ − 350 else {
8
+ − 351 mortality.probability = temperature * 0.0008 + 0.03;
5
+ − 352 }
8
+ − 353 return(mortality.probability);
6
+ − 354 }
+ − 355
38
+ − 356 parse_input_data = function(input_ytd, input_norm, num_days_ytd) {
43
+ − 357 if (is.null(input_ytd)) {
+ − 358 # We're analysing only the 30 year normals data, so create an empty
+ − 359 # data frame for containing temperature data after it is converted
+ − 360 # from the 30 year normals format to the year-to-date format.
+ − 361 temperature_data_frame = data.frame(matrix(ncol=6, nrow=0));
+ − 362 colnames(temperature_data_frame) = c("LATITUDE", "LONGITUDE", "DATE", "DOY", "TMIN", "TMAX");
+ − 363 # Base all dates on the current date since 30 year
+ − 364 # normals data does not include any dates.
+ − 365 year = format(Sys.Date(), "%Y");
+ − 366 start_date = paste(year, "01", "01", sep="-");
+ − 367 end_date = paste(year, "12", "31", sep="-");
+ − 368 # Set invalid start and end DOY.
+ − 369 start_doy_ytd = 0;
+ − 370 end_doy_ytd = 0;
+ − 371 } else {
+ − 372 # Read the input_ytd temperature datafile into a data frame.
+ − 373 # The input_ytd data has the following 6 columns:
+ − 374 # LATITUDE, LONGITUDE, DATE, DOY, TMIN, TMAX
+ − 375 temperature_data_frame = read.csv(file=input_ytd, header=T, strip.white=TRUE, stringsAsFactors=FALSE, sep=",");
+ − 376 # Set the temperature_data_frame column names for access.
+ − 377 colnames(temperature_data_frame) = c("LATITUDE", "LONGITUDE", "DATE", "DOY", "TMIN", "TMAX");
+ − 378 # Get the start date.
+ − 379 start_date = temperature_data_frame$DATE[1];
+ − 380 end_date = temperature_data_frame$DATE[num_days_ytd];
+ − 381 # Extract the year from the start date.
+ − 382 date_str = format(start_date);
+ − 383 date_str_items = strsplit(date_str, "-")[[1]];
+ − 384 year = date_str_items[1];
+ − 385 # Save the first DOY to later check if start_date is Jan 1.
+ − 386 start_doy_ytd = as.integer(temperature_data_frame$DOY[1]);
+ − 387 end_doy_ytd = as.integer(temperature_data_frame$DOY[num_days_ytd]);
+ − 388 }
38
+ − 389 # See if we're in a leap year.
+ − 390 is_leap_year = is_leap_year(start_date);
39
+ − 391 # Get the number of days in the year.
38
+ − 392 total_days = get_total_days(is_leap_year);
+ − 393 # Read the input_norm temperature datafile into a data frame.
+ − 394 # The input_norm data has the following 10 columns:
+ − 395 # STATIONID, LATITUDE, LONGITUDE, ELEV_M, NAME, ST, MMDD, DOY, TMIN, TMAX
+ − 396 norm_data_frame = read.csv(file=input_norm, header=T, strip.white=TRUE, stringsAsFactors=FALSE, sep=",");
+ − 397 # Set the norm_data_frame column names for access.
+ − 398 colnames(norm_data_frame) = c("STATIONID", "LATITUDE","LONGITUDE", "ELEV_M", "NAME", "ST", "MMDD", "DOY", "TMIN", "TMAX");
+ − 399 # All normals data includes Feb 29 which is row 60 in
+ − 400 # the data, so delete that row if we're not in a leap year.
+ − 401 if (!is_leap_year) {
+ − 402 norm_data_frame = norm_data_frame[-c(60),];
6
+ − 403 }
43
+ − 404 if (is.null(input_ytd)) {
+ − 405 # Convert the 30 year normals data to the year-to-date format.
+ − 406 for (i in 1:total_days) {
+ − 407 temperature_data_frame[i,] = get_next_normals_row(norm_data_frame, year, is_leap_year, i);
38
+ − 408 }
43
+ − 409 } else {
+ − 410 # Merge the year-to-date data with the 30 year normals data.
+ − 411 if (start_doy_ytd > 1) {
+ − 412 # The year-to-date data starts after Jan 1, so create a
+ − 413 # temporary data frame to contain the 30 year normals data
+ − 414 # from Jan 1 to the date immediately prior to start_date.
+ − 415 tmp_data_frame = temperature_data_frame[FALSE,];
+ − 416 for (i in 1:start_doy_ytd-1) {
+ − 417 tmp_data_frame[i,] = get_next_normals_row(norm_data_frame, year, is_leap_year, i);
+ − 418 }
+ − 419 # Next merge the temporary data frame with the year-to-date data frame.
+ − 420 temperature_data_frame = rbind(tmp_data_frame, temperature_data_frame);
+ − 421 }
+ − 422 # Define the next row for the year-to-date data from the 30 year normals data.
+ − 423 first_normals_append_row = end_doy_ytd + 1;
+ − 424 # Append the 30 year normals data to the year-to-date data.
+ − 425 for (i in first_normals_append_row:total_days) {
+ − 426 temperature_data_frame[i,] = get_next_normals_row(norm_data_frame, year, is_leap_year, i);
+ − 427 }
38
+ − 428 }
+ − 429 # Add a column containing the daylight length for each day.
+ − 430 temperature_data_frame = add_daylight_length(temperature_data_frame, total_days);
41
+ − 431 return(list(temperature_data_frame, start_date, end_date, start_doy_ytd, end_doy_ytd, is_leap_year, total_days));
5
+ − 432 }
+ − 433
34
+ − 434 render_chart = function(ticks, date_labels, chart_type, plot_std_error, insect, location, latitude, start_date, end_date, days, maxval,
39
+ − 435 replications, life_stage, group, group_std_error, group2=NULL, group2_std_error=NULL, group3=NULL, group3_std_error=NULL,
+ − 436 life_stages_adult=NULL, life_stages_nymph=NULL) {
10
+ − 437 if (chart_type=="pop_size_by_life_stage") {
+ − 438 if (life_stage=="Total") {
+ − 439 title = paste(insect, ": Reps", replications, ":", life_stage, "Pop :", location, ": Lat", latitude, ":", start_date, "-", end_date, sep=" ");
+ − 440 legend_text = c("Egg", "Nymph", "Adult");
+ − 441 columns = c(4, 2, 1);
35
+ − 442 plot(days, group, main=title, type="l", ylim=c(0, maxval), axes=FALSE, lwd=2, xlab="", ylab="", cex=3, cex.lab=3, cex.axis=3, cex.main=3);
10
+ − 443 legend("topleft", legend_text, lty=c(1, 1, 1), col=columns, cex=3);
+ − 444 lines(days, group2, lwd=2, lty=1, col=2);
+ − 445 lines(days, group3, lwd=2, lty=1, col=4);
38
+ − 446 axis(side=1, at=ticks, labels=date_labels, las=2, font.axis=3, xpd=TRUE, cex=3, cex.lab=3, cex.axis=3, cex.main=3);
35
+ − 447 axis(side=2, font.axis=3, xpd=TRUE, cex=3, cex.lab=3, cex.axis=3, cex.main=3);
10
+ − 448 if (plot_std_error=="yes") {
+ − 449 # Standard error for group.
+ − 450 lines(days, group+group_std_error, lty=2);
+ − 451 lines(days, group-group_std_error, lty=2);
+ − 452 # Standard error for group2.
+ − 453 lines(days, group2+group2_std_error, col=2, lty=2);
+ − 454 lines(days, group2-group2_std_error, col=2, lty=2);
+ − 455 # Standard error for group3.
+ − 456 lines(days, group3+group3_std_error, col=4, lty=2);
+ − 457 lines(days, group3-group3_std_error, col=4, lty=2);
+ − 458 }
+ − 459 } else {
+ − 460 if (life_stage=="Egg") {
+ − 461 title = paste(insect, ": Reps", replications, ":", life_stage, "Pop :", location, ": Lat", latitude, ":", start_date, "-", end_date, sep=" ");
+ − 462 legend_text = c(life_stage);
15
+ − 463 columns = c(4);
10
+ − 464 } else if (life_stage=="Nymph") {
16
+ − 465 stage = paste(life_stages_nymph, "Nymph Pop :", sep=" ");
10
+ − 466 title = paste(insect, ": Reps", replications, ":", stage, location, ": Lat", latitude, ":", start_date, "-", end_date, sep=" ");
16
+ − 467 legend_text = c(paste(life_stages_nymph, life_stage, sep=" "));
10
+ − 468 columns = c(2);
+ − 469 } else if (life_stage=="Adult") {
+ − 470 stage = paste(life_stages_adult, "Adult Pop", sep=" ");
+ − 471 title = paste(insect, ": Reps", replications, ":", stage, location, ": Lat", latitude, ":", start_date, "-", end_date, sep=" ");
+ − 472 legend_text = c(paste(life_stages_adult, life_stage, sep=" "));
+ − 473 columns = c(1);
+ − 474 }
35
+ − 475 plot(days, group, main=title, type="l", ylim=c(0, maxval), axes=FALSE, lwd=2, xlab="", ylab="", cex=3, cex.lab=3, cex.axis=3, cex.main=3);
10
+ − 476 legend("topleft", legend_text, lty=c(1), col="black", cex=3);
38
+ − 477 axis(side=1, at=ticks, labels=date_labels, las=2, font.axis=3, xpd=TRUE, cex=3, cex.lab=3, cex.axis=3, cex.main=3);
35
+ − 478 axis(side=2, font.axis=3, xpd=TRUE, cex=3, cex.lab=3, cex.axis=3, cex.main=3);
10
+ − 479 if (plot_std_error=="yes") {
+ − 480 # Standard error for group.
+ − 481 lines(days, group+group_std_error, lty=2);
+ − 482 lines(days, group-group_std_error, lty=2);
+ − 483 }
+ − 484 }
+ − 485 } else if (chart_type=="pop_size_by_generation") {
+ − 486 if (life_stage=="Total") {
+ − 487 title_str = ": Total Pop by Gen :";
+ − 488 } else if (life_stage=="Egg") {
+ − 489 title_str = ": Egg Pop by Gen :";
+ − 490 } else if (life_stage=="Nymph") {
16
+ − 491 title_str = paste(":", life_stages_nymph, "Nymph Pop by Gen", ":", sep=" ");
10
+ − 492 } else if (life_stage=="Adult") {
+ − 493 title_str = paste(":", life_stages_adult, "Adult Pop by Gen", ":", sep=" ");
+ − 494 }
+ − 495 title = paste(insect, ": Reps", replications, title_str, location, ": Lat", latitude, ":", start_date, "-", end_date, sep=" ");
8
+ − 496 legend_text = c("P", "F1", "F2");
+ − 497 columns = c(1, 2, 4);
36
+ − 498 plot(days, group, main=title, type="l", ylim=c(0, maxval), axes=FALSE, lwd=2, xlab="", ylab="", cex=3, cex.lab=3, cex.axis=3, cex.main=3);
10
+ − 499 legend("topleft", legend_text, lty=c(1, 1, 1), col=columns, cex=3);
+ − 500 lines(days, group2, lwd=2, lty=1, col=2);
+ − 501 lines(days, group3, lwd=2, lty=1, col=4);
38
+ − 502 axis(side=1, at=ticks, labels=date_labels, las=2, font.axis=3, xpd=TRUE, cex=3, cex.lab=3, cex.axis=3, cex.main=3);
35
+ − 503 axis(side=2, font.axis=3, xpd=TRUE, cex=3, cex.lab=3, cex.axis=3, cex.main=3);
10
+ − 504 if (plot_std_error=="yes") {
+ − 505 # Standard error for group.
+ − 506 lines(days, group+group_std_error, lty=2);
+ − 507 lines(days, group-group_std_error, lty=2);
+ − 508 # Standard error for group2.
+ − 509 lines(days, group2+group2_std_error, col=2, lty=2);
+ − 510 lines(days, group2-group2_std_error, col=2, lty=2);
+ − 511 # Standard error for group3.
+ − 512 lines(days, group3+group3_std_error, col=4, lty=2);
+ − 513 lines(days, group3-group3_std_error, col=4, lty=2);
+ − 514 }
5
+ − 515 }
+ − 516 }
+ − 517
10
+ − 518 # Determine if we're plotting generations separately.
+ − 519 if (opt$plot_generations_separately=="yes") {
+ − 520 plot_generations_separately = TRUE;
+ − 521 } else {
+ − 522 plot_generations_separately = FALSE;
+ − 523 }
38
+ − 524 # Display the total number of days in the Galaxy history item blurb.
+ − 525 cat("Year-to-date number of days: ", opt$num_days_ytd, "\n");
+ − 526
39
+ − 527 # Parse the inputs.
+ − 528 data_list = parse_input_data(opt$input_ytd, opt$input_norm, opt$num_days_ytd);
+ − 529 temperature_data_frame = data_list[[1]];
+ − 530 # Information needed for plots.
+ − 531 start_date = data_list[[2]];
41
+ − 532 end_date = data_list[[3]];
+ − 533 start_doy_ytd = data_list[[4]];
+ − 534 end_doy_ytd = data_list[[5]];
+ − 535 is_leap_year = data_list[[6]];
+ − 536 total_days = data_list[[7]];
39
+ − 537 total_days_vector = c(1:total_days);
38
+ − 538
31
+ − 539 # Create copies of the temperature data for generations P, F1 and F2 if we're plotting generations separately.
+ − 540 if (plot_generations_separately) {
+ − 541 temperature_data_frame_P = data.frame(temperature_data_frame);
+ − 542 temperature_data_frame_F1 = data.frame(temperature_data_frame);
+ − 543 temperature_data_frame_F2 = data.frame(temperature_data_frame);
+ − 544 }
38
+ − 545
+ − 546 # Get the ticks date labels for plots.
39
+ − 547 ticks_and_labels = get_x_axis_ticks_and_labels(temperature_data_frame, total_days, start_doy_ytd, end_doy_ytd);
34
+ − 548 ticks = c(unlist(ticks_and_labels[1]));
+ − 549 date_labels = c(unlist(ticks_and_labels[2]));
10
+ − 550 # All latitude values are the same, so get the value for plots from the first row.
8
+ − 551 latitude = temperature_data_frame$LATITUDE[1];
38
+ − 552
20
+ − 553 # Determine the specified life stages for processing.
10
+ − 554 # Split life_stages into a list of strings for plots.
+ − 555 life_stages_str = as.character(opt$life_stages);
+ − 556 life_stages = strsplit(life_stages_str, ",")[[1]];
38
+ − 557
10
+ − 558 # Determine the data we need to generate for plotting.
+ − 559 process_eggs = FALSE;
+ − 560 process_nymphs = FALSE;
20
+ − 561 process_young_nymphs = FALSE;
+ − 562 process_old_nymphs = FALSE;
+ − 563 process_total_nymphs = FALSE;
10
+ − 564 process_adults = FALSE;
23
+ − 565 process_previttelogenic_adults = FALSE;
+ − 566 process_vittelogenic_adults = FALSE;
20
+ − 567 process_diapausing_adults = FALSE;
+ − 568 process_total_adults = FALSE;
10
+ − 569 for (life_stage in life_stages) {
+ − 570 if (life_stage=="Total") {
+ − 571 process_eggs = TRUE;
+ − 572 process_nymphs = TRUE;
+ − 573 process_adults = TRUE;
+ − 574 } else if (life_stage=="Egg") {
+ − 575 process_eggs = TRUE;
+ − 576 } else if (life_stage=="Nymph") {
+ − 577 process_nymphs = TRUE;
+ − 578 } else if (life_stage=="Adult") {
+ − 579 process_adults = TRUE;
+ − 580 }
+ − 581 }
20
+ − 582 if (process_nymphs) {
+ − 583 # Split life_stages_nymph into a list of strings for plots.
+ − 584 life_stages_nymph_str = as.character(opt$life_stages_nymph);
+ − 585 life_stages_nymph = strsplit(life_stages_nymph_str, ",")[[1]];
23
+ − 586 for (life_stage_nymph in life_stages_nymph) {
20
+ − 587 if (life_stage_nymph=="Young") {
+ − 588 process_young_nymphs = TRUE;
+ − 589 } else if (life_stage_nymph=="Old") {
+ − 590 process_old_nymphs = TRUE;
+ − 591 } else if (life_stage_nymph=="Total") {
+ − 592 process_total_nymphs = TRUE;
+ − 593 }
+ − 594 }
+ − 595 }
16
+ − 596 if (process_adults) {
+ − 597 # Split life_stages_adult into a list of strings for plots.
+ − 598 life_stages_adult_str = as.character(opt$life_stages_adult);
+ − 599 life_stages_adult = strsplit(life_stages_adult_str, ",")[[1]];
23
+ − 600 for (life_stage_adult in life_stages_adult) {
+ − 601 if (life_stage_adult=="Pre-vittelogenic") {
+ − 602 process_previttelogenic_adults = TRUE;
24
+ − 603 } else if (life_stage_adult=="Vittelogenic") {
23
+ − 604 process_vittelogenic_adults = TRUE;
20
+ − 605 } else if (life_stage_adult=="Diapausing") {
+ − 606 process_diapausing_adults = TRUE;
+ − 607 } else if (life_stage_adult=="Total") {
+ − 608 process_total_adults = TRUE;
+ − 609 }
+ − 610 }
16
+ − 611 }
6
+ − 612 # Initialize matrices.
10
+ − 613 if (process_eggs) {
38
+ − 614 Eggs.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
10
+ − 615 }
23
+ − 616 if (process_young_nymphs | process_total_nymphs) {
38
+ − 617 YoungNymphs.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
20
+ − 618 }
23
+ − 619 if (process_old_nymphs | process_total_nymphs) {
38
+ − 620 OldNymphs.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
10
+ − 621 }
23
+ − 622 if (process_previttelogenic_adults | process_total_adults) {
38
+ − 623 Previttelogenic.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
23
+ − 624 }
+ − 625 if (process_vittelogenic_adults | process_total_adults) {
38
+ − 626 Vittelogenic.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
23
+ − 627 }
+ − 628 if (process_diapausing_adults | process_total_adults) {
38
+ − 629 Diapausing.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
10
+ − 630 }
38
+ − 631 newborn.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
+ − 632 adult.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
+ − 633 death.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
10
+ − 634 if (plot_generations_separately) {
+ − 635 # P is Parental, or overwintered adults.
38
+ − 636 P.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
10
+ − 637 # F1 is the first field-produced generation.
38
+ − 638 F1.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
10
+ − 639 # F2 is the second field-produced generation.
38
+ − 640 F2.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
10
+ − 641 if (process_eggs) {
38
+ − 642 P_eggs.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
+ − 643 F1_eggs.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
+ − 644 F2_eggs.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
10
+ − 645 }
20
+ − 646 if (process_young_nymphs) {
38
+ − 647 P_young_nymphs.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
+ − 648 F1_young_nymphs.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
+ − 649 F2_young_nymphs.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
20
+ − 650 }
+ − 651 if (process_old_nymphs) {
38
+ − 652 P_old_nymphs.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
+ − 653 F1_old_nymphs.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
+ − 654 F2_old_nymphs.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
20
+ − 655 }
+ − 656 if (process_total_nymphs) {
38
+ − 657 P_total_nymphs.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
+ − 658 F1_total_nymphs.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
+ − 659 F2_total_nymphs.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
10
+ − 660 }
23
+ − 661 if (process_previttelogenic_adults) {
38
+ − 662 P_previttelogenic_adults.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
+ − 663 F1_previttelogenic_adults.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
+ − 664 F2_previttelogenic_adults.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
23
+ − 665 }
+ − 666 if (process_vittelogenic_adults) {
38
+ − 667 P_vittelogenic_adults.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
+ − 668 F1_vittelogenic_adults.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
+ − 669 F2_vittelogenic_adults.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
23
+ − 670 }
+ − 671 if (process_diapausing_adults) {
38
+ − 672 P_diapausing_adults.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
+ − 673 F1_diapausing_adults.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
+ − 674 F2_diapausing_adults.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
23
+ − 675 }
+ − 676 if (process_total_adults) {
38
+ − 677 P_total_adults.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
+ − 678 F1_total_adults.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
+ − 679 F2_total_adults.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
10
+ − 680 }
+ − 681 }
+ − 682 # Total population.
38
+ − 683 population.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
5
+ − 684
6
+ − 685 # Process replications.
18
+ − 686 for (current_replication in 1:opt$replications) {
6
+ − 687 # Start with the user-defined number of insects per replication.
8
+ − 688 num_insects = opt$insects_per_replication;
6
+ − 689 # Generation, Stage, degree-days, T, Diapause.
8
+ − 690 vector.ini = c(0, 3, 0, 0, 0);
10
+ − 691 # Replicate to create a matrix where the columns are
+ − 692 # Generation, Stage, degree-days, T, Diapause and the
+ − 693 # rows are the initial number of insects per replication.
8
+ − 694 vector.matrix = rep(vector.ini, num_insects);
10
+ − 695 # Complete transposed matrix for the population, so now
+ − 696 # the rows are Generation, Stage, degree-days, T, Diapause
8
+ − 697 vector.matrix = base::t(matrix(vector.matrix, nrow=5));
5
+ − 698 # Time series of population size.
10
+ − 699 if (process_eggs) {
38
+ − 700 Eggs = rep(0, total_days);
10
+ − 701 }
23
+ − 702 if (process_young_nymphs | process_total_nymphs) {
38
+ − 703 YoungNymphs = rep(0, total_days);
23
+ − 704 }
+ − 705 if (process_old_nymphs | process_total_nymphs) {
38
+ − 706 OldNymphs = rep(0, total_days);
10
+ − 707 }
23
+ − 708 if (process_previttelogenic_adults | process_total_adults) {
38
+ − 709 Previttelogenic = rep(0, total_days);
23
+ − 710 }
+ − 711 if (process_vittelogenic_adults | process_total_adults) {
38
+ − 712 Vittelogenic = rep(0, total_days);
23
+ − 713 }
+ − 714 if (process_diapausing_adults | process_total_adults) {
38
+ − 715 Diapausing = rep(0, total_days);
10
+ − 716 }
38
+ − 717 N.newborn = rep(0, total_days);
+ − 718 N.adult = rep(0, total_days);
+ − 719 N.death = rep(0, total_days);
+ − 720 overwintering_adult.population = rep(0, total_days);
+ − 721 first_generation.population = rep(0, total_days);
+ − 722 second_generation.population = rep(0, total_days);
10
+ − 723 if (plot_generations_separately) {
+ − 724 # P is Parental, or overwintered adults.
+ − 725 # F1 is the first field-produced generation.
+ − 726 # F2 is the second field-produced generation.
+ − 727 if (process_eggs) {
38
+ − 728 P.egg = rep(0, total_days);
+ − 729 F1.egg = rep(0, total_days);
+ − 730 F2.egg = rep(0, total_days);
10
+ − 731 }
20
+ − 732 if (process_young_nymphs) {
38
+ − 733 P.young_nymph = rep(0, total_days);
+ − 734 F1.young_nymph = rep(0, total_days);
+ − 735 F2.young_nymph = rep(0, total_days);
20
+ − 736 }
+ − 737 if (process_old_nymphs) {
38
+ − 738 P.old_nymph = rep(0, total_days);
+ − 739 F1.old_nymph = rep(0, total_days);
+ − 740 F2.old_nymph = rep(0, total_days);
20
+ − 741 }
+ − 742 if (process_total_nymphs) {
38
+ − 743 P.total_nymph = rep(0, total_days);
+ − 744 F1.total_nymph = rep(0, total_days);
+ − 745 F2.total_nymph = rep(0, total_days);
10
+ − 746 }
23
+ − 747 if (process_previttelogenic_adults) {
38
+ − 748 P.previttelogenic_adult = rep(0, total_days);
+ − 749 F1.previttelogenic_adult = rep(0, total_days);
+ − 750 F2.previttelogenic_adult = rep(0, total_days);
23
+ − 751 }
+ − 752 if (process_vittelogenic_adults) {
38
+ − 753 P.vittelogenic_adult = rep(0, total_days);
+ − 754 F1.vittelogenic_adult = rep(0, total_days);
+ − 755 F2.vittelogenic_adult = rep(0, total_days);
23
+ − 756 }
+ − 757 if (process_diapausing_adults) {
38
+ − 758 P.diapausing_adult = rep(0, total_days);
+ − 759 F1.diapausing_adult = rep(0, total_days);
+ − 760 F2.diapausing_adult = rep(0, total_days);
23
+ − 761 }
+ − 762 if (process_total_adults) {
38
+ − 763 P.total_adult = rep(0, total_days);
+ − 764 F1.total_adult = rep(0, total_days);
+ − 765 F2.total_adult = rep(0, total_days);
10
+ − 766 }
+ − 767 }
8
+ − 768 total.population = NULL;
38
+ − 769 averages.day = rep(0, total_days);
+ − 770 # All the days included in the input_ytd temperature dataset.
+ − 771 for (row in 1:total_days) {
5
+ − 772 # Get the integer day of the year for the current row.
8
+ − 773 doy = temperature_data_frame$DOY[row];
5
+ − 774 # Photoperiod in the day.
8
+ − 775 photoperiod = temperature_data_frame$DAYLEN[row];
38
+ − 776 temp.profile = get_temperature_at_hour(latitude, temperature_data_frame, row, total_days);
8
+ − 777 mean.temp = temp.profile[1];
+ − 778 averages.temp = temp.profile[2];
+ − 779 averages.day[row] = averages.temp;
5
+ − 780 # Trash bin for death.
8
+ − 781 death.vector = NULL;
5
+ − 782 # Newborn.
8
+ − 783 birth.vector = NULL;
5
+ − 784 # All individuals.
6
+ − 785 for (i in 1:num_insects) {
+ − 786 # Individual record.
8
+ − 787 vector.individual = vector.matrix[i,];
6
+ − 788 # Adjustment for late season mortality rate (still alive?).
5
+ − 789 if (latitude < 40.0) {
8
+ − 790 post.mortality = 1;
+ − 791 day.kill = 300;
5
+ − 792 }
+ − 793 else {
8
+ − 794 post.mortality = 2;
+ − 795 day.kill = 250;
5
+ − 796 }
6
+ − 797 if (vector.individual[2] == 0) {
5
+ − 798 # Egg.
8
+ − 799 death.probability = opt$egg_mortality * mortality.egg(mean.temp);
5
+ − 800 }
6
+ − 801 else if (vector.individual[2] == 1 | vector.individual[2] == 2) {
18
+ − 802 # Nymph.
8
+ − 803 death.probability = opt$nymph_mortality * mortality.nymph(mean.temp);
5
+ − 804 }
6
+ − 805 else if (vector.individual[2] == 3 | vector.individual[2] == 4 | vector.individual[2] == 5) {
+ − 806 # Adult.
5
+ − 807 if (doy < day.kill) {
8
+ − 808 death.probability = opt$adult_mortality * mortality.adult(mean.temp);
5
+ − 809 }
+ − 810 else {
+ − 811 # Increase adult mortality after fall equinox.
8
+ − 812 death.probability = opt$adult_mortality * post.mortality * mortality.adult(mean.temp);
5
+ − 813 }
+ − 814 }
6
+ − 815 # Dependent on temperature and life stage?
8
+ − 816 u.d = runif(1);
6
+ − 817 if (u.d < death.probability) {
8
+ − 818 death.vector = c(death.vector, i);
6
+ − 819 }
5
+ − 820 else {
6
+ − 821 # End of diapause.
+ − 822 if (vector.individual[1] == 0 && vector.individual[2] == 3) {
27
+ − 823 # Overwintering adult (pre-vittelogenic).
6
+ − 824 if (photoperiod > opt$photoperiod && vector.individual[3] > 68 && doy < 180) {
5
+ − 825 # Add 68C to become fully reproductively matured.
+ − 826 # Transfer to vittelogenic.
8
+ − 827 vector.individual = c(0, 4, 0, 0, 0);
+ − 828 vector.matrix[i,] = vector.individual;
5
+ − 829 }
+ − 830 else {
27
+ − 831 # Add average temperature for current day.
8
+ − 832 vector.individual[3] = vector.individual[3] + averages.temp;
5
+ − 833 # Add 1 day in current stage.
8
+ − 834 vector.individual[4] = vector.individual[4] + 1;
+ − 835 vector.matrix[i,] = vector.individual;
5
+ − 836 }
+ − 837 }
6
+ − 838 if (vector.individual[1] != 0 && vector.individual[2] == 3) {
27
+ − 839 # Not overwintering adult (pre-vittelogenic).
8
+ − 840 current.gen = vector.individual[1];
6
+ − 841 if (vector.individual[3] > 68) {
5
+ − 842 # Add 68C to become fully reproductively matured.
+ − 843 # Transfer to vittelogenic.
8
+ − 844 vector.individual = c(current.gen, 4, 0, 0, 0);
+ − 845 vector.matrix[i,] = vector.individual;
5
+ − 846 }
+ − 847 else {
6
+ − 848 # Add average temperature for current day.
8
+ − 849 vector.individual[3] = vector.individual[3] + averages.temp;
5
+ − 850 # Add 1 day in current stage.
8
+ − 851 vector.individual[4] = vector.individual[4] + 1;
+ − 852 vector.matrix[i,] = vector.individual;
5
+ − 853 }
+ − 854 }
6
+ − 855 # Oviposition -- where population dynamics comes from.
+ − 856 if (vector.individual[2] == 4 && vector.individual[1] == 0 && mean.temp > 10) {
5
+ − 857 # Vittelogenic stage, overwintering generation.
6
+ − 858 if (vector.individual[4] == 0) {
5
+ − 859 # Just turned in vittelogenic stage.
8
+ − 860 num_insects.birth = round(runif(1, 2 + opt$min_clutch_size, 8 + opt$max_clutch_size));
5
+ − 861 }
+ − 862 else {
+ − 863 # Daily probability of birth.
8
+ − 864 p.birth = opt$oviposition * 0.01;
+ − 865 u1 = runif(1);
5
+ − 866 if (u1 < p.birth) {
8
+ − 867 num_insects.birth = round(runif(1, 2, 8));
5
+ − 868 }
+ − 869 }
6
+ − 870 # Add average temperature for current day.
8
+ − 871 vector.individual[3] = vector.individual[3] + averages.temp;
5
+ − 872 # Add 1 day in current stage.
8
+ − 873 vector.individual[4] = vector.individual[4] + 1;
+ − 874 vector.matrix[i,] = vector.individual;
6
+ − 875 if (num_insects.birth > 0) {
5
+ − 876 # Add new birth -- might be in different generations.
8
+ − 877 new.gen = vector.individual[1] + 1;
5
+ − 878 # Egg profile.
8
+ − 879 new.individual = c(new.gen, 0, 0, 0, 0);
+ − 880 new.vector = rep(new.individual, num_insects.birth);
5
+ − 881 # Update batch of egg profile.
8
+ − 882 new.vector = t(matrix(new.vector, nrow=5));
5
+ − 883 # Group with total eggs laid in that day.
8
+ − 884 birth.vector = rbind(birth.vector, new.vector);
5
+ − 885 }
+ − 886 }
6
+ − 887 # Oviposition -- for generation 1.
+ − 888 if (vector.individual[2] == 4 && vector.individual[1] == 1 && mean.temp > 12.5 && doy < 222) {
5
+ − 889 # Vittelogenic stage, 1st generation
6
+ − 890 if (vector.individual[4] == 0) {
5
+ − 891 # Just turned in vittelogenic stage.
8
+ − 892 num_insects.birth = round(runif(1, 2+opt$min_clutch_size, 8+opt$max_clutch_size));
5
+ − 893 }
+ − 894 else {
+ − 895 # Daily probability of birth.
8
+ − 896 p.birth = opt$oviposition * 0.01;
+ − 897 u1 = runif(1);
5
+ − 898 if (u1 < p.birth) {
8
+ − 899 num_insects.birth = round(runif(1, 2, 8));
5
+ − 900 }
+ − 901 }
6
+ − 902 # Add average temperature for current day.
8
+ − 903 vector.individual[3] = vector.individual[3] + averages.temp;
5
+ − 904 # Add 1 day in current stage.
8
+ − 905 vector.individual[4] = vector.individual[4] + 1;
+ − 906 vector.matrix[i,] = vector.individual;
6
+ − 907 if (num_insects.birth > 0) {
5
+ − 908 # Add new birth -- might be in different generations.
8
+ − 909 new.gen = vector.individual[1] + 1;
5
+ − 910 # Egg profile.
8
+ − 911 new.individual = c(new.gen, 0, 0, 0, 0);
+ − 912 new.vector = rep(new.individual, num_insects.birth);
5
+ − 913 # Update batch of egg profile.
8
+ − 914 new.vector = t(matrix(new.vector, nrow=5));
5
+ − 915 # Group with total eggs laid in that day.
8
+ − 916 birth.vector = rbind(birth.vector, new.vector);
5
+ − 917 }
+ − 918 }
6
+ − 919 # Egg to young nymph.
+ − 920 if (vector.individual[2] == 0) {
+ − 921 # Add average temperature for current day.
8
+ − 922 vector.individual[3] = vector.individual[3] + averages.temp;
6
+ − 923 if (vector.individual[3] >= (68+opt$young_nymph_accumulation)) {
+ − 924 # From egg to young nymph, degree-days requirement met.
8
+ − 925 current.gen = vector.individual[1];
5
+ − 926 # Transfer to young nymph stage.
8
+ − 927 vector.individual = c(current.gen, 1, 0, 0, 0);
5
+ − 928 }
+ − 929 else {
+ − 930 # Add 1 day in current stage.
8
+ − 931 vector.individual[4] = vector.individual[4] + 1;
5
+ − 932 }
8
+ − 933 vector.matrix[i,] = vector.individual;
5
+ − 934 }
6
+ − 935 # Young nymph to old nymph.
+ − 936 if (vector.individual[2] == 1) {
+ − 937 # Add average temperature for current day.
8
+ − 938 vector.individual[3] = vector.individual[3] + averages.temp;
6
+ − 939 if (vector.individual[3] >= (250+opt$old_nymph_accumulation)) {
+ − 940 # From young to old nymph, degree_days requirement met.
8
+ − 941 current.gen = vector.individual[1];
5
+ − 942 # Transfer to old nym stage.
8
+ − 943 vector.individual = c(current.gen, 2, 0, 0, 0);
5
+ − 944 if (photoperiod < opt$photoperiod && doy > 180) {
8
+ − 945 vector.individual[5] = 1;
5
+ − 946 } # Prepare for diapausing.
+ − 947 }
+ − 948 else {
+ − 949 # Add 1 day in current stage.
8
+ − 950 vector.individual[4] = vector.individual[4] + 1;
5
+ − 951 }
8
+ − 952 vector.matrix[i,] = vector.individual;
6
+ − 953 }
27
+ − 954 # Old nymph to adult: pre-vittelogenic or diapausing?
6
+ − 955 if (vector.individual[2] == 2) {
+ − 956 # Add average temperature for current day.
8
+ − 957 vector.individual[3] = vector.individual[3] + averages.temp;
6
+ − 958 if (vector.individual[3] >= (200+opt$adult_accumulation)) {
+ − 959 # From old to adult, degree_days requirement met.
8
+ − 960 current.gen = vector.individual[1];
6
+ − 961 if (vector.individual[5] == 0) {
+ − 962 # Previttelogenic.
8
+ − 963 vector.individual = c(current.gen, 3, 0, 0, 0);
5
+ − 964 }
+ − 965 else {
+ − 966 # Diapausing.
8
+ − 967 vector.individual = c(current.gen, 5, 0, 0, 1);
5
+ − 968 }
+ − 969 }
+ − 970 else {
+ − 971 # Add 1 day in current stage.
8
+ − 972 vector.individual[4] = vector.individual[4] + 1;
5
+ − 973 }
8
+ − 974 vector.matrix[i,] = vector.individual;
5
+ − 975 }
6
+ − 976 # Growing of diapausing adult (unimportant, but still necessary).
+ − 977 if (vector.individual[2] == 5) {
8
+ − 978 vector.individual[3] = vector.individual[3] + averages.temp;
+ − 979 vector.individual[4] = vector.individual[4] + 1;
+ − 980 vector.matrix[i,] = vector.individual;
5
+ − 981 }
+ − 982 } # Else if it is still alive.
+ − 983 } # End of the individual bug loop.
6
+ − 984
+ − 985 # Number of deaths.
8
+ − 986 num_insects.death = length(death.vector);
6
+ − 987 if (num_insects.death > 0) {
+ − 988 # Remove record of dead.
8
+ − 989 vector.matrix = vector.matrix[-death.vector,];
5
+ − 990 }
6
+ − 991 # Number of births.
8
+ − 992 num_insects.newborn = length(birth.vector[,1]);
+ − 993 vector.matrix = rbind(vector.matrix, birth.vector);
5
+ − 994 # Update population size for the next day.
8
+ − 995 num_insects = num_insects - num_insects.death + num_insects.newborn;
5
+ − 996
10
+ − 997 # Aggregate results by day. Due to multiple transpose calls
+ − 998 # on vector.matrix above, the columns of vector.matrix
+ − 999 # are now Generation, Stage, degree-days, T, Diapause,
+ − 1000 if (process_eggs) {
+ − 1001 # For egg population size, column 2 (Stage), must be 0.
+ − 1002 Eggs[row] = sum(vector.matrix[,2]==0);
+ − 1003 }
23
+ − 1004 if (process_young_nymphs | process_total_nymphs) {
10
+ − 1005 # For young nymph population size, column 2 (Stage) must be 1.
+ − 1006 YoungNymphs[row] = sum(vector.matrix[,2]==1);
20
+ − 1007 }
23
+ − 1008 if (process_old_nymphs | process_total_nymphs) {
10
+ − 1009 # For old nymph population size, column 2 (Stage) must be 2.
+ − 1010 OldNymphs[row] = sum(vector.matrix[,2]==2);
+ − 1011 }
23
+ − 1012 if (process_previttelogenic_adults | process_total_adults) {
+ − 1013 # For pre-vittelogenic population size, column 2 (Stage) must be 3.
+ − 1014 Previttelogenic[row] = sum(vector.matrix[,2]==3);
+ − 1015 }
+ − 1016 if (process_vittelogenic_adults | process_total_adults) {
+ − 1017 # For vittelogenic population size, column 2 (Stage) must be 4.
24
+ − 1018 Vittelogenic[row] = sum(vector.matrix[,2]==4);
23
+ − 1019 }
+ − 1020 if (process_diapausing_adults | process_total_adults) {
10
+ − 1021 # For diapausing population size, column 2 (Stage) must be 5.
+ − 1022 Diapausing[row] = sum(vector.matrix[,2]==5);
+ − 1023 }
5
+ − 1024
6
+ − 1025 # Newborn population size.
8
+ − 1026 N.newborn[row] = num_insects.newborn;
6
+ − 1027 # Adult population size.
8
+ − 1028 N.adult[row] = sum(vector.matrix[,2]==3) + sum(vector.matrix[,2]==4) + sum(vector.matrix[,2]==5);
6
+ − 1029 # Dead population size.
8
+ − 1030 N.death[row] = num_insects.death;
6
+ − 1031
8
+ − 1032 total.population = c(total.population, num_insects);
6
+ − 1033
10
+ − 1034 # For overwintering adult (P) population
+ − 1035 # size, column 1 (Generation) must be 0.
8
+ − 1036 overwintering_adult.population[row] = sum(vector.matrix[,1]==0);
10
+ − 1037 # For first field generation (F1) population
+ − 1038 # size, column 1 (Generation) must be 1.
8
+ − 1039 first_generation.population[row] = sum(vector.matrix[,1]==1);
10
+ − 1040 # For second field generation (F2) population
+ − 1041 # size, column 1 (Generation) must be 2.
8
+ − 1042 second_generation.population[row] = sum(vector.matrix[,1]==2);
5
+ − 1043
10
+ − 1044 if (plot_generations_separately) {
+ − 1045 if (process_eggs) {
18
+ − 1046 # For egg life stage of generation P population size,
10
+ − 1047 # column 1 (generation) is 0 and column 2 (Stage) is 0.
+ − 1048 P.egg[row] = sum(vector.matrix[,1]==0 & vector.matrix[,2]==0);
+ − 1049 # For egg life stage of generation F1 population size,
+ − 1050 # column 1 (generation) is 1 and column 2 (Stage) is 0.
+ − 1051 F1.egg[row] = sum(vector.matrix[,1]==1 & vector.matrix[,2]==0);
+ − 1052 # For egg life stage of generation F2 population size,
+ − 1053 # column 1 (generation) is 2 and column 2 (Stage) is 0.
+ − 1054 F2.egg[row] = sum(vector.matrix[,1]==2 & vector.matrix[,2]==0);
+ − 1055 }
20
+ − 1056 if (process_young_nymphs) {
+ − 1057 # For young nymph life stage of generation P population
+ − 1058 # size, the following combination is required:
+ − 1059 # - column 1 (Generation) is 0 and column 2 (Stage) is 1 (Young nymph)
+ − 1060 P.young_nymph[row] = sum(vector.matrix[,1]==0 & vector.matrix[,2]==1);
+ − 1061 # For young nymph life stage of generation F1 population
+ − 1062 # size, the following combination is required:
+ − 1063 # - column 1 (Generation) is 1 and column 2 (Stage) is 1 (Young nymph)
+ − 1064 F1.young_nymph[row] = sum(vector.matrix[,1]==1 & vector.matrix[,2]==1);
+ − 1065 # For young nymph life stage of generation F2 population
+ − 1066 # size, the following combination is required:
+ − 1067 # - column 1 (Generation) is 2 and column 2 (Stage) is 1 (Young nymph)
+ − 1068 F2.young_nymph[row] = sum(vector.matrix[,1]==2 & vector.matrix[,2]==1);
+ − 1069 }
+ − 1070 if (process_old_nymphs) {
+ − 1071 # For old nymph life stage of generation P population
+ − 1072 # size, the following combination is required:
+ − 1073 # - column 1 (Generation) is 0 and column 2 (Stage) is 2 (Old nymph)
+ − 1074 P.old_nymph[row] = sum(vector.matrix[,1]==0 & vector.matrix[,2]==2);
+ − 1075 # For old nymph life stage of generation F1 population
+ − 1076 # size, the following combination is required:
+ − 1077 # - column 1 (Generation) is 1 and column 2 (Stage) is 2 (Old nymph)
+ − 1078 F1.old_nymph[row] = sum(vector.matrix[,1]==1 & vector.matrix[,2]==2);
+ − 1079 # For old nymph life stage of generation F2 population
+ − 1080 # size, the following combination is required:
+ − 1081 # - column 1 (Generation) is 2 and column 2 (Stage) is 2 (Old nymph)
+ − 1082 F2.old_nymph[row] = sum(vector.matrix[,1]==2 & vector.matrix[,2]==2);
+ − 1083 }
+ − 1084 if (process_total_nymphs) {
+ − 1085 # For total nymph life stage of generation P population
10
+ − 1086 # size, one of the following combinations is required:
+ − 1087 # - column 1 (Generation) is 0 and column 2 (Stage) is 1 (Young nymph)
+ − 1088 # - column 1 (Generation) is 0 and column 2 (Stage) is 2 (Old nymph)
20
+ − 1089 P.total_nymph[row] = sum((vector.matrix[,1]==0 & vector.matrix[,2]==1) | (vector.matrix[,1]==0 & vector.matrix[,2]==2));
+ − 1090 # For total nymph life stage of generation F1 population
10
+ − 1091 # size, one of the following combinations is required:
+ − 1092 # - column 1 (Generation) is 1 and column 2 (Stage) is 1 (Young nymph)
+ − 1093 # - column 1 (Generation) is 1 and column 2 (Stage) is 2 (Old nymph)
20
+ − 1094 F1.total_nymph[row] = sum((vector.matrix[,1]==1 & vector.matrix[,2]==1) | (vector.matrix[,1]==1 & vector.matrix[,2]==2));
+ − 1095 # For total nymph life stage of generation F2 population
10
+ − 1096 # size, one of the following combinations is required:
+ − 1097 # - column 1 (Generation) is 2 and column 2 (Stage) is 1 (Young nymph)
+ − 1098 # - column 1 (Generation) is 2 and column 2 (Stage) is 2 (Old nymph)
20
+ − 1099 F2.total_nymph[row] = sum((vector.matrix[,1]==2 & vector.matrix[,2]==1) | (vector.matrix[,1]==2 & vector.matrix[,2]==2));
10
+ − 1100 }
23
+ − 1101 if (process_previttelogenic_adults) {
+ − 1102 # For previttelogenic adult life stage of generation P population
+ − 1103 # size, the following combination is required:
+ − 1104 # - column 1 (Generation) is 0 and column 2 (Stage) is 3 (Pre-vittelogenic)
+ − 1105 P.previttelogenic_adult[row] = sum(vector.matrix[,1]==0 & vector.matrix[,2]==3);
+ − 1106 # For previttelogenic adult life stage of generation F1 population
+ − 1107 # size, the following combination is required:
+ − 1108 # - column 1 (Generation) is 1 and column 2 (Stage) is 3 (Pre-vittelogenic)
+ − 1109 F1.previttelogenic_adult[row] = sum(vector.matrix[,1]==1 & vector.matrix[,2]==3);
+ − 1110 # For previttelogenic adult life stage of generation F2 population
+ − 1111 # size, the following combination is required:
+ − 1112 # - column 1 (Generation) is 2 and column 2 (Stage) is 3 (Pre-vittelogenic)
+ − 1113 F2.previttelogenic_adult[row] = sum(vector.matrix[,1]==2 & vector.matrix[,2]==3);
+ − 1114 }
+ − 1115 if (process_vittelogenic_adults) {
+ − 1116 # For vittelogenic adult life stage of generation P population
+ − 1117 # size, the following combination is required:
24
+ − 1118 # - column 1 (Generation) is 0 and column 2 (Stage) is 4 (Vittelogenic)
23
+ − 1119 P.vittelogenic_adult[row] = sum(vector.matrix[,1]==0 & vector.matrix[,2]==4);
+ − 1120 # For vittelogenic adult life stage of generation F1 population
+ − 1121 # size, the following combination is required:
24
+ − 1122 # - column 1 (Generation) is 1 and column 2 (Stage) is 4 (Vittelogenic)
23
+ − 1123 F1.vittelogenic_adult[row] = sum(vector.matrix[,1]==1 & vector.matrix[,2]==4);
+ − 1124 # For vittelogenic adult life stage of generation F2 population
+ − 1125 # size, the following combination is required:
24
+ − 1126 # - column 1 (Generation) is 2 and column 2 (Stage) is 4 (Vittelogenic)
23
+ − 1127 F2.vittelogenic_adult[row] = sum(vector.matrix[,1]==2 & vector.matrix[,2]==4);
+ − 1128 }
+ − 1129 if (process_diapausing_adults) {
+ − 1130 # For diapausing adult life stage of generation P population
+ − 1131 # size, the following combination is required:
10
+ − 1132 # - column 1 (Generation) is 0 and column 2 (Stage) is 5 (Diapausing)
23
+ − 1133 P.diapausing_adult[row] = sum(vector.matrix[,1]==0 & vector.matrix[,2]==5);
+ − 1134 # For diapausing adult life stage of generation F1 population
+ − 1135 # size, the following combination is required:
+ − 1136 # - column 1 (Generation) is 1 and column 2 (Stage) is 5 (Diapausing)
+ − 1137 F1.diapausing_adult[row] = sum(vector.matrix[,1]==1 & vector.matrix[,2]==5);
+ − 1138 # For diapausing adult life stage of generation F2 population
+ − 1139 # size, the following combination is required:
+ − 1140 # - column 1 (Generation) is 2 and column 2 (Stage) is 5 (Diapausing)
+ − 1141 F2.diapausing_adult[row] = sum(vector.matrix[,1]==2 & vector.matrix[,2]==5);
+ − 1142 }
+ − 1143 if (process_total_adults) {
+ − 1144 # For total adult life stage of generation P population
10
+ − 1145 # size, one of the following combinations is required:
23
+ − 1146 # - column 1 (Generation) is 0 and column 2 (Stage) is 3 (Pre-vittelogenic)
24
+ − 1147 # - column 1 (Generation) is 0 and column 2 (Stage) is 4 (Vittelogenic)
23
+ − 1148 # - column 1 (Generation) is 0 and column 2 (Stage) is 5 (Diapausing)
+ − 1149 P.total_adult[row] = sum((vector.matrix[,1]==0 & vector.matrix[,2]==3) | (vector.matrix[,1]==0 & vector.matrix[,2]==4) | (vector.matrix[,1]==0 & vector.matrix[,2]==5));
+ − 1150 # For total adult life stage of generation F1 population
+ − 1151 # size, one of the following combinations is required:
+ − 1152 # - column 1 (Generation) is 1 and column 2 (Stage) is 3 (Pre-vittelogenic)
24
+ − 1153 # - column 1 (Generation) is 1 and column 2 (Stage) is 4 (Vittelogenic)
10
+ − 1154 # - column 1 (Generation) is 1 and column 2 (Stage) is 5 (Diapausing)
23
+ − 1155 F1.total_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));
+ − 1156 # For total adult life stage of generation F2 population
10
+ − 1157 # size, one of the following combinations is required:
23
+ − 1158 # - column 1 (Generation) is 2 and column 2 (Stage) is 3 (Pre-vittelogenic)
24
+ − 1159 # - column 1 (Generation) is 2 and column 2 (Stage) is 4 (Vittelogenic)
10
+ − 1160 # - column 1 (Generation) is 2 and column 2 (Stage) is 5 (Diapausing)
23
+ − 1161 F2.total_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));
10
+ − 1162 }
+ − 1163 }
38
+ − 1164 } # End of days specified in the input_ytd temperature data.
5
+ − 1165
8
+ − 1166 averages.cum = cumsum(averages.day);
5
+ − 1167
6
+ − 1168 # Define the output values.
10
+ − 1169 if (process_eggs) {
18
+ − 1170 Eggs.replications[,current_replication] = Eggs;
10
+ − 1171 }
23
+ − 1172 if (process_young_nymphs | process_total_nymphs) {
18
+ − 1173 YoungNymphs.replications[,current_replication] = YoungNymphs;
20
+ − 1174 }
23
+ − 1175 if (process_old_nymphs | process_total_nymphs) {
18
+ − 1176 OldNymphs.replications[,current_replication] = OldNymphs;
10
+ − 1177 }
23
+ − 1178 if (process_previttelogenic_adults | process_total_adults) {
+ − 1179 Previttelogenic.replications[,current_replication] = Previttelogenic;
+ − 1180 }
+ − 1181 if (process_vittelogenic_adults | process_total_adults) {
24
+ − 1182 Vittelogenic.replications[,current_replication] = Vittelogenic;
23
+ − 1183 }
+ − 1184 if (process_diapausing_adults | process_total_adults) {
18
+ − 1185 Diapausing.replications[,current_replication] = Diapausing;
10
+ − 1186 }
18
+ − 1187 newborn.replications[,current_replication] = N.newborn;
+ − 1188 adult.replications[,current_replication] = N.adult;
+ − 1189 death.replications[,current_replication] = N.death;
10
+ − 1190 if (plot_generations_separately) {
+ − 1191 # P is Parental, or overwintered adults.
18
+ − 1192 P.replications[,current_replication] = overwintering_adult.population;
10
+ − 1193 # F1 is the first field-produced generation.
18
+ − 1194 F1.replications[,current_replication] = first_generation.population;
10
+ − 1195 # F2 is the second field-produced generation.
18
+ − 1196 F2.replications[,current_replication] = second_generation.population;
10
+ − 1197 if (process_eggs) {
18
+ − 1198 P_eggs.replications[,current_replication] = P.egg;
+ − 1199 F1_eggs.replications[,current_replication] = F1.egg;
+ − 1200 F2_eggs.replications[,current_replication] = F2.egg;
10
+ − 1201 }
20
+ − 1202 if (process_young_nymphs) {
+ − 1203 P_young_nymphs.replications[,current_replication] = P.young_nymph;
+ − 1204 F1_young_nymphs.replications[,current_replication] = F1.young_nymph;
+ − 1205 F2_young_nymphs.replications[,current_replication] = F2.young_nymph;
+ − 1206 }
+ − 1207 if (process_old_nymphs) {
+ − 1208 P_old_nymphs.replications[,current_replication] = P.old_nymph;
+ − 1209 F1_old_nymphs.replications[,current_replication] = F1.old_nymph;
+ − 1210 F2_old_nymphs.replications[,current_replication] = F2.old_nymph;
+ − 1211 }
+ − 1212 if (process_total_nymphs) {
+ − 1213 P_total_nymphs.replications[,current_replication] = P.total_nymph;
+ − 1214 F1_total_nymphs.replications[,current_replication] = F1.total_nymph;
+ − 1215 F2_total_nymphs.replications[,current_replication] = F2.total_nymph;
10
+ − 1216 }
23
+ − 1217 if (process_previttelogenic_adults) {
+ − 1218 P_previttelogenic_adults.replications[,current_replication] = P.previttelogenic_adult;
+ − 1219 F1_previttelogenic_adults.replications[,current_replication] = F1.previttelogenic_adult;
+ − 1220 F2_previttelogenic_adults.replications[,current_replication] = F2.previttelogenic_adult;
+ − 1221 }
+ − 1222 if (process_vittelogenic_adults) {
+ − 1223 P_vittelogenic_adults.replications[,current_replication] = P.vittelogenic_adult;
+ − 1224 F1_vittelogenic_adults.replications[,current_replication] = F1.vittelogenic_adult;
+ − 1225 F2_vittelogenic_adults.replications[,current_replication] = F2.vittelogenic_adult;
+ − 1226 }
+ − 1227 if (process_diapausing_adults) {
+ − 1228 P_diapausing_adults.replications[,current_replication] = P.diapausing_adult;
+ − 1229 F1_diapausing_adults.replications[,current_replication] = F1.diapausing_adult;
+ − 1230 F2_diapausing_adults.replications[,current_replication] = F2.diapausing_adult;
+ − 1231 }
+ − 1232 if (process_total_adults) {
+ − 1233 P_total_adults.replications[,current_replication] = P.total_adult;
+ − 1234 F1_total_adults.replications[,current_replication] = F1.total_adult;
+ − 1235 F2_total_adults.replications[,current_replication] = F2.total_adult;
10
+ − 1236 }
+ − 1237 }
18
+ − 1238 population.replications[,current_replication] = total.population;
+ − 1239 # End processing replications.
5
+ − 1240 }
+ − 1241
10
+ − 1242 if (process_eggs) {
+ − 1243 # Mean value for eggs.
+ − 1244 eggs = apply(Eggs.replications, 1, mean);
27
+ − 1245 temperature_data_frame = append_vector(temperature_data_frame, eggs, "EGG");
10
+ − 1246 # Standard error for eggs.
+ − 1247 eggs.std_error = apply(Eggs.replications, 1, sd) / sqrt(opt$replications);
27
+ − 1248 temperature_data_frame = append_vector(temperature_data_frame, eggs.std_error, "EGGSE");
10
+ − 1249 }
+ − 1250 if (process_nymphs) {
+ − 1251 # Calculate nymph populations for selected life stage.
16
+ − 1252 for (life_stage_nymph in life_stages_nymph) {
28
+ − 1253 if (life_stage_nymph=="Young") {
16
+ − 1254 # Mean value for young nymphs.
+ − 1255 young_nymphs = apply(YoungNymphs.replications, 1, mean);
27
+ − 1256 temperature_data_frame = append_vector(temperature_data_frame, young_nymphs, "YOUNGNYMPH");
16
+ − 1257 # Standard error for young nymphs.
+ − 1258 young_nymphs.std_error = apply(YoungNymphs.replications / sqrt(opt$replications), 1, sd);
27
+ − 1259 temperature_data_frame = append_vector(temperature_data_frame, young_nymphs.std_error, "YOUNGNYMPHSE");
18
+ − 1260 } else if (life_stage_nymph=="Old") {
16
+ − 1261 # Mean value for old nymphs.
+ − 1262 old_nymphs = apply(OldNymphs.replications, 1, mean);
27
+ − 1263 temperature_data_frame = append_vector(temperature_data_frame, old_nymphs, "OLDNYMPH");
16
+ − 1264 # Standard error for old nymphs.
+ − 1265 old_nymphs.std_error = apply(OldNymphs.replications / sqrt(opt$replications), 1, sd);
27
+ − 1266 temperature_data_frame = append_vector(temperature_data_frame, old_nymphs.std_error, "OLDNYMPHSE");
28
+ − 1267 } else if (life_stage_nymph=="Total") {
+ − 1268 # Mean value for all nymphs.
+ − 1269 total_nymphs = apply((YoungNymphs.replications+OldNymphs.replications), 1, mean);
+ − 1270 temperature_data_frame = append_vector(temperature_data_frame, total_nymphs, "TOTALNYMPH");
+ − 1271 # Standard error for all nymphs.
+ − 1272 total_nymphs.std_error = apply((YoungNymphs.replications+OldNymphs.replications) / sqrt(opt$replications), 1, sd);
+ − 1273 temperature_data_frame = append_vector(temperature_data_frame, total_nymphs.std_error, "TOTALNYMPHSE");
16
+ − 1274 }
10
+ − 1275 }
+ − 1276 }
+ − 1277 if (process_adults) {
+ − 1278 # Calculate adult populations for selected life stage.
16
+ − 1279 for (life_stage_adult in life_stages_adult) {
28
+ − 1280 if (life_stage_adult == "Pre-vittelogenic") {
23
+ − 1281 # Mean value for previttelogenic adults.
+ − 1282 previttelogenic_adults = apply(Previttelogenic.replications, 1, mean);
27
+ − 1283 temperature_data_frame = append_vector(temperature_data_frame, previttelogenic_adults, "PRE-VITADULT");
23
+ − 1284 # Standard error for previttelogenic adults.
+ − 1285 previttelogenic_adults.std_error = apply(Previttelogenic.replications, 1, sd) / sqrt(opt$replications);
27
+ − 1286 temperature_data_frame = append_vector(temperature_data_frame, previttelogenic_adults.std_error, "PRE-VITADULTSE");
18
+ − 1287 } else if (life_stage_adult == "Vittelogenic") {
23
+ − 1288 # Mean value for vittelogenic adults.
24
+ − 1289 vittelogenic_adults = apply(Vittelogenic.replications, 1, mean);
27
+ − 1290 temperature_data_frame = append_vector(temperature_data_frame, vittelogenic_adults, "VITADULT");
23
+ − 1291 # Standard error for vittelogenic adults.
24
+ − 1292 vittelogenic_adults.std_error = apply(Vittelogenic.replications, 1, sd) / sqrt(opt$replications);
27
+ − 1293 temperature_data_frame = append_vector(temperature_data_frame, vittelogenic_adults.std_error, "VITADULTSE");
18
+ − 1294 } else if (life_stage_adult == "Diapausing") {
23
+ − 1295 # Mean value for vittelogenic adults.
16
+ − 1296 diapausing_adults = apply(Diapausing.replications, 1, mean);
27
+ − 1297 temperature_data_frame = append_vector(temperature_data_frame, diapausing_adults, "DIAPAUSINGADULT");
23
+ − 1298 # Standard error for vittelogenic adults.
16
+ − 1299 diapausing_adults.std_error = apply(Diapausing.replications, 1, sd) / sqrt(opt$replications);
27
+ − 1300 temperature_data_frame = append_vector(temperature_data_frame, diapausing_adults.std_error, "DIAPAUSINGADULTSE");
28
+ − 1301 } else if (life_stage_adult=="Total") {
+ − 1302 # Mean value for all adults.
+ − 1303 total_adults = apply((Previttelogenic.replications+Vittelogenic.replications+Diapausing.replications), 1, mean);
+ − 1304 temperature_data_frame = append_vector(temperature_data_frame, total_adults, "TOTALADULT");
+ − 1305 # Standard error for all adults.
+ − 1306 total_adults.std_error = apply((Previttelogenic.replications+Vittelogenic.replications+Diapausing.replications), 1, sd) / sqrt(opt$replications);
+ − 1307 temperature_data_frame = append_vector(temperature_data_frame, total_adults.std_error, "TOTALADULTSE");
16
+ − 1308 }
10
+ − 1309 }
+ − 1310 }
5
+ − 1311
10
+ − 1312 if (plot_generations_separately) {
20
+ − 1313 m_se = get_mean_and_std_error(P.replications, F1.replications, F2.replications);
+ − 1314 P = m_se[[1]];
+ − 1315 P.std_error = m_se[[2]];
+ − 1316 F1 = m_se[[3]];
+ − 1317 F1.std_error = m_se[[4]];
+ − 1318 F2 = m_se[[5]];
+ − 1319 F2.std_error = m_se[[6]];
10
+ − 1320 if (process_eggs) {
20
+ − 1321 m_se = get_mean_and_std_error(P_eggs.replications, F1_eggs.replications, F2_eggs.replications);
+ − 1322 P_eggs = m_se[[1]];
+ − 1323 P_eggs.std_error = m_se[[2]];
31
+ − 1324 temperature_data_frame_P = append_vector(temperature_data_frame_P, P_eggs, "EGG-P");
+ − 1325 temperature_data_frame_P = append_vector(temperature_data_frame_P, P_eggs.std_error, "EGG-P-SE");
20
+ − 1326 F1_eggs = m_se[[3]];
+ − 1327 F1_eggs.std_error = m_se[[4]];
31
+ − 1328 temperature_data_frame_F1 = append_vector(temperature_data_frame_F1, F1_eggs, "EGG-F1");
+ − 1329 temperature_data_frame_F1 = append_vector(temperature_data_frame_F1, F1_eggs.std_error, "EGG-F1-SE");
20
+ − 1330 F2_eggs = m_se[[5]];
+ − 1331 F2_eggs.std_error = m_se[[6]];
31
+ − 1332 temperature_data_frame_F2 = append_vector(temperature_data_frame_F2, F2_eggs, "EGG-F2");
+ − 1333 temperature_data_frame_F2 = append_vector(temperature_data_frame_F2, F2_eggs.std_error, "EGG-F2-SE");
20
+ − 1334 }
+ − 1335 if (process_young_nymphs) {
+ − 1336 m_se = get_mean_and_std_error(P_young_nymphs.replications, F1_young_nymphs.replications, F2_young_nymphs.replications);
+ − 1337 P_young_nymphs = m_se[[1]];
+ − 1338 P_young_nymphs.std_error = m_se[[2]];
31
+ − 1339 temperature_data_frame_P = append_vector(temperature_data_frame_P, P_young_nymphs, "YOUNGNYMPH-P");
+ − 1340 temperature_data_frame_P = append_vector(temperature_data_frame_P, P_young_nymphs.std_error, "YOUNGNYMPH-P-SE");
20
+ − 1341 F1_young_nymphs = m_se[[3]];
+ − 1342 F1_young_nymphs.std_error = m_se[[4]];
31
+ − 1343 temperature_data_frame_F1 = append_vector(temperature_data_frame_F1, F1_young_nymphs, "YOUNGNYMPH-F1");
+ − 1344 temperature_data_frame_F1 = append_vector(temperature_data_frame_F1, F1_young_nymphs.std_error, "YOUNGNYMPH-F1-SE");
20
+ − 1345 F2_young_nymphs = m_se[[5]];
+ − 1346 F2_young_nymphs.std_error = m_se[[6]];
31
+ − 1347 temperature_data_frame_F2 = append_vector(temperature_data_frame_F2, F2_young_nymphs, "YOUNGNYMPH-F2");
+ − 1348 temperature_data_frame_F2 = append_vector(temperature_data_frame_F2, F2_young_nymphs.std_error, "YOUNGNYMPH-F2-SE");
10
+ − 1349 }
20
+ − 1350 if (process_old_nymphs) {
+ − 1351 m_se = get_mean_and_std_error(P_old_nymphs.replications, F1_old_nymphs.replications, F2_old_nymphs.replications);
+ − 1352 P_old_nymphs = m_se[[1]];
+ − 1353 P_old_nymphs.std_error = m_se[[2]];
31
+ − 1354 temperature_data_frame_P = append_vector(temperature_data_frame_P, P_old_nymphs, "OLDNYMPH-P");
+ − 1355 temperature_data_frame_P = append_vector(temperature_data_frame_P, P_old_nymphs.std_error, "OLDNYMPH-P-SE");
20
+ − 1356 F1_old_nymphs = m_se[[3]];
+ − 1357 F1_old_nymphs.std_error = m_se[[4]];
31
+ − 1358 temperature_data_frame_F1 = append_vector(temperature_data_frame_F1, F1_old_nymphs, "OLDNYMPH-F1");
+ − 1359 temperature_data_frame_F1 = append_vector(temperature_data_frame_F1, F1_old_nymphs.std_error, "OLDNYMPH-F1-SE");
20
+ − 1360 F2_old_nymphs = m_se[[5]];
+ − 1361 F2_old_nymphs.std_error = m_se[[6]];
31
+ − 1362 temperature_data_frame_F2 = append_vector(temperature_data_frame_F2, F2_old_nymphs, "OLDNYMPH-F2");
+ − 1363 temperature_data_frame_F2 = append_vector(temperature_data_frame_F2, F2_old_nymphs.std_error, "OLDNYMPH-F2-SE");
20
+ − 1364 }
+ − 1365 if (process_total_nymphs) {
+ − 1366 m_se = get_mean_and_std_error(P_total_nymphs.replications, F1_total_nymphs.replications, F2_total_nymphs.replications);
+ − 1367 P_total_nymphs = m_se[[1]];
+ − 1368 P_total_nymphs.std_error = m_se[[2]];
31
+ − 1369 temperature_data_frame_P = append_vector(temperature_data_frame_P, P_total_nymphs, "TOTALNYMPH-P");
+ − 1370 temperature_data_frame_P = append_vector(temperature_data_frame_P, P_total_nymphs.std_error, "TOTALNYMPH-P-SE");
20
+ − 1371 F1_total_nymphs = m_se[[3]];
+ − 1372 F1_total_nymphs.std_error = m_se[[4]];
31
+ − 1373 temperature_data_frame_F1 = append_vector(temperature_data_frame_F1, F1_total_nymphs, "TOTALNYMPH-F1");
+ − 1374 temperature_data_frame_F1 = append_vector(temperature_data_frame_F1, F1_total_nymphs.std_error, "TOTALNYMPH-F1-SE");
20
+ − 1375 F2_total_nymphs = m_se[[5]];
+ − 1376 F2_total_nymphs.std_error = m_se[[6]];
31
+ − 1377 temperature_data_frame_F2 = append_vector(temperature_data_frame_F2, F2_total_nymphs, "TOTALNYMPH-F2");
+ − 1378 temperature_data_frame_F2 = append_vector(temperature_data_frame_F2, F2_total_nymphs.std_error, "TOTALNYMPH-F2-SE");
10
+ − 1379 }
23
+ − 1380 if (process_previttelogenic_adults) {
+ − 1381 m_se = get_mean_and_std_error(P_previttelogenic_adults.replications, F1_previttelogenic_adults.replications, F2_previttelogenic_adults.replications);
+ − 1382 P_previttelogenic_adults = m_se[[1]];
+ − 1383 P_previttelogenic_adults.std_error = m_se[[2]];
31
+ − 1384 temperature_data_frame_P = append_vector(temperature_data_frame_P, P_previttelogenic_adults, "PRE-VITADULT-P");
+ − 1385 temperature_data_frame_P = append_vector(temperature_data_frame_P, P_previttelogenic_adults.std_error, "PRE-VITADULT-P-SE");
23
+ − 1386 F1_previttelogenic_adults = m_se[[3]];
+ − 1387 F1_previttelogenic_adults.std_error = m_se[[4]];
31
+ − 1388 temperature_data_frame_F1 = append_vector(temperature_data_frame_F1, F1_previttelogenic_adults, "PRE-VITADULT-F1");
+ − 1389 temperature_data_frame_F1 = append_vector(temperature_data_frame_F1, F1_previttelogenic_adults.std_error, "PRE-VITADULT-F1-SE");
23
+ − 1390 F2_previttelogenic_adults = m_se[[5]];
+ − 1391 F2_previttelogenic_adults.std_error = m_se[[6]];
31
+ − 1392 temperature_data_frame_F2 = append_vector(temperature_data_frame_F2, F2_previttelogenic_adults, "PRE-VITADULT-F2");
+ − 1393 temperature_data_frame_F2 = append_vector(temperature_data_frame_F2, F2_previttelogenic_adults.std_error, "PRE-VITADULT-F2-SE");
23
+ − 1394 }
+ − 1395 if (process_vittelogenic_adults) {
+ − 1396 m_se = get_mean_and_std_error(P_vittelogenic_adults.replications, F1_vittelogenic_adults.replications, F2_vittelogenic_adults.replications);
+ − 1397 P_vittelogenic_adults = m_se[[1]];
+ − 1398 P_vittelogenic_adults.std_error = m_se[[2]];
31
+ − 1399 temperature_data_frame_P = append_vector(temperature_data_frame_P, P_vittelogenic_adults, "VITADULT-P");
+ − 1400 temperature_data_frame_P = append_vector(temperature_data_frame_P, P_vittelogenic_adults.std_error, "VITADULT-P-SE");
23
+ − 1401 F1_vittelogenic_adults = m_se[[3]];
+ − 1402 F1_vittelogenic_adults.std_error = m_se[[4]];
31
+ − 1403 temperature_data_frame_F1 = append_vector(temperature_data_frame_F1, F1_vittelogenic_adults, "VITADULT-F1");
+ − 1404 temperature_data_frame_F1 = append_vector(temperature_data_frame_F1, F1_vittelogenic_adults.std_error, "VITADULT-F1-SE");
23
+ − 1405 F2_vittelogenic_adults = m_se[[5]];
+ − 1406 F2_vittelogenic_adults.std_error = m_se[[6]];
31
+ − 1407 temperature_data_frame_F2 = append_vector(temperature_data_frame_F2, F2_vittelogenic_adults, "VITADULT-F2");
+ − 1408 temperature_data_frame_F2 = append_vector(temperature_data_frame_F2, F2_vittelogenic_adults.std_error, "VITADULT-F2-SE");
23
+ − 1409 }
+ − 1410 if (process_diapausing_adults) {
+ − 1411 m_se = get_mean_and_std_error(P_diapausing_adults.replications, F1_diapausing_adults.replications, F2_diapausing_adults.replications);
+ − 1412 P_diapausing_adults = m_se[[1]];
+ − 1413 P_diapausing_adults.std_error = m_se[[2]];
31
+ − 1414 temperature_data_frame_P = append_vector(temperature_data_frame_P, P_diapausing_adults, "DIAPAUSINGADULT-P");
+ − 1415 temperature_data_frame_P = append_vector(temperature_data_frame_P, P_diapausing_adults.std_error, "DIAPAUSINGADULT-P-SE");
23
+ − 1416 F1_diapausing_adults = m_se[[3]];
+ − 1417 F1_diapausing_adults.std_error = m_se[[4]];
31
+ − 1418 temperature_data_frame_F1 = append_vector(temperature_data_frame_F1, F1_diapausing_adults, "DIAPAUSINGADULT-F1");
+ − 1419 temperature_data_frame_F1 = append_vector(temperature_data_frame_F1, F1_diapausing_adults.std_error, "DIAPAUSINGADULT-F1-SE");
23
+ − 1420 F2_diapausing_adults = m_se[[5]];
+ − 1421 F2_diapausing_adults.std_error = m_se[[6]];
31
+ − 1422 temperature_data_frame_F2 = append_vector(temperature_data_frame_F2, F2_diapausing_adults, "DIAPAUSINGADULT-F2");
+ − 1423 temperature_data_frame_F2 = append_vector(temperature_data_frame_F2, F2_diapausing_adults.std_error, "DIAPAUSINGADULT-F2-SE");
23
+ − 1424 }
+ − 1425 if (process_total_adults) {
+ − 1426 m_se = get_mean_and_std_error(P_total_adults.replications, F1_total_adults.replications, F2_total_adults.replications);
+ − 1427 P_total_adults = m_se[[1]];
+ − 1428 P_total_adults.std_error = m_se[[2]];
31
+ − 1429 temperature_data_frame_P = append_vector(temperature_data_frame_P, P_total_adults, "TOTALADULT-P");
+ − 1430 temperature_data_frame_P = append_vector(temperature_data_frame_P, P_total_adults.std_error, "TOTALADULT-P-SE");
23
+ − 1431 F1_total_adults = m_se[[3]];
+ − 1432 F1_total_adults.std_error = m_se[[4]];
31
+ − 1433 temperature_data_frame_F1 = append_vector(temperature_data_frame_F1, F1_total_adults, "TOTALADULT-F1");
+ − 1434 temperature_data_frame_F1 = append_vector(temperature_data_frame_F1, F1_total_adults.std_error, "TOTALADULT-F1-SE");
23
+ − 1435 F2_total_adults = m_se[[5]];
+ − 1436 F2_total_adults.std_error = m_se[[6]];
31
+ − 1437 temperature_data_frame_F2 = append_vector(temperature_data_frame_F2, F2_total_adults, "TOTALADULT-F2");
+ − 1438 temperature_data_frame_F2 = append_vector(temperature_data_frame_F2, F2_total_adults.std_error, "TOTALADULT-F2-SE");
10
+ − 1439 }
+ − 1440 }
6
+ − 1441
31
+ − 1442 # Save the analyzed data for combined generations.
34
+ − 1443 file_path = paste("output_data_dir", "04_combined_generations.csv", sep="/");
+ − 1444 write.csv(temperature_data_frame, file=file_path, row.names=F);
31
+ − 1445 if (plot_generations_separately) {
+ − 1446 # Save the analyzed data for generation P.
34
+ − 1447 file_path = paste("output_data_dir", "01_generation_P.csv", sep="/");
+ − 1448 write.csv(temperature_data_frame_P, file=file_path, row.names=F);
31
+ − 1449 # Save the analyzed data for generation F1.
34
+ − 1450 file_path = paste("output_data_dir", "02_generation_F1.csv", sep="/");
+ − 1451 write.csv(temperature_data_frame_F1, file=file_path, row.names=F);
31
+ − 1452 # Save the analyzed data for generation F2.
34
+ − 1453 file_path = paste("output_data_dir", "03_generation_F2.csv", sep="/");
+ − 1454 write.csv(temperature_data_frame_F2, file=file_path, row.names=F);
31
+ − 1455 }
5
+ − 1456
10
+ − 1457 if (plot_generations_separately) {
15
+ − 1458 for (life_stage in life_stages) {
10
+ − 1459 if (life_stage == "Egg") {
+ − 1460 # Start PDF device driver.
+ − 1461 dev.new(width=20, height=30);
19
+ − 1462 file_path = get_file_path(life_stage, "egg_pop_by_generation.pdf")
10
+ − 1463 pdf(file=file_path, width=20, height=30, bg="white");
+ − 1464 par(mar=c(5, 6, 4, 4), mfrow=c(3, 1));
+ − 1465 # Egg population size by generation.
18
+ − 1466 maxval = max(P_eggs+F1_eggs+F2_eggs) + 100;
38
+ − 1467 render_chart(ticks, date_labels, "pop_size_by_generation", opt$plot_std_error, opt$insect, opt$location, latitude,
+ − 1468 start_date, end_date, total_days_vector, maxval, opt$replications, life_stage, group=P_eggs, group_std_error=P_eggs.std_error,
+ − 1469 group2=F1_eggs, group2_std_error=F1_eggs.std_error, group3=F2_eggs, group3_std_error=F2_eggs.std_error);
10
+ − 1470 # Turn off device driver to flush output.
+ − 1471 dev.off();
+ − 1472 } else if (life_stage == "Nymph") {
16
+ − 1473 for (life_stage_nymph in life_stages_nymph) {
+ − 1474 # Start PDF device driver.
+ − 1475 dev.new(width=20, height=30);
19
+ − 1476 file_path = get_file_path(life_stage, "nymph_pop_by_generation.pdf", life_stage_nymph=life_stage_nymph)
16
+ − 1477 pdf(file=file_path, width=20, height=30, bg="white");
+ − 1478 par(mar=c(5, 6, 4, 4), mfrow=c(3, 1));
20
+ − 1479 if (life_stage_nymph=="Young") {
+ − 1480 # Young nymph population size by generation.
+ − 1481 maxval = max(P_young_nymphs+F1_young_nymphs+F2_young_nymphs) + 100;
+ − 1482 group = P_young_nymphs;
+ − 1483 group_std_error = P_young_nymphs.std_error;
+ − 1484 group2 = F1_young_nymphs;
+ − 1485 group2_std_error = F1_young_nymphs.std_error;
+ − 1486 group3 = F2_young_nymphs;
+ − 1487 group3_std_error = F2_young_nymphs.std_error;
+ − 1488 } else if (life_stage_nymph=="Old") {
+ − 1489 # Total nymph population size by generation.
+ − 1490 maxval = max(P_old_nymphs+F1_old_nymphs+F2_old_nymphs) + 100;
+ − 1491 group = P_old_nymphs;
+ − 1492 group_std_error = P_old_nymphs.std_error;
+ − 1493 group2 = F1_old_nymphs;
+ − 1494 group2_std_error = F1_old_nymphs.std_error;
+ − 1495 group3 = F2_old_nymphs;
+ − 1496 group3_std_error = F2_old_nymphs.std_error;
+ − 1497 } else if (life_stage_nymph=="Total") {
+ − 1498 # Total nymph population size by generation.
+ − 1499 maxval = max(P_total_nymphs+F1_total_nymphs+F2_total_nymphs) + 100;
+ − 1500 group = P_total_nymphs;
+ − 1501 group_std_error = P_total_nymphs.std_error;
+ − 1502 group2 = F1_total_nymphs;
+ − 1503 group2_std_error = F1_total_nymphs.std_error;
+ − 1504 group3 = F2_total_nymphs;
+ − 1505 group3_std_error = F2_total_nymphs.std_error;
+ − 1506 }
38
+ − 1507 render_chart(ticks, date_labels, "pop_size_by_generation", opt$plot_std_error, opt$insect, opt$location, latitude,
+ − 1508 start_date, end_date, total_days_vector, maxval, opt$replications, life_stage, group=group, group_std_error=group_std_error,
+ − 1509 group2=group2, group2_std_error=group2_std_error, group3=group3, group3_std_error=group3_std_error, life_stages_nymph=life_stage_nymph);
16
+ − 1510 # Turn off device driver to flush output.
+ − 1511 dev.off();
+ − 1512 }
10
+ − 1513 } else if (life_stage == "Adult") {
16
+ − 1514 for (life_stage_adult in life_stages_adult) {
+ − 1515 # Start PDF device driver.
+ − 1516 dev.new(width=20, height=30);
19
+ − 1517 file_path = get_file_path(life_stage, "adult_pop_by_generation.pdf", life_stage_adult=life_stage_adult)
16
+ − 1518 pdf(file=file_path, width=20, height=30, bg="white");
+ − 1519 par(mar=c(5, 6, 4, 4), mfrow=c(3, 1));
23
+ − 1520 if (life_stage_adult=="Pre-vittelogenic") {
+ − 1521 # Pre-vittelogenic adult population size by generation.
+ − 1522 maxval = max(P_previttelogenic_adults+F1_previttelogenic_adults+F2_previttelogenic_adults) + 100;
+ − 1523 group = P_previttelogenic_adults;
+ − 1524 group_std_error = P_previttelogenic_adults.std_error;
+ − 1525 group2 = F1_previttelogenic_adults;
+ − 1526 group2_std_error = F1_previttelogenic_adults.std_error;
+ − 1527 group3 = F2_previttelogenic_adults;
+ − 1528 group3_std_error = F2_previttelogenic_adults.std_error;
+ − 1529 } else if (life_stage_adult=="Vittelogenic") {
+ − 1530 # Vittelogenic adult population size by generation.
+ − 1531 maxval = max(P_vittelogenic_adults+F1_vittelogenic_adults+F2_vittelogenic_adults) + 100;
+ − 1532 group = P_vittelogenic_adults;
+ − 1533 group_std_error = P_vittelogenic_adults.std_error;
+ − 1534 group2 = F1_vittelogenic_adults;
+ − 1535 group2_std_error = F1_vittelogenic_adults.std_error;
+ − 1536 group3 = F2_vittelogenic_adults;
+ − 1537 group3_std_error = F2_vittelogenic_adults.std_error;
+ − 1538 } else if (life_stage_adult=="Diapausing") {
+ − 1539 # Diapausing adult population size by generation.
+ − 1540 maxval = max(P_diapausing_adults+F1_diapausing_adults+F2_diapausing_adults) + 100;
+ − 1541 group = P_diapausing_adults;
+ − 1542 group_std_error = P_diapausing_adults.std_error;
+ − 1543 group2 = F1_diapausing_adults;
+ − 1544 group2_std_error = F1_diapausing_adults.std_error;
+ − 1545 group3 = F2_diapausing_adults;
+ − 1546 group3_std_error = F2_diapausing_adults.std_error;
+ − 1547 } else if (life_stage_adult=="Total") {
+ − 1548 # Total adult population size by generation.
+ − 1549 maxval = max(P_total_adults+F1_total_adults+F2_total_adults) + 100;
+ − 1550 group = P_total_adults;
+ − 1551 group_std_error = P_total_adults.std_error;
+ − 1552 group2 = F1_total_adults;
+ − 1553 group2_std_error = F1_total_adults.std_error;
+ − 1554 group3 = F2_total_adults;
+ − 1555 group3_std_error = F2_total_adults.std_error;
+ − 1556 }
38
+ − 1557 render_chart(ticks, date_labels, "pop_size_by_generation", opt$plot_std_error, opt$insect, opt$location, latitude,
+ − 1558 start_date, end_date, total_days_vector, maxval, opt$replications, life_stage, group=group, group_std_error=group_std_error,
+ − 1559 group2=group2, group2_std_error=group2_std_error, group3=group3, group3_std_error=group3_std_error, life_stages_adult=life_stage_adult);
16
+ − 1560 # Turn off device driver to flush output.
+ − 1561 dev.off();
+ − 1562 }
10
+ − 1563 } else if (life_stage == "Total") {
+ − 1564 # Start PDF device driver.
18
+ − 1565 # Name collection elements so that they
+ − 1566 # are displayed in logical order.
10
+ − 1567 dev.new(width=20, height=30);
19
+ − 1568 file_path = get_file_path(life_stage, "total_pop_by_generation.pdf")
10
+ − 1569 pdf(file=file_path, width=20, height=30, bg="white");
+ − 1570 par(mar=c(5, 6, 4, 4), mfrow=c(3, 1));
+ − 1571 # Total population size by generation.
18
+ − 1572 maxval = max(P+F1+F2) + 100;
38
+ − 1573 render_chart(ticks, date_labels, "pop_size_by_generation", opt$plot_std_error, opt$insect, opt$location, latitude,
+ − 1574 start_date, end_date, total_days_vector, maxval, opt$replications, life_stage, group=P, group_std_error=P.std_error,
+ − 1575 group2=F1, group2_std_error=F1.std_error, group3=F2, group3_std_error=F2.std_error);
10
+ − 1576 # Turn off device driver to flush output.
+ − 1577 dev.off();
+ − 1578 }
15
+ − 1579 }
10
+ − 1580 } else {
+ − 1581 for (life_stage in life_stages) {
+ − 1582 if (life_stage == "Egg") {
+ − 1583 # Start PDF device driver.
+ − 1584 dev.new(width=20, height=30);
19
+ − 1585 file_path = get_file_path(life_stage, "egg_pop.pdf")
10
+ − 1586 pdf(file=file_path, width=20, height=30, bg="white");
+ − 1587 par(mar=c(5, 6, 4, 4), mfrow=c(3, 1));
+ − 1588 # Egg population size.
18
+ − 1589 maxval = max(eggs+eggs.std_error) + 100;
38
+ − 1590 render_chart(ticks, date_labels, "pop_size_by_life_stage", opt$plot_std_error, opt$insect, opt$location, latitude,
+ − 1591 start_date, end_date, total_days_vector, maxval, opt$replications, life_stage, group=eggs, group_std_error=eggs.std_error);
10
+ − 1592 # Turn off device driver to flush output.
+ − 1593 dev.off();
+ − 1594 } else if (life_stage == "Nymph") {
16
+ − 1595 for (life_stage_nymph in life_stages_nymph) {
+ − 1596 # Start PDF device driver.
+ − 1597 dev.new(width=20, height=30);
19
+ − 1598 file_path = get_file_path(life_stage, "nymph_pop.pdf", life_stage_nymph=life_stage_nymph)
16
+ − 1599 pdf(file=file_path, width=20, height=30, bg="white");
+ − 1600 par(mar=c(5, 6, 4, 4), mfrow=c(3, 1));
+ − 1601 if (life_stage_nymph=="Total") {
+ − 1602 # Total nymph population size.
+ − 1603 group = total_nymphs;
+ − 1604 group_std_error = total_nymphs.std_error;
+ − 1605 } else if (life_stage_nymph=="Young") {
+ − 1606 # Young nymph population size.
+ − 1607 group = young_nymphs;
+ − 1608 group_std_error = young_nymphs.std_error;
+ − 1609 } else if (life_stage_nymph=="Old") {
+ − 1610 # Old nymph population size.
+ − 1611 group = old_nymphs;
+ − 1612 group_std_error = old_nymphs.std_error;
+ − 1613 }
18
+ − 1614 maxval = max(group+group_std_error) + 100;
38
+ − 1615 render_chart(ticks, date_labels, "pop_size_by_life_stage", opt$plot_std_error, opt$insect, opt$location, latitude,
+ − 1616 start_date, end_date, total_days_vector, maxval, opt$replications, life_stage, group=group, group_std_error=group_std_error,
+ − 1617 life_stages_nymph=life_stage_nymph);
16
+ − 1618 # Turn off device driver to flush output.
+ − 1619 dev.off();
+ − 1620 }
10
+ − 1621 } else if (life_stage == "Adult") {
16
+ − 1622 for (life_stage_adult in life_stages_adult) {
+ − 1623 # Start PDF device driver.
+ − 1624 dev.new(width=20, height=30);
19
+ − 1625 file_path = get_file_path(life_stage, "adult_pop.pdf", life_stage_adult=life_stage_adult)
16
+ − 1626 pdf(file=file_path, width=20, height=30, bg="white");
+ − 1627 par(mar=c(5, 6, 4, 4), mfrow=c(3, 1));
+ − 1628 if (life_stage_adult=="Total") {
+ − 1629 # Total adult population size.
+ − 1630 group = total_adults;
+ − 1631 group_std_error = total_adults.std_error
+ − 1632 } else if (life_stage_adult=="Pre-vittelogenic") {
+ − 1633 # Pre-vittelogenic adult population size.
+ − 1634 group = previttelogenic_adults;
+ − 1635 group_std_error = previttelogenic_adults.std_error
+ − 1636 } else if (life_stage_adult=="Vittelogenic") {
+ − 1637 # Vittelogenic adult population size.
+ − 1638 group = vittelogenic_adults;
+ − 1639 group_std_error = vittelogenic_adults.std_error
+ − 1640 } else if (life_stage_adult=="Diapausing") {
+ − 1641 # Diapausing adult population size.
+ − 1642 group = diapausing_adults;
+ − 1643 group_std_error = diapausing_adults.std_error
+ − 1644 }
18
+ − 1645 maxval = max(group+group_std_error) + 100;
38
+ − 1646 render_chart(ticks, date_labels, "pop_size_by_life_stage", opt$plot_std_error, opt$insect, opt$location, latitude,
+ − 1647 start_date, end_date, total_days_vector, maxval, opt$replications, life_stage, group=group, group_std_error=group_std_error,
+ − 1648 life_stages_adult=life_stage_adult);
16
+ − 1649 # Turn off device driver to flush output.
+ − 1650 dev.off();
+ − 1651 }
10
+ − 1652 } else if (life_stage == "Total") {
+ − 1653 # Start PDF device driver.
+ − 1654 dev.new(width=20, height=30);
19
+ − 1655 file_path = get_file_path(life_stage, "total_pop.pdf")
10
+ − 1656 pdf(file=file_path, width=20, height=30, bg="white");
+ − 1657 par(mar=c(5, 6, 4, 4), mfrow=c(3, 1));
+ − 1658 # Total population size.
18
+ − 1659 maxval = max(eggs+eggs.std_error, total_nymphs+total_nymphs.std_error, total_adults+total_adults.std_error) + 100;
38
+ − 1660 render_chart(ticks, date_labels, "pop_size_by_life_stage", opt$plot_std_error, opt$insect, opt$location, latitude,
+ − 1661 start_date, end_date, total_days_vector, maxval, opt$replications, life_stage, group=total_adults, group_std_error=total_adults.std_error,
+ − 1662 group2=total_nymphs, group2_std_error=total_nymphs.std_error, group3=eggs, group3_std_error=eggs.std_error);
10
+ − 1663 # Turn off device driver to flush output.
+ − 1664 dev.off();
+ − 1665 }
+ − 1666 }
+ − 1667 }
