Finding PHASE between two sinusoids with FFT/DFT (ArduinoFFT w/ C)
You can't do better than this for accuracy in my experience.

A Two Bin Solution

Select a frame size of 1 1/2 cycles or 2 1/2 cycles.
You only need to calculate two bins, either 1 and 2, or 2 and 3.  No need for a full DFT.  So it is very efficient as well.  By placing the frequency squarely near the center of the two bins you minimize noise effects and avoid needing to deal with near bin frequencies which call for a different set of equations for the fitting vector.
The alternative form is for when the frequency is really close to a bin.

DFT Bin Value Formulas for Pure Real Tones
An Alternative Form of the Pure Real Tone DFT Bin Value Formula

There are other tricks to do it point by point.  I'll let others provide those.