#!/usr/bin/env Rscript

suppressPackageStartupMessages(library("optparse"))

options_list <- list(
    make_option(c("-i", "--input"), action="store", help="Input dataset")
    make_option(c("-o", "--output"), action="store", help="Output dataset")
)

parser <- OptionParser(usage="%prog [options] file", options_list)
args <- parse_args(parser, positional_arguments=TRUE)
opt <- args$options

#########################################
daylength=function(L){
# from Forsythe 1995
p=0.8333
dl<-NULL
for (i in 1:365) {
theta<-0.2163108+2*atan(0.9671396*tan(0.00860*(i-186)))
phi<-asin(0.39795*cos(theta))
dl[i]<-24-24/pi*acos((sin(p*pi/180)+sin(L*pi/180)*sin(phi))/(cos(L*pi/180)*cos(phi)))
}
dl   # return a vector of daylength in 365 days
}
#########################################

#########################################
# source("daylength.R")
hourtemp=function(L,date){
# L=37.5 specify this in main program
threshold<-12.7  # base development threshold for BMSB
# threshold2<-threshold/24 degree hour accumulation
#expdata<-tempdata[1:365,11:13] # Use daily max, min, mean
dnp<-expdata[date,2]  # daily minimum
dxp<-expdata[date,3]  # daily maximum
dmean<-0.5*(dnp+dxp)
#if (dmean>0) {
#dnp<-dnp-k1*dmean  
#dxp<-dxp+k2*dmean 
#} else {
#dnp<-dnp+k1*dmean  
#dxp<-dxp-k2*dmean
#}
dd<-0  # initialize degree day accumulation

if (dxp<threshold) {dd<-0} else
{
dlprofile<-daylength(L)  # extract daylength data for entire year
T<-NULL  # initialize hourly temperature
dh<-NULL #initialize degree hour vector
# date<-200
y<-dlprofile[date]  # calculate daylength in given date
z<-24-y     # night length
a<-1.86     # lag coefficient
b<-2.20     # night coefficient
#tempdata<-read.csv("tempdata.csv") #import raw data set
# Should be outside function otherwise its redundant
risetime<-12-y/2      # sunrise time
settime<-12+y/2       # sunset time
ts<-(dxp-dnp)*sin(pi*(settime-5)/(y+2*a))+dnp
for (i in 1:24){      
  if (i>risetime && i<settime) {
    m<-i-5  # number of hours after Tmin until sunset
    T[i]=(dxp-dnp)*sin(pi*m/(y+2*a))+dnp
    if (T[i]<8.4) {dh[i]<-0} else
    {dh[i]<-T[i]-8.4}
    } else
  if (i>settime){ 
    n<-i-settime
    T[i]=dnp+(ts-dnp)*exp(-b*n/z)
    if (T[i]<8.4) {dh[i]<-0} else
    {dh[i]<-T[i]-8.4}
    } else
           {
    n<-i+24-settime
    T[i]=dnp+(ts-dnp)*exp(-b*n/z)
    if (T[i]<8.4) {dh[i]<-0} else
    {dh[i]<-T[i]-8.4}
    }
}
dd<-sum(dh)/24
}
return=c(dmean,dd)
return
}
#########################################


#########################################
mortality.egg=function(temperature){
if (temperature<12.7) {
 mort.prob=0.8} else 
{mort.prob=0.8-temperature/40
if (mort.prob<0) {mort.prob=0.01}
}
return=mort.prob
return
}
#########################################


#########################################
mortality.nymph=function(temperature){
if (temperature<12.7) {
 mort.prob=0.03} else 
{mort.prob=temperature*0.0008+0.03}
return=mort.prob
return
}
#########################################

#########################################
mortality.adult=function(temperature){
if (temperature<12.7) {
 mort.prob=0.002} else 
{mort.prob=temperature*0.0005+0.02}
return=mort.prob
return
}
#########################################

# model initialization
# setwd(“/home/lunarmouse/Dropbox/Nelson's project/") 
# PLEASE CHANGE TO YOUR OWN DIRECTORY!!!
# PLEASE LOAD BSMB FUNCTIONS FIRST!!!

n<-1000 # start with 1000 individuals
# Generation, Stage, DD, T, Diapause
vec.ini<-c(0,3,0,0,0)
# overwintering, previttelogenic,DD=0, T=0, no-diapause
vec.mat<-rep(vec.ini,n)
vec.mat<-t(matrix(vec.mat,nrow=5)) # complete matrix for the population
L<-35.58 # latitude for Asheville NC
ph.p<-daylength(L) # complete photoperiod profile in a year, requires daylength function

#load("asheville2014.Rdat") # load temperature data@location/year
load(opt$input) # load temperature data@location/year
tot.pop<-NULL # time series of population size
gen0.pop<-rep(0,365) # gen.0 pop size
gen1.pop<-rep(0,365)
gen2.pop<-rep(0,365)
S0<-S1<-S2<-S3<-S4<-S5<-rep(0,365)
g0.adult<-g1.adult<-g2.adult<-rep(0,365)
N.newborn<-N.death<-N.adult<-rep(0,365)
dd.day<-rep(0,365)

ptm <- proc.time() # start tick

for (day in 1:365) { # all the day
photoperiod<-ph.p[day] # photoperiod in the day
temp.profile<-hourtemp(L,day)
mean.temp<-temp.profile[1]
dd.temp<-temp.profile[2]
dd.day[day]<-dd.temp
death.vec<-NULL # trash bin for death
birth.vec<-NULL # new born
#n<-length(vec.mat[,1]) # population size at previous day

  for (i in 1:n) { # all individual
  vec.ind<-vec.mat[i,] # find individual record
  
# first of all, still alive?  
if(vec.ind[2]==0){ # egg
death.prob=mortality.egg(mean.temp)
 } else if (vec.ind[2]==1 | vec.ind[2]==2) {
death.prob=mortality.nymph(mean.temp)
}  else if (vec.ind[2]==3 | vec.ind[2]==4 | vec.ind[2]==5) { # for adult
 if (day<120 && day>270) {death.prob=0.33*mortality.adult(mean.temp)
  } else { 
death.prob=mortality.adult(mean.temp)} # reduce adult mortality after fall equinox
}


#(or dependent on temperature and life stage?)
u.d<-runif(1)
if (u.d<death.prob) {
death.vec<-c(death.vec,i)} else # aggregrate index of dead bug
  {
# event 1 end of diapause
if (vec.ind[1]==0 && vec.ind[2]==3) { # overwintering adult (previttelogenic)
 if (photoperiod>13.5 && vec.ind[3]>77 && day<180) { # add 77C to become fully reproductively matured
    vec.ind<-c(0,4,0,0,0) # transfer to vittelogenic
	vec.mat[i,]<-vec.ind
	
 } else {
    vec.ind[3]<-vec.ind[3]+dd.temp # add to DD
    vec.ind[4]<-vec.ind[4]+1 # add 1 day in current stage
	vec.mat[i,]<-vec.ind
 }
}

if (vec.ind[1]!=0 && vec.ind[2]==3) { # NOT overwintering adult (previttelogenic)
    current.gen<-vec.ind[1]
 if (vec.ind[3]>77) { # add 77C to become fully reproductively matured
    vec.ind<-c(current.gen,4,0,0,0) # transfer to vittelogenic
	vec.mat[i,]<-vec.ind
 } else {
    vec.ind[3]<-vec.ind[3]+dd.temp # add to DD
    vec.ind[4]<-vec.ind[4]+1 # add 1 day in current stage
	vec.mat[i,]<-vec.ind
 }
}

 
 
# event 2 oviposition -- where population dynamics comes from
if (vec.ind[2]==4 && vec.ind[1]==0 && mean.temp>10) { # vittelogenic stage, overwintering generation
   if (vec.ind[4]==0) { # just turned in vittelogenic stage
   n.birth=round(runif(1,2,8))} else{
   p.birth=0.01 # daily probability of birth
   u1<-runif(1)
      if (u1<p.birth) {n.birth=round(runif(1,2,8))}
   }
    vec.ind[3]<-vec.ind[3]+dd.temp # add to DD
    vec.ind[4]<-vec.ind[4]+1 # add 1 day in current stage
	vec.mat[i,]<-vec.ind
     if (n.birth>0) { # add new birth -- might be in different generations
	  new.gen<-vec.ind[1]+1 # generation +1
	  new.ind<-c(new.gen,0,0,0,0) # egg profile
	  new.vec<-rep(new.ind,n.birth)
	  new.vec<-t(matrix(new.vec,nrow=5)) # update batch of egg profile
	  birth.vec<-rbind(birth.vec,new.vec) # group with total eggs laid in that day
	 }
}

# event 2 oviposition -- for gen 1.
if (vec.ind[2]==4 && vec.ind[1]==1 && mean.temp>12.5 && day<222) { # vittelogenic stage, 1st generation
   if (vec.ind[4]==0) { # just turned in vittelogenic stage
   n.birth=round(runif(1,2,8))} else{
   p.birth=0.01 # daily probability of birth
   u1<-runif(1)
      if (u1<p.birth) {n.birth=round(runif(1,2,8))}
   }
    vec.ind[3]<-vec.ind[3]+dd.temp # add to DD
    vec.ind[4]<-vec.ind[4]+1 # add 1 day in current stage
	vec.mat[i,]<-vec.ind
     if (n.birth>0) { # add new birth -- might be in different generations
	  new.gen<-vec.ind[1]+1 # generation +1
	  new.ind<-c(new.gen,0,0,0,0) # egg profile
	  new.vec<-rep(new.ind,n.birth)
	  new.vec<-t(matrix(new.vec,nrow=5)) # update batch of egg profile
	  birth.vec<-rbind(birth.vec,new.vec) # group with total eggs laid in that day
	 }
}



# event 3 development (with diapause determination)
  # event 3.1 egg development to young nymph (vec.ind[2]=0 -> egg)
if (vec.ind[2]==0) { # egg stage
   vec.ind[3]<-vec.ind[3]+dd.temp # add to DD
    if (vec.ind[3]>=68) { # from egg to young nymph, DD requirement met
	    current.gen<-vec.ind[1]
		vec.ind<-c(current.gen,1,0,0,0) # transfer to young nym stage
	} else {
	vec.ind[4]<-vec.ind[4]+1 # add 1 day in current stage
	}
	vec.mat[i,]<-vec.ind
}
  
  # event 3.2 young nymph to old nymph (vec.ind[2]=1 -> young nymph: determines diapause)
if (vec.ind[2]==1) { # young nymph stage
      vec.ind[3]<-vec.ind[3]+dd.temp # add to DD
   if (vec.ind[3]>=250) { # from young to old nymph, DD requirement met
	    current.gen<-vec.ind[1]
		vec.ind<-c(current.gen,2,0,0,0) # transfer to old nym stage
		  if (photoperiod<13.5 && day > 180) {vec.ind[5]<-1} # prepare for diapausing
	} else {
	vec.ind[4]<-vec.ind[4]+1 # add 1 day in current stage
	}
    vec.mat[i,]<-vec.ind
}  
  
  
  # event 3.3 old nymph to adult: previttelogenic or diapausing?
if (vec.ind[2]==2) { # old nymph stage
    vec.ind[3]<-vec.ind[3]+dd.temp # add to DD
       if (vec.ind[3]>=200) { # from old to adult, DD requirement met
       current.gen<-vec.ind[1]
         if (vec.ind[5]==0) { # non-diapausing adult -- previttelogenic
          vec.ind<-c(current.gen,3,0,0,0)
        } else { # diapausing 
          vec.ind<-c(current.gen,5,0,0,1)
        }		
   } else {
   vec.ind[4]<-vec.ind[4]+1 # add 1 day in current stage
   }
   vec.mat[i,]<-vec.ind
}

# event 4 growing of diapausing adult (unimportant, but still necessary)## 
if (vec.ind[2]==5) {
  vec.ind[3]<-vec.ind[3]+dd.temp
  vec.ind[4]<-vec.ind[4]+1
  vec.mat[i,]<-vec.ind
}

 } # else if it is still alive

} # end of the individual bug loop

# find how many died
n.death<-length(death.vec)
if (n.death>0) {
vec.mat<-vec.mat[-death.vec, ]}
# remove record of dead
# find how many new born  
n.newborn<-length(birth.vec[,1])
vec.mat<-rbind(vec.mat,birth.vec)
# update population size for the next day
n<-n-n.death+n.newborn 

# aggregate results by day
tot.pop<-c(tot.pop,n) 
s0<-sum(vec.mat[,2]==0) #egg
s1<-sum(vec.mat[,2]==1) # young nymph
s2<-sum(vec.mat[,2]==2) # old nymph
s3<-sum(vec.mat[,2]==3) # previtellogenic
s4<-sum(vec.mat[,2]==4) # vitellogenic
s5<-sum(vec.mat[,2]==5) # diapausing
gen0<-sum(vec.mat[,1]==0) # overwintering adult
gen1<-sum(vec.mat[,1]==1) # first generation
gen2<-sum(vec.mat[,1]==2) # second generation
n.adult<-sum(vec.mat[,2]==3)+sum(vec.mat[,2]==4)+sum(vec.mat[,2]==5) # sum of all adults
gen0.pop[day]<-gen0 # gen.0 pop size
gen1.pop[day]<-gen1
gen2.pop[day]<-gen2
S0[day]<-s0
S1[day]<-s1
S2[day]<-s2
S3[day]<-s3
S4[day]<-s4
S5[day]<-s5
g0.adult[day]<-sum(vec.mat[,1]==0)
g1.adult[day]<-sum((vec.mat[,1]==1 & vec.mat[,2]==3) | (vec.mat[,1]==1 & vec.mat[,2]==4) | (vec.mat[,1]==1 & vec.mat[,2]==5))
g2.adult[day]<-sum((vec.mat[,1]==2 & vec.mat[,2]==3) | (vec.mat[,1]==2 & vec.mat[,2]==4) | (vec.mat[,1]==2 & vec.mat[,2]==5))



N.newborn[day]<-n.newborn
N.death[day]<-n.death
N.adult[day]<-n.adult
print(c(day,n,n.adult))
}   

proc.time() - ptm
dd.cum<-cumsum(dd.day)
save(dd.day,dd.cum,S0,S1,S2,S3,S4,S5,N.newborn,N.death,N.adult,tot.pop,gen0.pop,gen1.pop,gen2.pop,g0.adult,g1.adult,g2.adult,file="opt$output")