function riemann_sphere(z0)

hold on;

%------------------------------------------------------------
% Draw the sphere
%------------------------------------------------------------
global R;
R = 1;
[x y z] = ellipsoid(0,0,R,R,R,R);
mesh(x,y,z);
colormap(gray);
axis([-5 5 -5 5 -4 5]);
grid on;
xlabel('x');
ylabel('y');

%------------------------------------------------------------
% Draw the projection
%------------------------------------------------------------
plot(z0,'*r');
[x1 y1 z1] = riemann_map(z0)
plot3(x1,y1,z1,'*r');

% Inverse
w0 = 1/z0;
plot(w0,'*');
[x2 y2 z2] = riemann_map(w0)
plot3(x2,y2,z2,'*');
hold off;

%------------------------------------------------------------
% Easier to view on the plane
%------------------------------------------------------------
figure;
[x y z] = ellipsoid(0,0,0,R,R,0);
mesh(x,y,z);
hold on;
plot(z0,'*r');
plot(x1,y1,'*r');
plot(w0,'*');
plot(x2,y2,'*');
hold off;
grid on;

%------------------------------------------------------------
% Riemann mapping
%------------------------------------------------------------
function [x y z] = riemann_map(z0)

global R;

theta = atan(imag(z0)/real(z0));
s     = abs(z0);

r = 4*R^2*s/(s^2+4*R^2);
h = 2*R*s^2/(s^2+4*R^2);

x = r*cos(theta);
y = r*sin(theta);
z = h;