hold_fpop = function(pop, cap) max(0, (pop - cap) / (1 - fgr))
hold_lazy = function(pop, cap) max(0, ceiling((pop - cap) / 120000) * 120000)

# Strategies - functions that return freighter pop given total pop (and previous turns freighter pop)
baseline = function(pop, lf) 0
hold33 = function(pop, lf) hold_fpop(pop,  pcap * .33)
hold25 = function(pop, lf) hold_fpop(pop,  pcap * .25)
hold42 = function(pop, lf) hold_fpop(pop,  pcap * .42)
hold50 = function(pop, lf) hold_fpop(pop,  pcap * .5)
hold33_lazy = function(pop, lf) hold_lazy(pop, pcap*.35)
hold42_lazy = function(pop, lf) hold_lazy(pop,  pcap * .42)
hold50_lazy = function(pop, lf) hold_lazy(pop,  pcap * .50)


hold33_50 = function(pop, lf) {
  if (pop >= pcap * .5) return(hold_fpop(pop, pcap * .5))
  if (pop >= pcap * .33) return(hold_fpop(pop, pcap * .33))
  return(0)
}

hold33_50_lazy = function(pop, lf) {
  if (pop >= pcap * .5) return(hold_lazy(pop, pcap * .5))
  if (pop >= pcap * .33) return(hold_lazy(pop, pcap * .33))
  return(0)
}

hold33_50_lazier = function(pop, lf) {
  if (pop >= pcap * .5) return(max(lf, hold_lazy(pop, pcap * .5)))
  if (pop >= pcap * .33) return(hold_lazy(pop, pcap * .33))
  return(0)
}

hold33_grow_50 = function(pop, lf) {
  if (pop < pcap * .33) return(0)
  if (pop > pcap * .5) return(lf)
  return(hold_fpop(pop, pcap*.33))
}
hold33_grow_50_lazy = function(pop, lf) {
  if (pop < pcap * .33) return(0)
  if (pop > pcap * .5) return(lf)
  return(hold_lazy(pop, pcap*.33))
}


# Code for simulating growth
strategy = "hold33_50"
overflow = T
pgr = .19
obrm = T
divisor = 1000
fgr = .5 * pgr

sim = function(strategy, target) {
  pop = 125000
  cumres =0
  result = NULL
  lagfreight = 0
  for (y in 2400:2450) {

    row = data.frame(year=y, pop=pop)
    row$freight=ifelse(pop >= pcap, 0, strategy(pop, lagfreight))
    row$planet = pop - row$freight
    row$fgrowth = row$freight * fgr
    row$after.overflow = row$planet + ifelse(overflow, row$fgrowth, 0)
    row$cap = row$after.overflow / pcap
    row$crowd = ifelse(row$cap>0.25, 16/9 * (1-row$cap)^2,1)
    row$pgrowth = row$after.overflow * pgr * row$crowd
    row$res = floor(row$after.overflow / divisor)
    cumres = cumres + row$res
    row$cumres = cumres
    row$newpop = row$pop + row$fgrowth + row$pgrowth
    pop = min(pcap, row$newpop)
    lagfreight = row$freight
    result = rbind(result, row)
  }
  result
}

# Run the simulations

comparison = NULL; alldata = NULL
for (hab in c(1.0, 0.5)) {

  for (strategy in c("baseline", "hold25", "hold33", "hold42", "hold33_lazy", "hold42_lazy", "hold50_lazy",
                     "hold33_50", "hold33_50_lazy", "hold33_50_lazier")) {

    pcap = 1000000 * ifelse(obrm, 1.1, 1) * hab
    target = pcap

    result = sim(get(strategy), target)
    fullpop = if (max(result$pop >= target)) min(result$year[result$pop >= target]) else NA

    row = data.frame(hab=hab, strategy=strategy,  fullpop=fullpop,
                     res10=result$cumres[result$year == 2410],
                     res15=result$cumres[result$year == 2415],
                     res20=result$cumres[result$year == 2420],
                     res30=result$cumres[result$year == 2430])
    comparison = rbind(comparison, row)
    result$strategy = strategy
    result$hab = hab
    alldata = rbind(alldata, result)
  }
}
comparison %>% arrange(hab, strategy)

# Analyses and plots

alldata = alldata %>% mutate(popgrowth=pgrowth+fgrowth)

alldata %>% filter(strategy %in% c("baseline", "hold33", "hold33_lazy")) %>% ggplot + geom_line(aes(x=year, y=pop, colour=strategy))
alldata %>% filter(year < 2425, strategy %in% c("hold33", "hold33_50", "hold33_grow_50")) %>% ggplot + geom_line(aes(x=year, y=cumres, colour=strategy))
alldata %>% filter(year < 2425, strategy %in% c("hold33", "hold33_50", "hold33_grow_50")) %>% ggplot + geom_line(aes(x=year, y=res, colour=strategy))

alldata %>% filter(year < 2420, strategy %in% c("hold33_lazy", "hold42_lazy", "hold50_lazy", "hold33_50_lazier")) %>% ggplot + geom_line(aes(x=year, y=cumres, colour=strategy))
alldata %>% filter(year < 2420, strategy %in% c("hold33_lazy", "hold42_lazy", "hold50_lazy", "hold33_50_lazier")) %>% ggplot + geom_line(aes(x=year, y=res, colour=strategy))
alldata %>% filter(year < 2420, strategy %in% c("hold33_lazy", "hold42_lazy", "hold50_lazy", "hold33_50_lazier")) %>% ggplot + geom_line(aes(x=year, y=popgrowth, colour=strategy))
alldata %>% filter(year < 2420, strategy %in% c("hold33_lazy", "hold42_lazy", "hold50_lazy", "hold33_50_lazier")) %>% ggplot + geom_line(aes(x=year, y=pop, colour=strategy))

alldata %>% filter(year < 2420 & hab==0.5, strategy %in% c("hold33_lazy", "hold42_lazy", "hold50_lazy", "hold33_50_lazier")) %>% ggplot + geom_line(aes(x=year, y=cumres, colour=strategy))
alldata %>% filter(year < 2420 & hab==0.5, strategy %in% c("hold33_lazy", "hold42_lazy", "hold50_lazy", "hold33_50_lazier")) %>% ggplot + geom_line(aes(x=year, y=res, colour=strategy))

par(mfrow=c(1,2))
alldata %>% filter(hab==1.0, year < 2418, strategy %in% c("hold33", "hold50_lazy", "hold33_50_lazier")) %>% ggplot + geom_line(aes(x=year, y=cumres, colour=strategy))
alldata %>% filter(hab==1.0, year < 2418, strategy %in% c("hold33", "hold50_lazy", "hold33_50_lazier")) %>% ggplot + geom_line(aes(x=year, y=res, colour=strategy))


alldata %>% filter(year < 2425, strategy %in% c("hold33_lazy", "hold33_50_lazy", "hold33_50_lazier")) %>% ggplot + geom_line(aes(x=year, y=res, colour=strategy))


alldata %>% filter(year < 2425, strategy %in% c("baseline", "hold25", "hold33", "hold42")) %>% ggplot + geom_line(aes(x=year, y=cumres, colour=strategy))
alldata %>% filter(year < 2425, strategy %in% c("baseline", "hold25", "hold33", "hold42")) %>% ggplot + geom_line(aes(x=year, y=res, colour=strategy))

alldata %>% filter(strategy %in% c("baseline", "hold33", "hold33_lazy")) %>% ggplot + geom_line(aes(x=year, y=pop, colour=strategy))
alldata %>% filter(strategy %in% c("baseline", "hold33_lazy")) %>% ggplot + geom_line(aes(x=year, y=pop, colour=strategy))

alldata %>% filter(strategy == "hold33_grow_50_lazy")
alldata %>% filter(strategy %in% c("hold50_lazy"), hab==1.0, year < 2420)
alldata %>% filter(strategy == "hold33_lazy", hab=.5)
alldata %>% filter(strategy == "hold33_50")