For the phase between 2 blocks of sinusoidal data, block1 and block2:
(using classical phase detector mapping)
Phase in degrees= 90- 90*mean(sign(block1-mean(block1))* sign(block2-mean(block2)))
Note the above is for phase angles between 0 and 180 degrees.
For angles between 180 and 360 degrees:
Phase in degrees = 270 + 90*mean(sign(block1-mean(block1))* sign(block2-mean(block2)))
This approach converts the sine waves to normalized square waves (value of +/-1). If you multiple two square waves in the time domain, the time domain average will be linearly proportional to the phase between the square waves.
The plot shows truth vs an estimate generated from the scilab
code below.

scilab
code for implementing this idea
//8673
T=1000;
t = [0:T-1];
omega = 2*%pi*0.398374982349;
phi1 = 0;
x1 = sin(omega*t + phi1);
phase_est = [];
phases = -180:180;
for phase_in_degrees = phases,
phi2 = 2*%pi*phase_in_degrees/360;
x2 = sin(omega*t + phi2);
if (phase_in_degrees > 0)
phase_est = [phase_est; phase_in_degrees (90 - 90*mean(sign(x1 - mean(x1)).*sign(x2 - mean(x2))))]
else
phase_est = [phase_est; phase_in_degrees (-90 + 90*mean(sign(x1 - mean(x1)).*sign(x2 - mean(x2))))]
end
end
clf
subplot(211)
plot(phases, phase_est')
title('Plot of true phase and estimated phase')
subplot(212)
plot(phases, phase_est(:,1) - phase_est(:,2))
title('Plot of error between truth and estimate')