qtp1<-function(qn,k)
{
  # calculate frquency of q after n+1 generations 
  abo<-(1-k)*qn
  bel<-(1-k*qn)
  return(abo/bel)
}
 
qtn<-function(qn,k,n)
{
  # frequency of deleterious allele a after n generations; qn is the starting frequency
    qk<-qn
    for(i in 1:n)
    {
       qk<-qtp1(qk,k)
    }
    return(qk)
}
 
calculateD<-function(qn,k,n)
{
  # genetic load D ~ k*Sum(qn)
  D=0
  qk<-qn
  for(i in 1:n)
  {
    D<-D+qk*k
    qk<-qtp1(qk,k)
  }
  return(D)
}
 
plotqdecay<-function(startfreq)
{
  # plot the decay of the frequency q with a constant selection of k (=death rate)
  x<-seq(1,100)
  y<-sapply(x,qtn,qn=startfreq,k=0.1)
  plot(x,y,type="l")
}
 
plotload<-function()
{
  # frequency of deleterious allele
  x<-seq(0.0001,0.9999,0.01)
  y<-sapply(x,calculateD,k=0.001,n=10000)
  plot(x,y)
}