Estimating a fuzzy RDD (rddtools)

I am trying to estimate a fuzzy RDD, so far I have used the rddtools package (but would be open to any package that can help me with this) by means of 2SLS.

My first regression Formula would look like this:

Treat = Alpha + betaDummy + f(Force) + errorterm

and then

Growth = Theta + RhoTreat + f(Force) + errorterm

In both equations, f(force) needs to be a pth-order parametric polynomial function. In the first step, treatment needs to be estimated by means of a probit regression and then enters the final regression as usual as far as I can tell. I am quite lost with how to do this as I am afraid that I don't fully understand how much of the process is already automatically done by R.

I tried to change below's Code (procided in rddtools) but only received various Errors.

reg_bin_glm <- rdd_gen_reg(rdd_object=house_rdd, fun= glm, family=binomial(link='probit'))

Do I need two separate equations for this in R ?

I am an absolute beginner when it comes to this and am quite lost.
Thanks!