Implemented julian(x, a, b)

michele-summo · Apr 15, 2022 · 3af3414 · 3af3414
1 parent 8952418
commit 3af3414
Show file tree

Hide file tree

Showing 2 changed files with 20 additions and 1 deletion.
diff --git a/src/include/sffe_cmplx_gsl.h b/src/include/sffe_cmplx_gsl.h
@@ -13,7 +13,7 @@
 #include <gsl/gsl_complex.h>
 #include <gsl/gsl_complex_math.h>
 
-#define sffnctscount 78
+#define sffnctscount 79
 #define sfvarscount 6
 #define cmplxset(c, r, i) GSL_SET_COMPLEX(&c, r, i)
 #define real(c) GSL_REAL((c))
@@ -102,6 +102,8 @@ sfarg *sfcoti(sfarg *const p);
 sfarg *sftrunc(sfarg *const p);
 sfarg *sfsawtooth(sfarg *const p);
 sfarg *sftwave(sfarg *const p);
+/* other functions */
+sfarg *sfjulian(sfarg *const p); // |z^a|*e^(i*b*arg(z))
 
 /*const eval*/
 void sfcPI(sfNumber *cnst);

diff --git a/src/sffe/sffe_cmplx_gsl.cpp b/src/sffe/sffe_cmplx_gsl.cpp
@@ -102,6 +102,8 @@ const sffunction sfcmplxfunc[sffnctscount] = {
     {sfsawtooth, 1, "sawtooth\0"},
     {sftwave, 1, "twave\0"},
 
+    {sfjulian, 3, "julian\0"},
+
     {NULL, 1, "rad\0"},
     {NULL, 1, "deg\0"},
     {NULL, 1, "sign\0"},
@@ -676,6 +678,21 @@ sfarg *sftwave(sfarg *const p)
     return sfaram1(p);
 }
 
+sfarg *sfjulian(sfarg *const p)
+{
+    gsl_complex z = sfvalue(sfaram3(p));
+    double mag = gsl_complex_abs(
+                gsl_complex_pow(z, sfvalue(sfaram2(p))));
+    gsl_complex b = sfvalue(sfaram1(p));
+    double arg = gsl_complex_arg(z);
+    double byg = exp(-GSL_IMAG(b)*arg);
+    double bxg = arg * GSL_REAL(b);
+
+    GSL_REAL(sfvalue(p)) = mag * byg * cos(bxg);
+    GSL_IMAG(sfvalue(p)) = mag * byg * sin(bxg);
+    return sfaram3(p);
+}
+
 // const eval
 void sfcPI(sfNumber *cnst) { GSL_SET_COMPLEX(cnst, 4 * atan(1), 0); }