### find t of intersection point of circle with radius r with the bezier curve
library(bezier)
library(nleqslv)
## inits
y0 <- function(t) (1-t)^3
y1 <- function(t) 3*(1-t)^2*t
y2 <- function(t) 3*(1-t)*t^2
y3 <- function(t) t^3
r <- 0.2
P0x <- P0y <- P3y <- 0
P1x <- 0.4
P1y <- 0.9
P2x <- 0.8
P2y <- 0.4
P3x <- 1
## functions for x and y components of B(t)
b1 <- function(t) y0(t)*P0x + y1(t)*P1x + y2(t)*P2x + y3(t)*P3x
b2 <- function(t) y0(t)*P0y + y1(t)*P1y + y2(t)*P2y + y3(t)*P3y
## equation to solve is:
## b1^2 + b2^2 = r^2
eq <- function(t) b1(t)^2 + b2(t)^2 - r^2
## find root (with Broyden (or Newston))
t.r <- nleqslv(0.2, eq)
t.r$x
## [1] 0.0746813
## plot bezier
t <- seq(0, 1, length=100)
p <- matrix(c(P0x,P0y,P1x,P1y,P2x,P2y,P3x,P3y), ncol=2, byrow=TRUE)
be <- bezier(t, p)
plot(be, xlim=c(0,1), ylim=c(0,1), type="l")
## plot 1/4-circle
rx <- r * sin(seq(0, pi/2, length=100))
ry <- r * cos(seq(0, pi/2, length=100))
lines(rx,ry)
## plot intersection point
bb <- bezier(t.r$x, p)
points(bb)
