Which formulas are you referring to by "sin" and "true cos"? If you refer to the formulas on my previous post, you will see that I don't know of any formula for the triplex that produce abs(power) roots. Are you adding multiples of 2pi to both rotational angles theta and phi, or just phi, or are you doing something else?
Hi Paul, my previous post was hasty due to rushing out. See if this makes sense.
I only use the trig formula, basically as yours but with three variations of phi handling. I had difficulty with the trig version roots initially, because treating the roots independently is not quite reversible due to spherical coordinates being a double cover if phi goes beyond pi (which it does when calculating powers greater than 1). To get equivalence to a forward iteration method requires some parity correction by pi to get the root in the right hemisphere of theta. It took some experiment and thought to get the two correct variations working.
When I wrote "using Abs(power) in appropriate places," I meant in setting the rootlimit and number of branches, and the multiplier for dp, so those all handle negative powers as is. The trig formulas don't need any special handling for negative power.
In my code I've dropped the root validity checking as I've given up on fractional powers, but rootlimit is still used for
some angle choices because it works. The following is pseudocode for clarity. (Delphi uses := for assignment, and I have some optimizations such as precomputing invpower, and computing sin and cos in a single sincos function.)
So I have, in all cases:
NumBranches = ceil(abs(power));
Numchoices = sqr(numbranches);
// ktheta and kphi are independent choices, from 0 to numbranches-1
r = sqrt(sqr(x)+sqr(y)+sqr(z));
theta= theta + 2pi*(ktheta - floor((theta+rootlimit)/2pi);
then, for standard triplex:
theta= theta + kphi*pi; // a correction to get the parity of theta right
phi = arcsin(z/r);
if phi>0 then phi= phi - kphi*pi else phi= phi + kphi*pi;
// but for the inverse, z=-sin(phi) doesn't work. To get both hemispheres, need to compare new and old z
if newz*z>0 then newz= -newz;
or, for my variation of the z=cos(phi) triplex (the simple one with some flipped roots):
phi = arccos(z/r);
phi= phi + 2pi*(kphi - floor((phi+rootlimit)/2pi);
or, for the "true" z=cos(phi) triplex:
theta= theta + kphi*pi;
phi = arccos(z/r);
phi= phi + 2pi*kphi;
There may be a purer way to get the standard triplex working but I tried hard without success, so if you do find one, please post it.