diff --git a/src/libnptm/clu_subs.cpp b/src/libnptm/clu_subs.cpp
index 42d136b28844f9ddb2a13f1909e355d1e742918e..78bfc4af3adbbc82f73387b8125347f34e6eaacd 100644
--- a/src/libnptm/clu_subs.cpp
+++ b/src/libnptm/clu_subs.cpp
@@ -1473,21 +1473,25 @@ void polar(
   bool onx = (y == 0.0);
   bool ony = (x == 0.0);
   bool onz = (onx && ony);
+  double rhos = 0.0;
   double rho = 0.0;
   if (!onz) {
     if (!onx) {
       if (!ony) {
-	rho = sqrt(x * x + y * y);
+	rhos = x * x + y * y;
+	rho = sqrt(rhos);
 	cph = x / rho;
 	sph = y / rho;
 	// goes to 25
       } else { // label 20
+	rhos = y * y;
 	rho = (y > 0.0) ? y : -y;
 	cph = 0.0;
 	sph = (y > 0.0) ? 1.0 : -1.0;
 	// goes to 25
       }
     } else { // label 15
+      rhos = x * x;
       rho = (x > 0.0) ? x : -x;
       cph = (x > 0.0) ? 1.0 : -1.0;
       sph = 0.0;
@@ -1513,7 +1517,7 @@ void polar(
     }
   } else { // label 35
     if (!onz) {
-      r = sqrt(rho * rho + z * z);
+      r = sqrt(rhos + z * z);
       cth = z / r;
       sth = rho / r;
       // returns