compute_r.m

function r = compute_r(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;


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;




end