function [stat p]=circ_ttest_p_second(alpha_1,k_1,alpha_2,k_2)
%[stat p]=circ_ttest_np_second(alpha_1,k_1,alpha_2,k_2)
%parametric paired ttest for second order analysis
%see Zar,1999,, page 641
%written by Andre M. Cravo 07/Sep/2011

%alpha_1 and alpha_2 are vector with circular values
%k_1 and k_2 are the mean resultant lentgh for each

%here are just some tests to see if the script works
% alpha_1=circ_vmrnd(0,2,50);
% k_1=rand(50,1);
% alpha_2=circ_vmrnd(pi/2,2,50);
% k_2=rand(50,1);


 %equations 27.37, page 646)
X=(k_2.*(cos(alpha_2)))-(k_1.*(cos(alpha_1)));
Y=(k_2.*(sin(alpha_2)))-(k_1.*(sin(alpha_1)));
   
%equations 27.24 to 27.26 (page 638)
k=numel(alpha_1);
sum_xsquare=[sum((X.^2))]-[((sum(X))^2)/k];
sum_ysquare=[sum((Y.^2))]-[((sum(Y))^2)/k];
sum_xy=[sum(X.*Y)]-[(sum(X)*sum(Y))/k];

%equations 27.27 (page 639)
k_factor=(k*(k-2))/2;
mean_X=mean(X);
mean_Y=mean(Y);
mean_X_sq=mean_X^2;
mean_Y_sq=mean_Y^2;

F_num=[(mean_X_sq*sum_ysquare)-(2*mean_X*mean_Y*sum_xy)+(mean_Y_sq*sum_xsquare)];
F_den=[(sum_xsquare*sum_ysquare)-(sum_xy)^2];

stat=k_factor*(F_num/F_den);
 
%critical value one-tailed F test wit df=2, k-2, App22 (Zarr)
p = 1-fcdf(stat,2,k-2);

end