From b5f823f599894139345fb71eb8040f7058548adf Mon Sep 17 00:00:00 2001
From: rlar <rlar>
Date: Sat, 6 Aug 2011 17:20:58 +0000
Subject: [PATCH] cph(vec), rewrite, cleanup, and add missing test_cx_cph.c

---
 ChangeLog                      |  6 ++++
 src/maths/cmaths/Makefile.am   |  2 +-
 src/maths/cmaths/cmath1.c      | 18 ++++------
 src/maths/cmaths/test_cx_cph.c | 61 ++++++++++++++++++++++++++++++++++
 4 files changed, 74 insertions(+), 13 deletions(-)
 create mode 100644 src/maths/cmaths/test_cx_cph.c

diff --git a/ChangeLog b/ChangeLog
index 397f092fe..fce21fc8f 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,9 @@
+2011-08-06  Robert Larice
+ * src/maths/cmaths/Makefile.am ,
+ * src/maths/cmaths/cmath1.c ,
+ * src/maths/cmaths/test_cx_cph.c :
+   cph(vec), rewrite, cleanup, and add missing test_cx_cph.c
+
 2011-08-06  Holger Vogt
  * fteext.h, cmaths/makefile.am, cmath1.c, cmath1.h, parse.c:
    function cph(vec) delivers phase from a vector without jumps at +-PI
diff --git a/src/maths/cmaths/Makefile.am b/src/maths/cmaths/Makefile.am
index 0ad441c14..fbf2877f1 100644
--- a/src/maths/cmaths/Makefile.am
+++ b/src/maths/cmaths/Makefile.am
@@ -28,7 +28,7 @@ test_cx_ph_LDADD = \
 		$(TCL_LIB_SPEC)
 
 test_cx_cph_SOURCES = \
-               test_cx_ph.c
+               test_cx_cph.c
 
 test_cx_cph_LDADD = \
                 libcmaths.la \
diff --git a/src/maths/cmaths/cmath1.c b/src/maths/cmaths/cmath1.c
index c86b7c224..d0f26f265 100644
--- a/src/maths/cmaths/cmath1.c
+++ b/src/maths/cmaths/cmath1.c
@@ -76,23 +76,17 @@ cx_cph(void *data, short int type, int length, int *newlength, short int *newtyp
 {
     double *d = alloc_d(length);
     ngcomplex_t *cc = (ngcomplex_t *) data;
-    int i, j, jmin;
-    int period = 0;
-    double ph[3], phd[3];
+    int i;
 
     *newlength = length;
     *newtype = VF_REAL;
     if (type == VF_COMPLEX) {
-        d[0] = radtodeg(cph(&cc[0]));
+        double last_ph = cph(&cc[0]);
+        d[0] = radtodeg(last_ph);
         for (i = 1; i < length; i++) {
-            jmin = 0;
-            for (j = 0; j < 3; j++) {
-                ph[j] = radtodeg((period+j-1) * 2*M_PI + cph(&cc[i]));
-                phd[j] = abs(ph[j] - d[i-1]);
-                jmin = (phd[j] < phd[jmin]) ? j : jmin;
-            }
-            d[i] = ph[jmin];
-            period += jmin-1;
+            double ph = cph(&cc[i]);
+            last_ph = ph - (2*M_PI) * floor((ph - last_ph)/(2*M_PI) + 0.5);
+            d[i] = radtodeg(last_ph);
         }
     }
     /* Otherwise it is 0, but tmalloc zeros the stuff already. */
diff --git a/src/maths/cmaths/test_cx_cph.c b/src/maths/cmaths/test_cx_cph.c
new file mode 100644
index 000000000..61258d0ed
--- /dev/null
+++ b/src/maths/cmaths/test_cx_cph.c
@@ -0,0 +1,61 @@
+#include <stdio.h>
+#include <math.h>
+#include <float.h>
+
+#include <config.h>
+#include <memory.h>
+#include <dvec.h>
+#include <complex.h>
+
+#include "defines.h"
+#include "cmath.h"
+#include "cmath1.h"
+
+FILE *cp_err;
+
+static int eq_p(double a, double b, double eps)
+{
+    return fabs(a-b) < eps;
+}
+
+
+int
+main(void)
+{
+    ngcomplex_t *c = NULL;
+    double *d = NULL;
+    short int t1;
+    short int t2;
+    int n1, n2;
+    double eps = DBL_EPSILON;
+
+    cp_err = stderr;
+    n1 = 9;
+    t1 = VF_COMPLEX;
+    c = alloc_c(n1);
+
+    realpart(&c[0]) =  0.0; imagpart(&c[0]) = +1.0;  /* i^1 */
+    realpart(&c[1]) = -1.0; imagpart(&c[1]) =  0.0;  /* i^2 */
+    realpart(&c[2]) =  0.0; imagpart(&c[2]) = -1.0;  /* i^3 */
+    realpart(&c[3]) = +1.0; imagpart(&c[3]) =  0.0;  /* i^4 */
+    realpart(&c[4]) =  0.0; imagpart(&c[4]) = +1.0;  /* i^5 */
+    realpart(&c[5]) = +1.0; imagpart(&c[5]) =  0.0;  /* i^4 */
+    realpart(&c[6]) =  0.0; imagpart(&c[6]) = -1.0;  /* i^3 */
+    realpart(&c[7]) = -1.0; imagpart(&c[7]) =  0.0;  /* i^2 */
+    realpart(&c[8]) =  0.0; imagpart(&c[8]) = +1.0;  /* i^1 */
+
+    d = (double *) cx_cph((void *) c, t1, n1, &n2, &t2);
+
+    if ( eq_p(1*M_PI/2, d[0], eps) &&
+         eq_p(2*M_PI/2, d[1], eps) &&
+         eq_p(3*M_PI/2, d[2], eps) &&
+         eq_p(4*M_PI/2, d[3], eps) &&
+         eq_p(5*M_PI/2, d[4], eps) &&
+         eq_p(4*M_PI/2, d[5], eps) &&
+         eq_p(3*M_PI/2, d[6], eps) &&
+         eq_p(2*M_PI/2, d[7], eps) &&
+         eq_p(1*M_PI/2, d[8], eps) )
+	return 0;
+    else
+	return 1;
+}