Copy of aRts for Trinity modifications

git-svn-id: svn://anonsvn.kde.org/home/kde/branches/trinity/dependencies/arts@1070145 283d02a7-25f6-0310-bc7c-ecb5cbfe19da

1 files changed, 334 insertions, 0 deletions

diff --git a/flow/gsl/gslmathtest.c b/flow/gsl/gslmathtest.c
new file mode 100644
index 0000000..55130aa
--- /dev/null
+++ b/flow/gsl/gslmathtest.c
@@ -0,0 +1,334 @@
+/* GSL - Generic Sound Layer
+ * Copyright (C) 2001 Stefan Westerfeld and Tim Janik
+ *
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public
+ * License as published by the Free Software Foundation; either
+ * version 2 of the License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General
+ * Public License along with this program; if not, write to the
+ * Free Software Foundation, Inc., 59 Temple Place, Suite 330,
+ * Boston, MA 02111-1307, USA.
+ */
+#include <gsl/gslmath.h>
+#include <gsl/gslfilter.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+
+#define	PREC	"15"
+
+static void
+usage (char *s)
+{
+  printf ("usage: gslmathtest %s\n", s);
+  exit (1);
+}
+
+int
+main (int   argc,
+      char *argv[])
+{
+  gchar *arg;
+
+  if (argc < 2)
+    goto abort;
+
+  if (strcmp (argv[1], "rf") == 0)
+    {
+      double x, y, z;
+      if (argc != 5)
+	usage ("rf <x> <y> <z>");
+      x = atof (argv[2]);
+      y = atof (argv[3]);
+      z = atof (argv[4]);
+      
+      printf ("rf(%f, %f, %f) = %."PREC"f\n", x, y, z, gsl_ellip_rf (x, y, z));
+    }
+  else if (strcmp (argv[1], "F") == 0)
+    {
+      double phi, ak;
+      if (argc != 4)
+	usage ("F <phi> <ak>");
+      phi = atof (argv[2]);
+      ak = atof (argv[3]);
+      
+      printf ("F(%f, %f) = %."PREC"f\n", phi, ak, gsl_ellip_F (phi, ak));
+    }
+  else if (strcmp (argv[1], "sn") == 0)
+    {
+      double u, emmc;
+      if (argc != 4)
+	usage ("sn <u> <emmc>");
+      u = atof (argv[2]);
+      emmc = atof (argv[3]);
+      
+      printf ("sn(%f, %f) = %."PREC"f\n", u, emmc, gsl_ellip_sn (u, emmc));
+    }
+  else if (strcmp (argv[1], "snc") == 0)
+    {
+      GslComplex u, emmc;
+      if (argc != 6)
+	usage ("sn <u.re> <u.im> <emmc.re> <emmc.im>");
+      u.re = atof (argv[2]);
+      u.im = atof (argv[3]);
+      emmc.re = atof (argv[4]);
+      emmc.im = atof (argv[5]);
+      
+      printf ("snc(%s, %s) = %s\n",
+	      gsl_complex_str (u),
+	      gsl_complex_str (emmc),
+	      gsl_complex_str (gsl_complex_ellip_sn (u, emmc)));
+    }
+  else if (strcmp (argv[1], "sci_snc") == 0)
+    {
+      GslComplex u, k2;
+      if (argc != 6)
+	usage ("sci_sn <u.re> <u.im> <k2.re> <k2.im>");
+      u.re = atof (argv[2]);
+      u.im = atof (argv[3]);
+      k2.re = atof (argv[4]);
+      k2.im = atof (argv[5]);
+      
+      printf ("sci_snc(%s, %s) = %s\n",
+	      gsl_complex_str (u),
+	      gsl_complex_str (k2),
+	      gsl_complex_str (gsl_complex_ellip_sn (u, gsl_complex_sub (gsl_complex (1.0, 0), k2))));
+    }
+  else if (strcmp (argv[1], "asn") == 0)
+    {
+      double y, emmc;
+      if (argc != 4)
+	usage ("asn <y> <emmc>");
+      y = atof (argv[2]);
+      emmc = atof (argv[3]);
+      
+      printf ("asn(%f, %f) = %."PREC"f\n", y, emmc, gsl_ellip_asn (y, emmc));
+    }
+  else if (strcmp (argv[1], "asnc") == 0)
+    {
+      GslComplex y, emmc;
+      if (argc != 6)
+	usage ("asnc <y.re> <y.im> <emmc.re> <emmc.im>");
+      y.re = atof (argv[2]);
+      y.im = atof (argv[3]);
+      emmc.re = atof (argv[4]);
+      emmc.im = atof (argv[5]);
+      
+      printf ("asnc(%s, %s) = %s\n",
+	      gsl_complex_str (y), gsl_complex_str (emmc),
+	      gsl_complex_str (gsl_complex_ellip_asn (y, emmc)));
+      printf ("asn(%f, %f = %."PREC"f\n",
+	      y.re, emmc.re, gsl_ellip_asn (y.re, emmc.re));
+    }
+  else if (strcmp (argv[1], "sci_sn") == 0)
+    {
+      double u, k2;
+      if (argc != 4)
+	usage ("sci_sn <u> <k2>");
+      u = atof (argv[2]);
+      k2 = atof (argv[3]);
+      
+      printf ("sci_sn(%f, %f) = %."PREC"f\n", u, k2, gsl_ellip_sn (u, 1.0 - k2));
+    }
+  else if (strcmp (argv[1], "sci_asn") == 0)
+    {
+      double y, k2;
+      if (argc != 4)
+	usage ("sci_asn <y> <k2>");
+      y = atof (argv[2]);
+      k2 = atof (argv[3]);
+      
+      printf ("sci_asn(%f, %f) = %."PREC"f\n", y, k2, gsl_ellip_asn (y, 1.0 - k2));
+    }
+  else if (strcmp (argv[1], "sci_asnc") == 0)
+    {
+      GslComplex y, k2;
+      if (argc != 6)
+	usage ("sci_asnc <y.re> <y.im> <k2.re> <k2.im>");
+      y.re = atof (argv[2]);
+      y.im = atof (argv[3]);
+      k2.re = atof (argv[4]);
+      k2.im = atof (argv[5]);
+      
+      printf ("sci_asnc(%s, %s) = %s\n",
+	      gsl_complex_str (y), gsl_complex_str (k2),
+	      gsl_complex_str (gsl_complex_ellip_asn (y, gsl_complex_sub (gsl_complex (1.0, 0), k2))));
+      printf ("asn(%f, %f = %."PREC"f\n",
+	      y.re, k2.re, gsl_ellip_asn (y.re, 1.0 - k2.re));
+    }
+  else if (strcmp (argv[1], "sin") == 0)
+    {
+      GslComplex phi;
+      if (argc != 4)
+	usage ("sin <phi.re> <phi.im>");
+      phi.re = atof (argv[2]);
+      phi.im = atof (argv[3]);
+      
+      printf ("sin(%s) = %s\n",
+	      gsl_complex_str (phi),
+	      gsl_complex_str (gsl_complex_sin (phi)));
+    }
+  else if (strcmp (argv[1], "cos") == 0)
+    {
+      GslComplex phi;
+      if (argc != 4)
+	usage ("cos <phi.re> <phi.im>");
+      phi.re = atof (argv[2]);
+      phi.im = atof (argv[3]);
+      
+      printf ("cos(%s) = %s\n",
+	      gsl_complex_str (phi),
+	      gsl_complex_str (gsl_complex_cos (phi)));
+    }
+  else if (strcmp (argv[1], "tan") == 0)
+    {
+      GslComplex phi;
+      if (argc != 4)
+	usage ("tan <phi.re> <phi.im>");
+      phi.re = atof (argv[2]);
+      phi.im = atof (argv[3]);
+      
+      printf ("tan(%s) = %s\n",
+	      gsl_complex_str (phi),
+	      gsl_complex_str (gsl_complex_tan (phi)));
+    }
+  else if (strcmp (argv[1], "sinh") == 0)
+    {
+      GslComplex phi;
+      if (argc != 4)
+	usage ("sinh <phi.re> <phi.im>");
+      phi.re = atof (argv[2]);
+      phi.im = atof (argv[3]);
+      
+      printf ("sinh(%s) = %s\n",
+	      gsl_complex_str (phi),
+	      gsl_complex_str (gsl_complex_sinh (phi)));
+    }
+  else if (strcmp (argv[1], "cosh") == 0)
+    {
+      GslComplex phi;
+      if (argc != 4)
+	usage ("cosh <phi.re> <phi.im>");
+      phi.re = atof (argv[2]);
+      phi.im = atof (argv[3]);
+      
+      printf ("cosh(%s) = %s\n",
+	      gsl_complex_str (phi),
+	      gsl_complex_str (gsl_complex_cosh (phi)));
+    }
+  else if (strcmp (argv[1], "tanh") == 0)
+    {
+      GslComplex phi;
+      if (argc != 4)
+	usage ("tanh <phi.re> <phi.im>");
+      phi.re = atof (argv[2]);
+      phi.im = atof (argv[3]);
+      
+      printf ("tanh(%s) = %s\n",
+	      gsl_complex_str (phi),
+	      gsl_complex_str (gsl_complex_tanh (phi)));
+    }
+  else if (strcmp (argv[1], "t1") == 0)
+    {
+      guint order;
+      double f, e;
+      if (argc != 5)
+	usage ("t1 <order> <freq> <epsilon>");
+      order = atoi (argv[2]);
+      f = atof (argv[3]);
+      e = atof (argv[4]);
+      f *= GSL_PI / 2.;
+      e = gsl_trans_zepsilon2ss (e);
+      {
+	double a[order + 1], b[order + 1];
+	gsl_filter_tscheb1 (order, f, e, a, b);
+	g_print ("# Tschebyscheff Type1 order=%u freq=%f s^2epsilon=%f norm0=%f:\n",
+		 order, f, e,
+		 gsl_poly_eval (order, a, 1) / gsl_poly_eval (order, b, 1));
+	g_print ("H%u(z)=%s/%s\n", order,
+		 gsl_poly_str (order, a, "z"),
+		 gsl_poly_str (order, b, "z"));
+      }
+    }
+  else if (strcmp (argv[1], "t2") == 0)
+    {
+      guint order;
+      double fc, fr, e;
+      if (argc != 6)
+	usage ("t1 <order> <freqc> <freqr> <epsilon>");
+      order = atoi (argv[2]);
+      fc = atof (argv[3]);
+      fr = atof (argv[4]);
+      e = atof (argv[5]);
+      fc *= GSL_PI / 2.;
+      fr *= GSL_PI / 2.;
+      e = gsl_trans_zepsilon2ss (e);
+      {
+	double a[order + 1], b[order + 1];
+	gsl_filter_tscheb2 (order, fc, fr, e, a, b);
+	g_print ("# Tschebyscheff Type2 order=%u freq_c=%f freq_r=%f s^2epsilon=%f norm=%f:\n",
+		 order, fc, fr, e,
+		 gsl_poly_eval (order, a, 1) / gsl_poly_eval (order, b, 1));
+	g_print ("H%u(z)=%s/%s\n", order,
+		 gsl_poly_str (order, a, "z"),
+		 gsl_poly_str (order, b, "z"));
+      }
+    }
+  else if (strncmp (argv[1], "test", 4) == 0)
+    {
+      guint order;
+      arg = argv[1] + 4;
+      if (argc != argc)
+	usage ("test");
+      order = 2;
+      {
+	double a[100] = { 1, 2, 1 }, b[100] = { 1, -3./2., 0.5 };
+	g_print ("# Test order=%u  norm=%f:\n",
+		 order,
+		 gsl_poly_eval (order, a, 1) / gsl_poly_eval (order, b, 1));
+	g_print ("H%u(z)=%s/%s\n", order,
+		 gsl_poly_str (order, a, "z"),
+		 gsl_poly_str (order, b, "z"));
+	if (*arg)
+	  {
+	    GslComplex root, roots[100];
+	    guint i;
+
+	    if (*arg == 'r')
+	      {
+		g_print ("#roots:\n");
+		gsl_poly_complex_roots (order, a, roots);
+		for (i = 0; i < order; i++)
+		  {
+		    root = gsl_complex_div (gsl_complex (1, 0), roots[i]);
+		    g_print ("%+.14f %+.14f # %.14f\n", root.re, root.im, gsl_complex_abs (root));
+		  }
+	      }
+	    if (*arg == 'p')
+	      {
+		g_print ("#poles:\n");
+		gsl_poly_complex_roots (order, b, roots);
+		for (i = 0; i < order; i++)
+		  {
+		    root = gsl_complex_div (gsl_complex (1, 0), roots[i]);
+		    g_print ("%+.14f %+.14f # %.14f\n", root.re, root.im, gsl_complex_abs (root));
+		  }
+	      }
+	  }
+      }
+    }
+  else
+    {
+    abort:
+      usage ("{rf|F|sn|snc|sci_sn|sci_snc|asn|asnc|sci_asn|sci_asnc|sin(h)|cos(h)|tan(h)|t1|t2} ...");
+    }
+      
+  return 0;
+}