
function [s, c] = getrot(a, b, d);

% function [s, c] = getrot(a, b, d);
%
% calculate sine and cosine of a rotation to zero symmetric
% off-diagonal matrix elements; i.e., solve for s and c, in
% terms of a, b, and d, so that the off-diagonal elements of
% the product
%
%           | c -s | | a  b | | c  s |
%   J'A J = |      |*|      |*|      |
%           | s  c | | b  d | |-s  c |
% 
% are zero.

disc = sqrt(a^2 - 2*a*d + d^2 + 4*b^2);

t1 =  (a - d + disc) / (2*b);

t2 =  (a - d - disc) / (2*b);

if abs(t1) < abs(t2)
  t = t1;
else
  t = t2;
end

% return sine and cosine

c = 1 / sqrt(1 + t^2);

s = c * t;

