It looks like what you are looking for is a threshold model, also called segmented regression, or change-point model (also the latter refers more to time series settings). In a threshold model, parameters change depending on the threshold variable (I guess your $\theta$) are within some ranges.
Estimation
To simplify, let us denote by $x$ the covariate that is also threshold variable, by $x_0$ the transition/threshold point. To estimate your mode, you can do a conditional/concentrated estimation: the idea is that, conditional on a given $x_0$, you can use the standard procedure to estimate your parameters $\alpha_0(x_0)$, $\beta_0(x_0)$, and obtain a cost function (sum of squares, likelihood...), let's write it $SSR(x_0)$. Now you simply look for the $x_0$ that minimizes that cost function. Given that the regime indicator is a discontinuous function, you can not use standard algorithms. A simple solution is to run a grid search over possible values. As you are looking for multiple threshold, you will have a 2, 3-dimensional grid. To avoid the dimensionality problem, you could run a conditional grid search, looking first for only one threshold, there for a second given the first one, eventually iterating, etc.
In steps:
Write your cost function, given regression parameters $\alpha^1$, $\beta^1$,$\theta_0^1$, $\alpha^2$, $\beta^2$,$\theta_0^2$ and threshold value $x_0$.
Write a function optimizing step 1, with a given threshold value. I.e. find $\hat{\alpha}^1(x_0)$, $\hat{\beta}^1(x_0)$,$\hat{\theta}_0^1(x_0)$, $\hat{\alpha}^2(x_0)$, $\hat{\beta}^2(x_0)$,$\hat{\theta}_0^2(x_0)$, given $x_0$.
Select a range of candidates for $x_0$. Typically, one sort x, take unique values, and keep values from the 10% to 90% quantiles. You can use less than the full values obviously. use however existing ones (i.e., if only 2.1 and 2.2 are observed, using 2.01, 2.02 etc as candidates will give same result).
Run the function of step 2 for every possible candidates
Pick the x value that gives smallest cost function in step 4. This is your estimate of $\hat{x}_0$
Use the function in 2, given now the optimal $\hat{x}_0$. This leads the optimal parameters. $\hat{\alpha}^1(\hat{x}_0)$,
$\hat{\beta}^1(\hat{x}_0)$, etc.
Find below a small code as example, in R:
n <- 500
x <- sort(rnorm(n, mean=5, sd=2))
fo <- function(x, alpha, beta, theta) alpha -beta / cos(theta-x)
df <- data.frame(x=x, y=ifelse(x<5, fo(x, 1, 3, 5), fo(x, 2, 4, 5))+rnorm(n, sd=10))
plot(y~x, data=df, type="l")
## return sum of squares: simple and with regimes (2 here)
nls_ssr <- function(par, x, y) sum( (y-fo(x, par[1], par[2], par[3]))^2 )
nls_ssr_2reg <- function(par, x, y, thresh) {
ssr_reg1 <- nls_ssr(par[1:3],x[x<thresh], y[x<thresh])
ssr_reg2 <- nls_ssr(par[4:6],x[x>=thresh], y[x>=thresh])
ssr_reg1+ssr_reg2
}
## grid search for all regimes
x_search <- x[ (n*0.1):(n*0.9)]
ssrs <- sapply(x_search, function(thresh) optim(c(2,3,5, 2, 3, 5), nls_ssr_2reg, y=df$y, x=df$x, thresh=thresh)$value)
#3 select best, re-estimate
x_search[which.min(ssrs)]
optim(c(2,3,5, 2, 3, 5), nls_ssr_2reg, y=df$y, x=df$x, thresh=x_search[which.min(ssrs)])
This gives regression parameters values of 1.719521 3.147032 5.000069 1.136988 4.053732 4.999990, and threshold value 6.5, not too far from the true values 1,3,5,2,4,5 and 5
