function eyes = GetEyesGivenCV(I,Hcdeg,Vcdeg,Tcdeg,HDdeg,VDdeg,TDdeg)
% Given eye posture in cyclopean & vergence coordinates in deg, return an eyes
% structure containing everything we might need.

global L R CYC DELTA HALFDELTA ALL


eyes(ALL).I=I;
% Convert to radians:
eyes(CYC).H = Hcdeg * pi/180;
eyes(CYC).V = Vcdeg * pi/180;
eyes(CYC).T = Tcdeg * pi/180;
eyes(DELTA).H = HDdeg * pi/180;
eyes(DELTA).V = VDdeg * pi/180;
eyes(DELTA).T = TDdeg * pi/180;

angles = {'H' 'V' 'T'};
for j = 1:3
    % Get half-vergence
    eyes(HALFDELTA).(angles{j}) = eyes(DELTA).(angles{j})/2;
    % Get left-eye coordinate:
    eyes(L).(angles{j}) = eyes(CYC).(angles{j}) - eyes(HALFDELTA).(angles{j});
    eyes(R).(angles{j}) = eyes(CYC).(angles{j}) + eyes(HALFDELTA).(angles{j});
end

% Other things:
eyes(L).nodal = [I/2;0;0];
eyes(R).nodal = [-I/2;0;0];

% Left and right eye Rotation matrices:
for jeye=[L R]
    eyes(jeye).M = RotationMatrix3(eyes(jeye).V,[1;0;0]) * RotationMatrix3(eyes(jeye).H,[0;1;0]) * RotationMatrix3(eyes(jeye).T,[0;0;1]);
end
for j = 1:3
    eyes(CYC).M = 0.5*(eyes(L).M + eyes(R).M);
    eyes(DELTA).M = eyes(L).M - eyes(R).M;
    eyes(HALFDELTA).M = eyes(DELTA).M /2;
end
% Get rotation vectors:
for k=1:3
    for angle = [CYC DELTA HALFDELTA L R];
        eval(['eyes(angle).m' num2str(k) '= eyes(angle).M(:,k);']);
    end
end


function R = RotationMatrix3(theta,omega)
% function R = RotationMatrix3(theta,omega)
% constructs a 3d rotation matrix for a rotation through an angle theta (in radians) about an axis given by the vector omega
%
% Example: 
% suppose we have a vector v1 = [3;4;5]; and we want to rotate this vector through 30degrees about the axis [1 1 1].
% Then the rotated vector is given by v2 = R * v1; where R = RotationMatrix3(pi/6,[1 1 1]).
%
% Jenny Read 2/4/2003

if prod(size(theta))~=1
    disp('theta must be a scalar')
    return
end
if prod(size(omega))~=3
    disp('omega must be a 3-entry vector')
    return
end
%Normalise omega
omega = omega ./ sqrt(sum(omega.^2));
c = cos(theta);
s = sin(theta);
w1 = omega(1);
w2 = omega(2);
w3 = omega(3);
R = [    1+(1-c)*(-w3^2-w2^2), -s*w3+(1-c)*w2*w1,  s*w2+(1-c)*w3*w1 ; ...
        s*w3+(1-c)*w2*w1,    1+(1-c)*(-w3^2-w1^2), -s*w1+(1-c)*w3*w2 ; ...
        -s*w2+(1-c)*w3*w1,  s*w1+(1-c)*w3*w2,    1+(1-c)*(-w2^2-w1^2)];

