M = 1;
dQ = 0.1;
Q = [0:dQ:1-dQ];
horizon = M + sqrt(M*M - Q.*Q);

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

rr = transpose(r) * ones(size(Q));
QQ = transpose(ones(size(r))) * Q;

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

dzdr = sqrt((2*M*rr - QQ.*QQ) ./ (rr.*rr - 2*M*rr + QQ.*QQ));
dzdr = dzdr .* valid;
z = cumsum(dzdr * dr, 1);

surf(rr,QQ,z);
xlabel('r');
ylabel('Q');
zlabel('z');
