xkcd 135 #2

solve_xkcd135b <- function(h=0.001) {
    # start positions
    start_rtop <- c(0, sqrt(20^2-10^2))
    start_rright <- c(10, 0)
    start_rleft <- c(-10,0)
    start_me <- c(0, sqrt(3)/6*20)

    # survival time, given angle phi and stepsize h
    pos_me <- function(t, phi) start_me + c(sin(phi)*6*t, cos(phi)*6*t)
    velo_rtop <- function(t) pmin(4*t, 10) # wounded top raptor
    velo_rapt <- function(t) pmin(4*t, 25) # left & right
    dist_rapt <- function(from, to, velo) integrate(velo, from, to)$value
    surv_time <- function(phi, h, pos0, velo) {
        curr <- pos0
        t <- 0
        repeat {
            t <- t + h
            old <- curr
            direction <- c(pos_me(t-h, phi)-old)
            len <- dist_rapt(t-h, t, velo)
            if(dist(rbind(pos_me(t, phi), old))<len) return(t) # gotcha
            curr <- old + direction*len/sqrt(sum(direction^2))
        }
        return(curr)
    }
    surv3_time <- function(phi, h) min(c(
        surv_time(phi, h, start_rtop, velo_rtop),
        surv_time(phi, h, start_rleft, velo_rapt),
        surv_time(phi, h, start_rright, velo_rapt)
    ))

    # find the best angle
    res <- optimize(surv3_time, c(0, 2*pi), h=h, maximum=TRUE)
    res_angle <- res$maximum
    res_time <- res$objective

    # plot
    plot(1, col="white", xlim=c(-17, 11), ylim=c(-1, sqrt(20^2-10^2)+1), axes=FALSE, xlab="", ylab="")
    me <- Vectorize(pos_me, "t")(seq(0, res_time, length.out=100), res_angle)
    lines(me[1,], me[2,], col="red", lwd=2)

    pos_rapt <- function(t, phi, h, pos0, velo) {
        pos <- matrix(0, nrow=100, ncol=2)
        pos[1,] <- pos0
        t_vec <- seq(0, t, length.out=100)
        # browser(skipCalls=2)
        for (i in 2:100) {
            direction <- pos_me(t_vec[i-1], phi)-pos[i-1,]
            len <- dist_rapt(t_vec[i-1], t_vec[i], velo)
            pos[i,] <- pos[i-1,] + direction*len/sqrt(sum(direction^2))
        }
        return(pos)
    }
    rtop <- pos_rapt(res_time*0.97, res_angle, h, start_rtop, velo_rtop)
    rright <- pos_rapt(res_time*0.97, res_angle, h, start_rright, velo_rapt)
    rleft <- pos_rapt(res_time*0.97, res_angle, h, start_rleft, velo_rapt)

    lines(rtop[,1], rtop[,2], lwd=2, col="blue")
    lines(rright[,1], rright[,2], lwd=2, col="gray")
    lines(rleft[,1], rleft[,2], lwd=2, col="black")

    # answer
    return(res_angle)
}