generations <- 0:100
p <- 16/22
q <- 1 - p

AA <- p^2
Aa <- 2 * p * q
aa <- q^2

plot(generations, rep(p, length(generations)),
     type = "l",
     lwd = 2.5,
     col = "yellowgreen",  # "#3B6D11",
     ylim = c(0, 1),
     xlab = "generation",
     ylab = "allele frequency (p)",
     las = 1,
     cex = 1.5,
     cex.axis = 1.5
     )

polygon(c(generations[1], generations, generations[length(generations)]),
        c(0, rep(p, length(generations)), 0),
        col = "#a1d83a", border = NA)

polygon(c(generations[1], generations, generations[length(generations)]),
        c(1, rep(p, length(generations)), 1),
        col = "#fff371", border = NA)

lines(generations, rep(p, length(generations)),
      lwd = 2.5, col = "#3B6D11")

lines(generations, rep(AA, length(generations)),
      lwd = 1, lty = 2, col = "#444441")

lines(generations, rep(AA + Aa, length(generations)),
      lwd = 1, lty = 2, col = "#444441")

#text(max(generations) + 0.2, AA / 2,          "AA", adj = 0, cex = 1.2)
#text(max(generations) + 0.2, AA + Aa / 2,     "Aa", adj = 0, cex = 1.2)
#text(max(generations) + 0.2, AA + Aa + aa / 2,"aa", adj = 0, cex = 1.2)

cols <- hcl.colors(n_sims, palette = "Dark 3")

for (i in seq_len(n_sims)) {
  traj <- wf_sim(p, N, max(generations))
  lines(generations, traj,
        col = adjustcolor(cols[i], alpha.f = 0.9),
        lwd = 2)
}

npop=10; y=c()
for (i in 1:10000)
{
    x<-sample(c("A", "a"), size = npop, prob = c(p,q), replace = T)
    x<-sum(x=="A")/npop
    y=c(y,x)
}
hist(y)