### 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)