compute_F.m

function F = compute_F(alpha,B,q1,q2,q3)

% r = compute_r(alpha,B,q1,q2,q3);

v = [cos(alpha);sin(alpha)];

b1 = B(1,1);
b2 = B(2,2);
b3 = B(1,2);

%compute Q(v_alpha x v_alpha)
Q = q1.*v(1).^4 + q2.*v(2).^4 + 2.*q3.*v(1).^2.*v(2).^2;
% 
% TP = v*v';
% L = stVenant(TP,mu,lambda);

tmp1 = q1.*(v(1).^2+b1).^2 + q2.*(v(2).^2+b2).^2 + q3.*(sqrt(2)*v(1)*v(2)+b3).^2;
tmp2 = q1.*b1.^2 + q2.*b2.^2+ q3.*b3.^2;
L = 0.5*(tmp1-Q-tmp2) ;       %Polarization identity 


r = L./Q;


% F = r.^2.*Q - 2.*r.*trace(L'*B)
F = (r.^2).*Q - 2.*r.*L;


end