Mercurial > repos > greg > bmsb

comparison bmsb.R @ 37:d3a6b291e096 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author greg
date Sun, 18 Dec 2016 13:28:04 -0500
parents 7c40c2b303f1
children
comparison
equal deleted inserted replaced
36:956daea2c7fe 37:d3a6b291e096
24 24
25 parser <- OptionParser(usage="%prog [options] file", option_list=option_list) 25 parser <- OptionParser(usage="%prog [options] file", option_list=option_list)
26 args <- parse_args(parser, positional_arguments=TRUE) 26 args <- parse_args(parser, positional_arguments=TRUE)
27 opt <- args$options 27 opt <- args$options
28 28
29 data.input=function(loc, start.yr, temperature.dataset) 29 data.input=function(loc, year, temperature.dataset)
30 { 30 {
31 expdata <- matrix(rep(0, 365 * 3), nrow=365) 31 expdata <- matrix(rep(0, 365 * 3), nrow=365)
32 # replace 2004 with start. yr 32 namedat <- paste(loc, year, ".Rdat", sep="")
33 yr <- start.yr
34 namedat <- paste(loc, yr, ".Rdat", sep="")
35 temp.data <- read.csv(file=temperature.dataset, header=T) 33 temp.data <- read.csv(file=temperature.dataset, header=T)
36 34
37 expdata[,1] <- c(1:365) 35 expdata[,1] <- c(1:365)
36 # Minimum
37 expdata[,2] <- temp.data[c(1:365), 3]
38 # Maximum
39 expdata[,3] <- temp.data[c(1:365), 2]
38 save(expdata, file=namedat) 40 save(expdata, file=namedat)
39 namedat 41 namedat
40 } 42 }
41 43
42 daylength=function(L) 44 daylength=function(latitude)
43 { 45 {
44 # from Forsythe 1995 46 # from Forsythe 1995
45 p=0.8333 47 p=0.8333
46 dl <- NULL 48 dl <- NULL
47 for (i in 1:365) 49 for (i in 1:365)
48 { 50 {
49 theta <- 0.2163108 + 2 * atan(0.9671396 * tan(0.00860 * (i - 186))) 51 theta <- 0.2163108 + 2 * atan(0.9671396 * tan(0.00860 * (i - 186)))
50 phi <- asin(0.39795 * cos(theta)) 52 phi <- asin(0.39795 * cos(theta))
51 dl[i] <- 24 - 24 / pi * acos((sin(p * pi / 180) + sin(L * pi / 180) * sin(phi)) / (cos(L * pi / 180) * cos(phi))) 53 dl[i] <- 24 - 24 / pi * acos((sin(p * pi / 180) + sin(latitude * pi / 180) * sin(phi)) / (cos(latitude * pi / 180) * cos(phi)))
52 } 54 }
53 dl # return a vector of daylength in 365 days 55 dl # return a vector of daylength in 365 days
54 } 56 }
55 57
56 hourtemp=function(L, date, temperature_file_path) 58 hourtemp=function(latitude, date, temperature_file_path)
57 { 59 {
58 load(temperature_file_path) 60 load(temperature_file_path)
59 threshold <- 14.17 # base development threshold for BMSB 61 threshold <- 14.17 # base development threshold for BMSB
60 dnp <- expdata[date, 2] # daily minimum 62 dnp <- expdata[date, 2] # daily minimum
61 dxp <- expdata[date, 3] # daily maximum 63 dxp <- expdata[date, 3] # daily maximum
66 { 68 {
67 dd <- 0 69 dd <- 0
68 } 70 }
69 else 71 else
70 { 72 {
71 dlprofile <- daylength(L) # extract daylength data for entire year 73 dlprofile <- daylength(latitude) # extract daylength data for entire year
72 T <- NULL # initialize hourly temperature 74 T <- NULL # initialize hourly temperature
73 dh <- NULL #initialize degree hour vector 75 dh <- NULL #initialize degree hour vector
74 # date <- 200 76 # date <- 200
75 y <- dlprofile[date] # calculate daylength in given date 77 y <- dlprofile[date] # calculate daylength in given date
76 z <- 24 - y # night length 78 z <- 24 - y # night length
363 { 365 {
364 # vittelogenic stage, overwintering generation 366 # vittelogenic stage, overwintering generation
365 if (vec.ind[4] == 0) 367 if (vec.ind[4] == 0)
366 { 368 {
367 # just turned in vittelogenic stage 369 # just turned in vittelogenic stage
368 n.birth=round(runif(1, 2 + min.ovi.adj, 8 + max.ovi.adj)) 370 n.birth=round(runif(1, 2 + opt$min_clutch_size, 8 + opt$max_clutch_size))
369 } 371 }
370 else 372 else
371 { 373 {
372 # daily probability of birth 374 # daily probability of birth
373 p.birth = opt$oviposition * 0.01 375 p.birth = opt$oviposition * 0.01
402 { 404 {
403 # vittelogenic stage, 1st generation 405 # vittelogenic stage, 1st generation
404 if (vec.ind[4] == 0) 406 if (vec.ind[4] == 0)
405 { 407 {
406 # just turned in vittelogenic stage 408 # just turned in vittelogenic stage
407 n.birth=round(runif(1, 2 + min.ovi.adj, 8 + max.ovi.adj)) 409 n.birth=round(runif(1, 2 + opt$min_clutch_size, 8 + opt$max_clutch_size))
408 } 410 }
409 else 411 else
410 { 412 {
411 # daily probability of birth 413 # daily probability of birth
412 p.birth = opt$oviposition * 0.01 414 p.birth = opt$oviposition * 0.01
441 if (vec.ind[2] == 0) 443 if (vec.ind[2] == 0)
442 { 444 {
443 # egg stage 445 # egg stage
444 # add to DD 446 # add to DD
445 vec.ind[3] <- vec.ind[3] + dd.temp 447 vec.ind[3] <- vec.ind[3] + dd.temp
446 if (vec.ind[3] >= (68 + dd.adj1)) 448 if (vec.ind[3] >= (68 + opt$young_nymph_accum))
447 { 449 {
448 # from egg to young nymph, DD requirement met 450 # from egg to young nymph, DD requirement met
449 current.gen <- vec.ind[1] 451 current.gen <- vec.ind[1]
450 # transfer to young nym stage 452 # transfer to young nym stage
451 vec.ind <- c(current.gen, 1, 0, 0, 0) 453 vec.ind <- c(current.gen, 1, 0, 0, 0)
462 if (vec.ind[2] == 1) 464 if (vec.ind[2] == 1)
463 { 465 {
464 # young nymph stage 466 # young nymph stage
465 # add to DD 467 # add to DD
466 vec.ind[3] <- vec.ind[3] + dd.temp 468 vec.ind[3] <- vec.ind[3] + dd.temp
467 if (vec.ind[3] >= (250 +dd.adj2)) 469 if (vec.ind[3] >= (250 + opt$old_nymph_accum))
468 { 470 {
469 # from young to old nymph, DD requirement met 471 # from young to old nymph, DD requirement met
470 current.gen <- vec.ind[1] 472 current.gen <- vec.ind[1]
471 # transfer to old nym stage 473 # transfer to old nym stage
472 vec.ind <- c(current.gen, 2, 0, 0, 0) 474 vec.ind <- c(current.gen, 2, 0, 0, 0)
487 if (vec.ind[2] == 2) 489 if (vec.ind[2] == 2)
488 { 490 {
489 # old nymph stage 491 # old nymph stage
490 # add to DD 492 # add to DD
491 vec.ind[3] <- vec.ind[3] + dd.temp 493 vec.ind[3] <- vec.ind[3] + dd.temp
492 if (vec.ind[3] >= (200 + dd.adj3)) 494 if (vec.ind[3] >= (200 + opt$adult_accum))
493 { 495 {
494 # from old to adult, DD requirement met 496 # from old to adult, DD requirement met
495 current.gen <- vec.ind[1] 497 current.gen <- vec.ind[1]
496 if (vec.ind[5] == 0) 498 if (vec.ind[5] == 0)
497 { 499 {