Skip to content
Browse files

Merge branch 'bessel_array'

  • Loading branch information...
2 parents 95185d9 + 78296d7 commit 7a7ebf8e4865da694dd0bfb99e9527cadcb4b553 @leto committed Mar 25, 2016
Showing with 280 additions and 25 deletions.
  1. +11 −0 Changes
  2. +64 −21 pod/SF.pod
  3. +126 −0 swig/SF.i
  4. +26 −3 swig/gsl_typemaps.i
  5. +53 −1 t/SF.t
View
11 Changes
@@ -4,6 +4,17 @@
- Fixed error messages on Matrix set_col and set_row
methods (Jon Schutz)
+ - Implemented array-based Bessel functions:
+ gsl_sf_bessel_Jn_array
+ gsl_sf_bessel_Yn_array
+ gsl_sf_bessel_In_array
+ gsl_sf_bessel_In_scaled_array
+ gsl_sf_bessel_Kn_array
+ gsl_sf_bessel_Kn_scaled_array
+ gsl_sf_bessel_jl_array
+ gsl_sf_bessel_jl_steed_array
+ gsl_sf_bessel_yl_array
+ gsl_sf_bessel_il_scaled_array
=head1 v0.35 - April 2015
View
85 pod/SF.pod
@@ -766,11 +766,11 @@ These routines compute the regular cylindrical Bessel function of order n, J_n($
=over
-=item C<gsl_sf_bessel_Jn_array($nmin, $nmax, $x, $result_array)>
+=item C<gsl_sf_bessel_Jn_array($nmin, $nmax, $x)>
This routine computes the values of the regular cylindrical Bessel functions
-J_n($x) for n from $nmin to $nmax inclusive, storing the results in the array
-$result_array. The values are computed using recurrence relations for
+J_n($x) for n from $nmin to $nmax inclusive, returning the results in an array
+reference. The values are computed using recurrence relations for
efficiency, and therefore may differ slightly from the exact values.
=back
@@ -807,8 +807,13 @@ These routines compute the irregular cylindrical Bessel function of order $n, Y_
=over
-=item C<gsl_sf_bessel_Yn_array>
+=item C<gsl_sf_bessel_Yn_array($nmin, $nmax, $x)>
+This routine computes the values of the irregular cylindrical Bessel functions
+Y_n(x) for n from $nmin to $nmax inclusive, returning the results in an array
+reference. The domain of the function is $x>0. The values are computed using
+recurrence relations for efficiency, and therefore may differ slightly from the
+exact values.
=back
@@ -829,7 +834,7 @@ These routines compute the regular modified cylindrical Bessel function of zerot
=item C<gsl_sf_bessel_I1($x)>
-These routines compute the regular modified cylindrical Bessel function of first order, I_1(x).
+
=back
@@ -845,9 +850,13 @@ These routines compute the regular modified cylindrical Bessel function of order
=over
-=item C<gsl_sf_bessel_In_array>
-
+=item C<gsl_sf_bessel_In_array($nmin, $nmax, $x)>
+This routine computes the values of the regular modified cylindrical Bessel
+functions I_n(x) for n from $nmin to $nmax inclusive, returning the results in an
+array reference. The start of the range nmin must be positive or zero. The
+values are computed using recurrence relations for efficiency, and therefore may
+differ slightly from the exact values.
=back
@@ -883,9 +892,13 @@ These routines compute the scaled regular modified cylindrical Bessel function o
=over
-=item C<gsl_sf_bessel_In_scaled_array>
-
+=item C<gsl_sf_bessel_In_scaled_array($nmin, $nmax, $x)>
+This routine computes the values of the scaled regular cylindrical Bessel
+functions exp(-|$x|) I_n($x) for n from $nmin to $nmax inclusive, returning the
+results in an array reference. The start of the range nmin must be positive
+or zero. The values are computed using recurrence relations for efficiency, and
+therefore may differ slightly from the exact values.
=back
@@ -921,9 +934,13 @@ These routines compute the irregular modified cylindrical Bessel function of ord
=over
-=item C<gsl_sf_bessel_Kn_array>
-
+=item C<gsl_sf_bessel_Kn_array($nmin, $nmax, $x)>
+This routine computes the values of the Irregular Modified Cylindrical Bessel Functions
+K_n($x) for n from $nmin to $nmax inclusive, returning the results in an array
+reference. The values are computed using recurrence relations for
+efficiency, and therefore may differ slightly from the exact values.
+This function is only defined for positive $x and with throw an exception otherwise.
=back
@@ -959,8 +976,14 @@ These routines compute the scaled irregular modified cylindrical Bessel function
=over
-=item C<gsl_sf_bessel_Kn_scaled_array >
+=item C<gsl_sf_bessel_Kn_scaled_array($nmin, $max, $x) >
+This routine computes the values of the scaled irregular cylindrical Bessel
+functions exp(x) K_n(x) for n from nmin to nmax inclusive, storing the results
+in the array result_array. The start of the range nmin must be positive or zero.
+The domain of the function is x>0. The values are computed using recurrence
+relations for efficiency, and therefore may differ slightly from the exact
+values.
=back
@@ -1007,17 +1030,25 @@ These routines compute the scaled irregular modified cylindrical Bessel function
=over
-=item C<gsl_sf_bessel_jl_array>
-
+=item C<gsl_sf_bessel_jl_array($lmax, $x)>
+This routine computes the values of the regular spherical Bessel functions
+j_l(x) for l from 0 to C<$lmax> inclusive for C<$lmax> >= 0 and $x >= 0, returning
+the results in an array reference. The values are computed using recurrence
+relations for efficiency, and therefore may differ slightly from the exact
+values.
=back
=over
-=item C<gsl_sf_bessel_jl_steed_array>
-
+=item C<gsl_sf_bessel_jl_steed_array($lmax, $x)>
+This routine uses Steed’s method to compute the values of the regular spherical
+Bessel functions j_l(x) for l from 0 to $lmax inclusive for $lmax >= 0 and $x >= 0,
+storing the results in the array result_array. The Steed/Barnett algorithm is
+described in Comp. Phys. Comm. 21, 297 (1981). Steed’s method is more stable
+than the recurrence used in the other functions but is also slower.
=back
@@ -1063,9 +1094,12 @@ These routines compute the scaled irregular modified cylindrical Bessel function
=over
-=item C<gsl_sf_bessel_yl_array>
-
+=item C<gsl_sf_bessel_yl_array($lmax, $x)>
+This routine computes the values of the irregular spherical Bessel functions
+y_l(x) for l from 0 to $lmax inclusive for lmax >= 0, returning the results in an
+array reference. The values are computed using recurrence relations for
+efficiency, and therefore may differ slightly from the exact values.
=back
@@ -1111,9 +1145,13 @@ These routines compute the scaled irregular modified cylindrical Bessel function
=over
-=item C<gsl_sf_bessel_il_scaled_array>
-
+=item C<gsl_sf_bessel_il_scaled_array($lmax, $x)>
+This routine computes the values of the scaled regular modified spherical Bessel
+functions exp(-|x|) i_l(x) for l from 0 to $lmax inclusive for $lmax >= 0,
+returning the results in an array reference. The values are computed using
+recurrence relations for efficiency, and therefore may differ slightly from the
+exact values.
=back
@@ -1159,8 +1197,13 @@ These routines compute the scaled irregular modified cylindrical Bessel function
=over
-=item C<gsl_sf_bessel_kl_scaled_array>
+=item C<gsl_sf_bessel_kl_scaled_array($lmax, $x)>
+This routine computes the values of the scaled irregular modified spherical Bessel
+functions exp($x) k_l($x) for l from 0 to lmax inclusive for $lmax >= 0 and $x>0,
+returning the results in an array reference. The values are computed using
+recurrence relations for efficiency, and therefore may differ slightly from the
+exact values.
=back
View
126 swig/SF.i
@@ -5,6 +5,41 @@
%apply double *OUTPUT { double * sn, double * cn, double * dn, double * sgn };
+// rename wrappers to original
+
+%ignore gsl_sf_bessel_Jn_array;
+%rename (gsl_sf_bessel_Jn_array) gsl_sf_bessel_Jn_array_wrapper;
+array_wrapper* gsl_sf_bessel_Jn_array_wrapper(int nmin, int nmax, double x);
+%ignore gsl_sf_bessel_Yn_array;
+%rename (gsl_sf_bessel_Yn_array) gsl_sf_bessel_Yn_array_wrapper;
+array_wrapper* gsl_sf_bessel_Yn_array_wrapper(int nmin, int nmax, double x);
+%ignore gsl_sf_bessel_In_array;
+%rename (gsl_sf_bessel_In_array) gsl_sf_bessel_In_array_wrapper;
+array_wrapper* gsl_sf_bessel_In_array_wrapper(int nmin, int nmax, double x);
+%ignore gsl_sf_bessel_In_scaled_array;
+%rename (gsl_sf_bessel_In_scaled_array) gsl_sf_bessel_In_scaled_array_wrapper;
+array_wrapper* gsl_sf_bessel_In_scaled_array_wrapper(int nmin, int nmax, double x);
+%ignore gsl_sf_bessel_Kn_array;
+%rename (gsl_sf_bessel_Kn_array) gsl_sf_bessel_Kn_array_wrapper;
+array_wrapper* gsl_sf_bessel_Kn_array_wrapper(int nmin, int nmax, double x);
+%ignore gsl_sf_bessel_Kn_scaled_array;
+%rename (gsl_sf_bessel_Kn_scaled_array) gsl_sf_bessel_Kn_scaled_array_wrapper;
+array_wrapper* gsl_sf_bessel_Kn_scaled_array_wrapper(int nmin, int nmax, double x);
+%ignore gsl_sf_bessel_jl_array;
+%rename (gsl_sf_bessel_jl_array) gsl_sf_bessel_jl_array_wrapper;
+array_wrapper* gsl_sf_bessel_jl_array_wrapper(int lmax, double x);
+%ignore gsl_sf_bessel_jl_steed_array;
+%rename (gsl_sf_bessel_jl_steed_array) gsl_sf_bessel_jl_steed_array_wrapper;
+array_wrapper* gsl_sf_bessel_jl_steed_array_wrapper(int lmax, double x);
+%ignore gsl_sf_bessel_yl_array;
+%rename (gsl_sf_bessel_yl_array) gsl_sf_bessel_yl_array_wrapper;
+array_wrapper* gsl_sf_bessel_yl_array_wrapper(int lmax, double x);
+%ignore gsl_sf_bessel_il_scaled_array;
+%rename (gsl_sf_bessel_il_scaled_array) gsl_sf_bessel_il_scaled_array_wrapper;
+array_wrapper* gsl_sf_bessel_il_scaled_array_wrapper(int lmax, double x);
+%ignore gsl_sf_bessel_kl_scaled_array;
+%rename (gsl_sf_bessel_kl_scaled_array) gsl_sf_bessel_kl_scaled_array_wrapper;
+
%{
#include "gsl/gsl_types.h"
#include "gsl/gsl_version.h"
@@ -40,7 +75,98 @@
#include "gsl/gsl_sf_transport.h"
#include "gsl/gsl_sf_trig.h"
#include "gsl/gsl_sf_zeta.h"
+
+ array_wrapper* gsl_sf_bessel_Jn_array_wrapper(int nmin, int nmax, double x) {
+ int ret;
+ array_wrapper * wrapper = array_wrapper_alloc(nmax - nmin + 1, awDouble);
+ ret = gsl_sf_bessel_Jn_array(nmin,nmax,x, (double*)(wrapper->data));
+ if (ret != GSL_SUCCESS)
+ croak("Math::GSL::SF::gsl_sf_bessel_Jn_array returned %s", gsl_strerror(ret));
+ return wrapper;
+ }
+ array_wrapper* gsl_sf_bessel_Yn_array_wrapper(int nmin, int nmax, double x) {
+ int ret;
+ array_wrapper * wrapper = array_wrapper_alloc(nmax - nmin + 1, awDouble);
+ ret = gsl_sf_bessel_Yn_array(nmin,nmax,x, (double*)(wrapper->data));
+ if (ret != GSL_SUCCESS)
+ croak("Math::GSL::SF::gsl_sf_bessel_Yn_array returned %s", gsl_strerror(ret));
+ return wrapper;
+ }
+ array_wrapper* gsl_sf_bessel_In_array_wrapper(int nmin, int nmax, double x) {
+ int ret;
+ array_wrapper * wrapper = array_wrapper_alloc(nmax - nmin + 1, awDouble);
+ ret = gsl_sf_bessel_In_array(nmin,nmax,x, (double*)(wrapper->data));
+ if (ret != GSL_SUCCESS)
+ croak("Math::GSL::SF::gsl_sf_bessel_In_array returned %s", gsl_strerror(ret));
+ return wrapper;
+ }
+ array_wrapper* gsl_sf_bessel_In_scaled_array_wrapper(int nmin, int nmax, double x) {
+ int ret;
+ array_wrapper * wrapper = array_wrapper_alloc(nmax - nmin + 1, awDouble);
+ ret = gsl_sf_bessel_In_scaled_array(nmin,nmax,x, (double*)(wrapper->data));
+ if (ret != GSL_SUCCESS)
+ croak("Math::GSL::SF::gsl_sf_bessel_In_scaled_array returned %s", gsl_strerror(ret));
+ return wrapper;
+ }
+ array_wrapper* gsl_sf_bessel_Kn_array_wrapper(int nmin, int nmax, double x) {
+ int ret;
+ array_wrapper * wrapper = array_wrapper_alloc(nmax - nmin + 1, awDouble);
+ ret = gsl_sf_bessel_Kn_array(nmin,nmax,x, (double*)(wrapper->data));
+ if (ret != GSL_SUCCESS)
+ croak("Math::GSL::SF::gsl_sf_bessel_Kn_array returned %s", gsl_strerror(ret));
+ return wrapper;
+ }
+ array_wrapper* gsl_sf_bessel_Kn_scaled_array_wrapper(int nmin, int nmax, double x) {
+ int ret;
+ array_wrapper * wrapper = array_wrapper_alloc(nmax - nmin + 1, awDouble);
+ ret = gsl_sf_bessel_Kn_scaled_array(nmin,nmax,x, (double*)(wrapper->data));
+ if (ret != GSL_SUCCESS)
+ croak("Math::GSL::SF::gsl_sf_bessel_Kn_scaled_array returned %s", gsl_strerror(ret));
+ return wrapper;
+ }
+ array_wrapper* gsl_sf_bessel_jl_array_wrapper(int lmax, double x) {
+ int ret;
+ array_wrapper * wrapper = array_wrapper_alloc(lmax + 1, awDouble);
+ ret = gsl_sf_bessel_jl_array(lmax,x, (double*)(wrapper->data));
+ if (ret != GSL_SUCCESS)
+ croak("Math::GSL::SF::gsl_sf_bessel_jl_array returned %s", gsl_strerror(ret));
+ return wrapper;
+ }
+ array_wrapper* gsl_sf_bessel_jl_steed_array_wrapper(int lmax, double x) {
+ int ret;
+ array_wrapper * wrapper = array_wrapper_alloc(lmax + 1, awDouble);
+ ret = gsl_sf_bessel_jl_steed_array(lmax,x, (double*)(wrapper->data));
+ if (ret != GSL_SUCCESS)
+ croak("Math::GSL::SF::gsl_sf_bessel_jl_steed_array returned %s", gsl_strerror(ret));
+ return wrapper;
+ }
+ array_wrapper* gsl_sf_bessel_yl_array_wrapper(int lmax, double x) {
+ int ret;
+ array_wrapper * wrapper = array_wrapper_alloc(lmax + 1, awDouble);
+ ret = gsl_sf_bessel_yl_array(lmax,x, (double*)(wrapper->data));
+ if (ret != GSL_SUCCESS)
+ croak("Math::GSL::SF::gsl_sf_bessel_yl_array returned %s", gsl_strerror(ret));
+ return wrapper;
+ }
+ array_wrapper* gsl_sf_bessel_il_scaled_array_wrapper(int lmax, double x) {
+ int ret;
+ array_wrapper * wrapper = array_wrapper_alloc(lmax + 1, awDouble);
+ ret = gsl_sf_bessel_il_scaled_array(lmax,x, (double*)(wrapper->data));
+ if (ret != GSL_SUCCESS)
+ croak("Math::GSL::SF::gsl_sf_bessel_il_scaled_array returned %s", gsl_strerror(ret));
+ return wrapper;
+ }
+ array_wrapper* gsl_sf_bessel_kl_scaled_array_wrapper(int lmax, double x) {
+ int ret;
+ array_wrapper * wrapper = array_wrapper_alloc(lmax + 1, awDouble);
+ ret = gsl_sf_bessel_kl_scaled_array(lmax,x, (double*)(wrapper->data));
+ if (ret != GSL_SUCCESS)
+ croak("Math::GSL::SF::gsl_sf_bessel_kl_scaled_array returned %s", gsl_strerror(ret));
+ return wrapper;
+ }
+
%}
+
%include "gsl/gsl_types.h"
%include "gsl/gsl_version.h"
%include "gsl/gsl_mode.h"
View
29 swig/gsl_typemaps.i
@@ -11,7 +11,7 @@
%}
/*****************************
- * handle 'double const []' as an input array of doubles
+ * handle 'int const []' as an input array of doubles
* We allocate the C array at the beginning and free it at the end
*/
%typemap(in) int const [] {
@@ -167,6 +167,7 @@
}
}
+
%typemap(freearg) float [] {
// if ($1) free(((char*)$1)-sizeof(struct perl_array));
}
@@ -183,6 +184,7 @@
float *d1, float *d2
};
+
/*****************************
* handle 'size_t const []' as an input array of size_t
* We allocate the C array at the begining and free it at the end
@@ -220,6 +222,7 @@
const double * base, const double xrange[], const double yrange[] ,
const double * array , const double data2[], const double w[] ,
double *v,
+ double result_array[],
gsl_complex_packed_array data
};
@@ -351,6 +354,26 @@ void array_wrapper_free(array_wrapper * daw){
$2 = 0; /* apparently arrays are not null'd for multi-args? */
}
+%typemap(in) (int min, int nmax, double x, double result_array[]) {
+ AV* av;
+ int i;
+
+ if (!SvROK($input))
+ croak("Argument $argnum is not a reference.");
+ if (SvTYPE(SvRV($input)) != SVt_PVAV)
+ croak("Argument $argnum is not an array.");
+ av = (AV*)SvRV($input);
+ $1 = av_len(av) + 1;
+
+ $2 = malloc($1 * sizeof(double));
+ if ($4 == NULL)
+ croak("%%typemap(in) - can't malloc");
+
+ for (i = 0; i < $1; i++) {
+ $4[i] = (double) SvNV(* av_fetch(av, i, 0));
+ }
+}
+
%typemap(in) (size_t SIZE, const double ARRAY[]) {
AV* av;
int i;
@@ -364,7 +387,7 @@ void array_wrapper_free(array_wrapper * daw){
$2 = malloc($1 * sizeof(double));
if ($2 == NULL)
- croak("%typemap(in) (int , const double []) - can't malloc");
+ croak("%%typemap(in) (int , const double []) - can't malloc");
for (i = 0; i < $1; i++) {
$2[i] = (double) SvNV(* av_fetch(av, i, 0));
@@ -393,7 +416,7 @@ void array_wrapper_free(array_wrapper * daw){
$2 = malloc($1 * sizeof(unsigned int));
if ($2 == NULL)
- croak("%typemap(in) (int , unsigned int []) - can't malloc");
+ croak("%%typemap(in) (int , unsigned int []) - can't malloc");
for (i = 0; i < $1; i++) {
$2[i] = (unsigned int) SvUV(* av_fetch(av, i, 0));
View
54 t/SF.t
@@ -3,7 +3,7 @@ use strict;
use warnings;
use Math::GSL::Test qw/:all/;
use base q{Test::Class};
-use Test::Most tests => 1120;
+use Test::Most tests => 1132;
use Math::GSL qw/:all/;
use Math::GSL::Const qw/:all/;
use Math::GSL::Errno qw/:all/;
@@ -1232,4 +1232,56 @@ sub TEST_MATHIEU : Tests(8) {
verify_results($results, 'Math::GSL::SF');
}
+sub TEST_ZBESSEL_ARRAYS : Tests {
+ my $J = gsl_sf_bessel_Jn_array (0, 15, 0);
+ ok_similar( $J, [ 1, (0) x 15 ], "gsl_sf_bessel_Jn_array(0,15,0)");
+
+ $J = gsl_sf_bessel_Jn_array(0,1,1);
+ # values from http://functions.wolfram.com/webMathematica/FunctionEvaluation.jsp?name=BesselJ
+ ok_similar( $J, [ 0.765197686558, 0.440050585745 ], "gsl_sf_bessel_Jn_array(0,1,1)");
+
+ my $K = gsl_sf_bessel_Kn_array(0,1,1);
+ # values from http://functions.wolfram.com/webMathematica/FunctionEvaluation.jsp?name=BesselK
+ ok_similar( $K, [ 0.421024438241, 0.601907230197], "gsl_sf_bessel_Kn_array(0,1,1)");
+
+ $K = gsl_sf_bessel_Kn_scaled_array(0,1,1);
+ ok_similar( $K, [ exp(1)*0.421024438241, exp(1)*0.601907230197], "gsl_sf_bessel_Kn_scaled_array(0,1,1)");
+
+ my $j = gsl_sf_bessel_jl_array(1,1);
+ ok_similar($j, [sin(1), sin(1) - cos(1)], "gsl_sf_bessel_jl_array(0,1,1)");
+
+ $j = gsl_sf_bessel_jl_steed_array(1,1);
+ ok_similar($j, [sin(1), sin(1) - cos(1)], "gsl_sf_bessel_jl_steed_array(0,1,1)");
+
+ my $y = gsl_sf_bessel_yl_array(1,1);
+ ok_similar($y, [-cos(1), -cos(1)-sin(1)], "gsl_sf_bessel_yl_array(0,1,1)");
+
+ dies_ok( sub {
+ my $K = gsl_sf_bessel_Kn_array(0,1,0);
+ }, "gsl_sf_bessel_Kn_array(0,1,0) dies because it is not defined at 0");
+
+ my $I = gsl_sf_bessel_In_array(0,1,1);
+ # values from http://functions.wolfram.com/webMathematica/FunctionEvaluation.jsp?name=BesselI
+ ok_similar( $I, [ 1.26606587775,0.565159103992 ], "gsl_sf_bessel_In_array(0,1,1)");
+
+ my $Is = gsl_sf_bessel_In_scaled_array(0,1,1);
+
+ ok_similar( $Is, [ exp(-1)*1.26606587775,exp(-1)*0.565159103992 ], "gsl_sf_bessel_In_scaled_array(0,1,1)");
+
+ my $Y = gsl_sf_bessel_Yn_array(0,1,1);
+ # values from http://functions.wolfram.com/webMathematica/FunctionEvaluation.jsp?name=BesselY
+ ok_similar( $Y, [ 0.0882569642157, -0.781212821300], "gsl_sf_bessel_Yn_array(0,1,1)");
+
+ lives_ok(sub { my $il = gsl_sf_bessel_il_scaled_array(1,1) }, 'gsl_sf_bessel_il_scaled_array(1,1) lives ');
+
+ return;
+ # kl_scaled_array does not seem to be found by SWIG, although it is for gsl2.0
+ lives_ok(sub {
+ my $kl = gsl_sf_bessel_kl_scaled_array(1,1)
+ }, 'gsl_sf_bessel_kl_scaled_array(1,1) lives ');
+
+ my $kl = gsl_sf_bessel_kl_scaled_array(1,1);
+ die Dumper [ $kl ];
+}
+
Test::Class->runtests;

0 comments on commit 7a7ebf8

Please sign in to comment.
Something went wrong with that request. Please try again.