Hi bil,
Here is a function that does that job. It appears that you want to do the division (a-i)/(a+i) before taking the square root. To put the branch cut on the +x axis, use
sqrtbc(0,(a-i)/(a+i))
It's also possible to do sqrt(a-i) / sqrt(a+i) in which case the function below can send the <branch cut from the brance point at i> and the <branch cut from the branch point at -i> in two different directions if needed.
eps = 1e-10;
% demonstrate branch cut along -x axis
sqrtbc(pi,-4+eps*i) -sqrtbc(pi,-4-eps*i)
ans = 0.0000 + 4.0000i
% demonstrate branch cut along +x axis
sqrtbc(0,4+eps*i) -sqrtbc(0,4-eps*i)
ans = 4
% demonstrate branch cut along +y axis
sqrtbc(pi/2,4i+eps) -sqrtbc(pi/2,4i-eps)
ans = 2.8284 + 2.8284i
function y = sqrtbc(theta,zarg)
% sqrt function with branch cut in zarg from 0 to infinity along a ray
% at angle theta (in radians) measured from the +x axis in the usual way,
% with -pi<=theta<=pi. theta = pi is the usual square root.
% for zarg on the +x axis, sqrt behavior is conserved,
% i.e. sqrtbc(theta,zarg) is positive and real for any theta.
%
% y = sqrtbc(theta,zarg)
if theta==0
phi = pi;
else
phi = theta -pi*sign(theta);
end
y = exp(i*phi/2)*sqrt(zarg*exp(-i*phi));
% translations: sqrtbc(theta, z-b) has branch cut in the z plane from
% branch point z = b out to infinity, along a ray at angle theta.
%
% for the usual square root with branch cut along -x,
% the real part of sqrt(z) is positive (or 0) for all z.
% for the modified square root with branch cut along +x,
% the imaginary part of sqrt(z) is positive (or 0) for all z.
