I don’t see the pattern for the (etc). If x=x_c then from (1) we have
y_0 = f(x_c)\sin(\theta_1)\cos(\phi_1)
but from (2) we have
y_0 = f(x_c)[1 + \sin(\theta_1)\cos(\phi_1) - \sin(\theta_2)\cos(\phi_2)]
If the formula carries for the other components it looks like the magnitude is discontinuous?
I think your intuition is right to want to use a spline model (parameterized on the logit/probit scale if 0 \leq f \leq 1). And when you say “covariance”, do you mean that the error terms covary
y | x = f(x, \theta, \phi) + \epsilon
with \mathrm{cov}[\epsilon_i, \epsilon_j] \neq 0
Or simply that even though the errors are independent \mathrm{cov}[\epsilon_i, \epsilon_j] =0 the marginal covariance is not \mathrm{cov}[y_i, y_j] \neq 0?