- alpha = atan2(ys - tan(beta)*ze/ys, xs);
- delta = asin(sin(beta)*ye/ys + cos(beta)*ze);
-
- tmp = alpha - (SGD_2PI/24)*gst;
- if (tmp < -SGD_PI) {
- do tmp += SGD_2PI;
- while (tmp < -SGD_PI);
- } else if (tmp > SGD_PI) {
- do tmp -= SGD_2PI;
- while (tmp < -SGD_PI);
- }
+ double ra = atan2(ye, xs);
+ double dec = atan2(ze, sqrt(xs * xs + ye * ye));
+
+ tmp = ra - (SGD_2PI/24)*gst;
+
+ double signedPI = (tmp < 0.0) ? -SGD_PI : SGD_PI;
+ tmp = fmod(tmp+signedPI, SGD_2PI) - signedPI;