use new getArgument(re, im) dealing with numerical precision

Author: bourgesl

src/main/java/fr/jmmc/jmcs/util/NumberUtils.java

@@ -48,7 +48,12 @@
 public final class NumberUtils {
 
     /** Smallest positive number used in double comparisons (rounding). */
-    public final static double EPSILON = 1e-6d;
+    public final static double EPSILON = 1e-6;
+    /* PI and PI/2 constants */
+    public final static double PI = Math.PI;
+    public final static double PI_HALF = PI / 2.0;
+    /** Typical scale precision in re/im comparisons (rounding). */
+    public static final double ARG_EPSILON = Math.ulp(1.0);
     /* shared Double instances */
     /** shared Double = NaN instance */
     public final static Double DBL_NAN = Double.valueOf(Double.NaN);
@@ -61,6 +66,33 @@ public final class NumberUtils {
     /** scientific formatter */
     private final static NumberFormat _fmtScience = new DecimalFormat("0.0##E0");
 
+    public static double getArgumentInDegrees(final double re, final double im) {
+        return Math.toDegrees(getArgument(re, im));
+    }
+
+    public static double getArgument(final double re, final double im) {
+        // check |re| == 0
+        final double normRe = Math.abs(re);
+        if (normRe == 0.0) {
+            return (im >= 0.0) ? PI_HALF : -PI_HALF;
+        }
+        // check |im| == 0
+        final double normIm = Math.abs(im);
+        if (normIm == 0.0) {
+            return (re >= 0.0) ? 0.0 : PI;
+        }
+        final double epsilon = ARG_EPSILON * (normRe + normIm);
+        // check |re| << |im| or |im| << |re|
+        if (normIm <= epsilon) {
+            return (re >= 0.0) ? 0.0 : PI;
+        }
+        // check |im| << |re|
+        if (normRe <= epsilon) {
+            return (im >= 0.0) ? PI_HALF : -PI_HALF;
+        }
+        return Math.atan2(im, re);
+    }
+
     /**
      * Private constructor
      */

src/main/java/fr/jmmc/oitools/model/OIT3.java

@@ -19,6 +19,7 @@
  ******************************************************************************/
 package fr.jmmc.oitools.model;
 
+import fr.jmmc.jmcs.util.NumberUtils;
 import fr.jmmc.oitools.OIFitsConstants;
 import fr.jmmc.oitools.meta.ArrayColumnMeta;
 import static fr.jmmc.oitools.meta.CellMeta.NO_STR_VALUES;
@@ -437,13 +438,13 @@ public double[] getPosAngle() {
 
                 switch (j) {
                     case 1: // AB
-                        angle[i] = Math.atan2(u1coord[i], v1coord[i]);
+                        angle[i] = NumberUtils.getArgument(v1coord[i], u1coord[i]);
                         break;
                     case 2: // BC
-                        angle[i] = Math.atan2(u2coord[i], v2coord[i]);
+                        angle[i] = NumberUtils.getArgument(v2coord[i], u2coord[i]);
                         break;
                     case 3: // AC
-                        angle[i] = Math.atan2(u3, v3);
+                        angle[i] = NumberUtils.getArgument(v3, u3);
                         break;
                     default:
                         angle[i] = 0d;

src/main/java/fr/jmmc/oitools/model/OIVis.java

@@ -19,6 +19,7 @@
  ******************************************************************************/
 package fr.jmmc.oitools.model;
 
+import fr.jmmc.jmcs.util.NumberUtils;
 import fr.jmmc.oitools.OIFitsConstants;
 import fr.jmmc.oitools.meta.ArrayColumnMeta;
 import static fr.jmmc.oitools.meta.CellMeta.NO_STR_VALUES;
@@ -578,7 +579,7 @@ public double[] getPosAngle() {
             final double[] vcoord = getVCoord();
 
             for (int i = 0, j; i < nRows; i++) {
-                angle[i] = Math.toDegrees(Math.atan2(ucoord[i], vcoord[i]));
+                angle[i] = NumberUtils.getArgumentInDegrees(vcoord[i], ucoord[i]);
             }
 
             this.setColumnDerivedValue(OIFitsConstants.COLUMN_POS_ANGLE, angle);

src/main/java/fr/jmmc/oitools/model/OIVis2.java

@@ -19,6 +19,7 @@
  ******************************************************************************/
 package fr.jmmc.oitools.model;
 
+import fr.jmmc.jmcs.util.NumberUtils;
 import fr.jmmc.oitools.OIFitsConstants;
 import fr.jmmc.oitools.meta.ArrayColumnMeta;
 import static fr.jmmc.oitools.meta.CellMeta.NO_STR_VALUES;
@@ -391,7 +392,7 @@ public double[] getPosAngle() {
             final double[] vcoord = getVCoord();
 
             for (int i = 0, j; i < nRows; i++) {
-                angle[i] = Math.toDegrees(Math.atan2(ucoord[i], vcoord[i]));
+                angle[i] = NumberUtils.getArgumentInDegrees(vcoord[i], ucoord[i]);
             }
 
             this.setColumnDerivedValue(OIFitsConstants.COLUMN_POS_ANGLE, angle);

src/main/java/fr/jmmc/oitools/processing/BeamEstimator.java

@@ -3,6 +3,7 @@
  ******************************************************************************/
 package fr.jmmc.oitools.processing;
 
+import fr.jmmc.jmcs.util.NumberUtils;
 import fr.jmmc.oitools.model.IndexMask;
 import fr.jmmc.oitools.model.OIData;
 import fr.jmmc.oitools.model.OIFitsFile;
@@ -276,7 +277,7 @@ private static BeamInfo computeBeam(final double[] covMatrix) {
             final double rx = Math.toDegrees(SCALE_FACTOR / Math.sqrt(results.ev1)) * DEG_IN_MILLI_ARCSEC;
             final double ry = Math.toDegrees(SCALE_FACTOR / Math.sqrt(results.ev2)) * DEG_IN_MILLI_ARCSEC;
 
-            double angle = Math.toDegrees(Math.atan2(results.vec1[0], results.vec1[1]));
+            double angle = NumberUtils.getArgumentInDegrees(results.vec1[1], results.vec1[0]);
 
             // correct angle in range [-90; 90] as collinearity !
             if (angle < -90.0) {