M = 1;
Q = sqrt(1);
a = sqrt(0);
horizon = M + sqrt(M*M - Q*Q - a*a);
dr = 0.2;
r = [horizon+dr:dr:10*horizon];

delta = r.*r - 2*M*r + a*a + Q*Q;
sigmasq = (r.^2 + a^2).^2 - a^2 * delta;
pomega = sqrt(sigmasq) ./ r;

dzdpomega = sqrt((sigmasq ./ delta) .* (r.^3 ./ (r.^4 - a^2*M*r + a^2*Q^2)).^2 - 1);
dpomega = dr * (r.^4 - a^2*M*r + a^2*Q^2) ./ (sqrt(sigmasq) .* r.^2);

z = cumsum(dzdpomega .* dpomega);


theta = [0:pi/16:2*pi];

rr  = transpose(r) * ones(size(theta));
TT  = transpose(ones(size(r))) * theta;
zz  = transpose(z) * ones(size(theta));

rrr = transpose(r([length(r):-1:1])) * ones(size(theta));
rzz = transpose(z([length(z):-1:1])) * ones(size(theta));

rr = [ rrr; rr];
TT = [  TT; TT];
zz = [-rzz; zz];


pomegamatrix  = transpose(pomega)                        * ones(size(theta));
rpomegamatrix = transpose(pomega([length(pomega):-1:1])) * ones(size(theta));
pomegamatrix = [rpomegamatrix; pomegamatrix];

xx = pomegamatrix .* cos(TT);
yy = pomegamatrix .* sin(TT);

surf(xx,yy,zz);
%title('M = 1, Q = sqrt(1/2), a = 0');
