% We give two points that satisfy the requirements, i.e. that are within
% the unit circle and not on the same line through the origin

isOctave = exist('OCTAVE_VERSION', 'builtin');

% Set points A and B
A = [0.5 0.7]';
B = [-0.1 -0.8]';

if rank([A, B]) < 2
	error('Requirements not fulfilled')
end

%% --------------------------------------------------------------------------------------------------
% Project points A and B onto the unit sphere
% Note: for each of the points, there are two possibilities, but we have
% only to distinguish between projected points on the same side of the sphere
% and those on opposite sides.
x = A;
x(3) = sqrt(1 - A(1)^2 - A(2)^2);

y = B;
y(3) = sqrt(1 - B(1)^2 - B(2)^2);

z = B;
z(3) = - sqrt(1 - B(1)^2 - B(2)^2);

% Calculate normal vector n that is perpendicular to points x and y, i.e.
% the n is a normal vector of the plane through the origin and points x and y
n = cross(x,y);
n = n / norm(n);

% Determine the parameters of a 2d cirle in the 3 dimensional space
t1 = cross(n,x);
t2 = cross(n,t1);

% Determine points along the cirle in 3d
t = linspace(0,2*pi,200);
Y = t1*cos(t) + t2*sin(t);

% Do the same for the point on the opposite site of the sphere
n = cross(x,z);
n = n / norm(n);

% Determine the parameters of a 2d cirle in the 3 dimensional space
t1 = cross(n,x);
t2 = cross(n,t1);

% Determine points along the cirle in 3d
Z = t1*cos(t) + t2*sin(t);
%% --------------------------------------------------------------------------------------------------

% Make 2d plot
figure(1), clf

% Plot unit circle
plot(cos(t),sin(t),'k-'), hold on

axis(1.1*[-1 +1 -1 1])
axis equal

% Plot points A and B
plot(A(1), A(2), 'ko', 'MarkerFaceColor', 'r')
plot(B(1), B(2), 'ko', 'MarkerFaceColor', 'r')

% Plot origin
plot(0, 0, 'k+')

% Plot the projection of the circles in 2d
plot(Y(1,:),Y(2,:),'b--')
plot(Z(1,:),Z(2,:),'b--')

% Make 3d plot
figure(2), clf

% Plot unit sphere
[sX,sY,sZ] = sphere(50);
if isOctave
	h = mesh(sX,sY,sZ);
else
	h = surf(sX,sY,sZ);
	set(h, 'FaceAlpha', 0.5, 'EdgeColor', 'none', 'LineStyle', 'none')
end
hold on

% Plot points x, y and z
plot3(x(1),x(2),x(3), 'ko', 'MarkerFaceColor', 'r')
plot3(y(1),y(2),y(3), 'ko', 'MarkerFaceColor', 'r')
plot3(z(1),z(2),z(3), 'ko', 'MarkerFaceColor', 'r')

% Plot origin
plot3(0,0,0,'k+')
axis(1.1*[-1 +1 -1 1 -1 1])

% Plot the circles in 3d
plot3(Y(1,:),Y(2,:),Y(3,:))
plot3(Z(1,:),Z(2,:),Z(3,:))
