   function [s, c] = jacrot(A, p, q)
%% function [s, c] = jacrot(A, p, q)
%% determine jacobi rotation for rows/columns p & q
   if (A(p,q) == 0) 
       c = 1; s =0.0;
   else
       tau =  0.5*(A(q,q) - A(p,p))/A(p,q);
       %% choose smallest angle 
       if (tau>0.0)
	 t = 1.0/(tau + sqrt(1+tau^2));
       else
	 t = 1.0/(tau - sqrt(1+tau^2));
       end
       c = 1/sqrt(1+t^2);
       s = t*c;
   end 
   
