M = 1;
Q = 0;
da = 0.1;


pomegamatrix = [];
zmatrix = [];
for a = 0:da:sqrt(M*M-Q*Q),


horizon = M + sqrt(M*M - Q*Q - a*a);
dr = 0.001 * horizon;
r = [horizon+dr:dr:horizon+2000*dr];

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);
plot(pomega,z);
xlabel('\varpi');
ylabel('z(\varpi)');



pomegamatrix = [pomegamatrix ; pomega];
zmatrix = [zmatrix ; z];
end

plot(transpose(pomegamatrix),transpose(zmatrix));
xlabel('\varpi');
ylabel('z(\varpi)');
