M = 1;
Q = 0;
da = 0.1;
a = [0:da:sqrt(M^2-Q^2)-da];
horizon = M + sqrt(M*M - Q*Q - a.*a);

dr = 0.03;
r = [M:dr:5*M];

rr = transpose(r) * ones(size(a));
aa = transpose(ones(size(r))) * a;

horizonmatrix = transpose(ones(size(r))) * horizon;
valid = rr >= horizonmatrix+dr;
rr = rr .* valid + (horizonmatrix+dr) .* (1-valid);

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

dzdpomega = sqrt((sigmasq ./ delta) .* (rr.^3 ./ (rr.^4 - (aa.^2).*(M*rr) + (aa.^2)*Q^2)).^2 - 1);
dzdpomega = dzdpomega .* valid;

dpomega = dr * (rr.^4 - (aa.^2).*(M*rr) + (aa.^2)*Q^2) ./ (sqrt(sigmasq) .* rr.^2);

z = cumsum(dzdpomega .* dpomega, 1);


surf(pomega,aa,z);
xlabel('\varpi');
ylabel('a');
zlabel('z');
