Is there an algorithm to compute the phase for a single frequecy?
The problem can be formulated as (nonlinear) least-squares problem:

$$F(\phi) = \frac{1}{2}\sum^{n}_{i=1}\left[ A \cdot \sin(\omega i+\phi) - f_{i}(\omega) \right]^{2}$$

where $F(\phi)$ is the objective function to minimize with respect to $\phi$.

The derivative is very simple:

$$F'(\phi)=\sum^{n}_{i=1} A\cdot \cos(\omega i+\phi)  \left[A \cdot \sin(\omega i+\phi) - f_{i}(\omega)\right]$$

The above objective function can be minimized iteratively using Gradient descent method (first order approximation), Newton method, Gauss-Newton method or Levenberg-Marquardt method (second order approximation - $F''(\phi)$ need to be provided in these).

Obviously, the above objective function has multiple minima because of periodicity, hence some penalty term can be added to discriminate other minima (for example, adding $\phi^{2}$ to the model equation). But I think the optimization will just converge to the nearest minima and you can update the result subtracting $2\pi k, k\in N  $.