hankel.lisp 3.85 KB
Newer Older
1 2
;; Discrete Hankel Transforms.
;; Liam Healy, Sat Dec  8 2007 - 16:50
3
;; Time-stamp: <2011-05-26 12:37:29EDT hankel.lisp>
4
;;
5
;; Copyright 2007, 2008, 2009, 2011 Liam M. Healy
6 7 8 9 10 11 12 13 14 15 16 17 18 19
;; Distributed under the terms of the GNU General Public License
;;
;; This program is free software: you can redistribute it and/or modify
;; it under the terms of the GNU General Public License as published by
;; the Free Software Foundation, either version 3 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 General Public License for more details.
;;
;; You should have received a copy of the GNU General Public License
;; along with this program.  If not, see <http://www.gnu.org/licenses/>.
lhealy's avatar
lhealy committed
20
;; $Id$
21 22 23

(in-package :gsl)

24 25
;;; /usr/include/gsl/gsl_dht.h

26
;;; Everything compiles, but not tested -- need example.
27
;;; We don't use gsl_dht_new = gsl_dht_alloc + gsl_dht_init.
28

Liam Healy's avatar
Liam Healy committed
29 30 31 32
;;; Create the Hankel transform object and essential methods/functions.
(defmobject hankel "gsl_dht"
  ((size sizet))
  "discrete Hankel Transform"
Liam Healy's avatar
Liam Healy committed
33
  :documentation
34 35
  "Allocate a Discrete Hankel transform object of the given size and
   optionally initialize the transform for the given values of nu and x."
Liam Healy's avatar
Liam Healy committed
36 37
  :initialize-suffix "init"
  :initialize-args ((nu :double) (xmax :double)))
38

39 40 41
(defmfun apply-hankel
    (hankel array-in
     &optional
42
     (array-out (grid:make-foreign-array 'double-float :dimensions (grid:dimensions array-in))))
43
  "gsl_dht_apply"
44 45
  (((mpointer hankel) :pointer) ((grid:foreign-pointer array-in) :pointer)
   ((grid:foreign-pointer array-out) :pointer))
46 47 48
  :inputs (array-in)
  :outputs (array-out)
  :return (array-out)
liam's avatar
liam committed
49 50
  :documentation			; FDL
  "Apply the transform to the array array-in
51
   whose size is equal to the size of the transform.  The result is stored
liam's avatar
liam committed
52
   in the array array-out which must be of the same length.")
53

liam's avatar
liam committed
54
(defmfun sample-x-hankel (hankel n)
55
  "gsl_dht_x_sample"
56
  (((mpointer hankel) :pointer) (n :int))
57
  :c-return :double
liam's avatar
liam committed
58
  :documentation			; FDL
59
  "The value of the n-th sample point in the unit
60
  interval, (j_{\nu,n+1}/j_{\nu,M}) X, which are the
61
  points where the function f(t) is assumed to be sampled.")
62

liam's avatar
liam committed
63
(defmfun sample-k-hankel (hankel n)
64
  "gsl_dht_k_sample"
65
  (((mpointer hankel) :pointer) (n :int))
66
  :c-return :double
liam's avatar
liam committed
67 68
  :documentation			; FDL
  "The value of the n-th sample point in k-space, j_{\nu,n+1}/X}.")
69 70 71


(save-test hankel					; tests from dht/test.c
72
 (grid:copy-to (apply-hankel (make-hankel 3 1.0d0 1.0d0) #m(1.0d0 2.0d0 3.0d0)))
73 74 75 76 77 78 79 80 81
 ;; Exact, forward-inverse transform should be accurate to 2e-5
 (let ((hank (make-hankel 3 1.0d0 1.0d0)))
   (copy
    (elt* (expt (bessel-zero-j1 4) 2)	; inverse transformation
	  (apply-hankel
	   hank (apply-hankel hank #m(1.0d0 2.0d0 3.0d0))))
    'array))
 ;; Simple
 (let ((hank (make-hankel 128 0.0d0 100.0d0))
82
       (in (grid:make-foreign-array 'double-float :dimensions 128)))
83
   (loop for n from 0 below 128
84
      do (setf (grid:aref in n)
85
	       (/ (1+ (expt (sample-x-hankel hank n) 2)))))
86
   (grid:copy-to (apply-hankel hank in)))
87 88 89
 ;; Integrate[ x exp(-x) J_1(a x), {x,0,Inf}] = a F(3/2, 2; 2; -a^2)
 ;; expected accuracy only about 2%
 (let ((hank (make-hankel 128 1.0d0 20.0d0))
90
       (in (grid:make-foreign-array 'double-float :dimensions 128)))
91
   (loop for n from 0 below 128
92
      do (setf (grid:aref in n) (exp (- (sample-x-hankel hank n)))))
93
   (grid:copy-to (apply-hankel hank in)))
94 95
 ;; Integrate[ x^2 (1-x^2) J_1(a x), {x,0,1}] = 2/a^2 J_3(a)
 (let ((hank (make-hankel 128 1.0d0 1.0d0))
96
       (in (grid:make-foreign-array 'double-float :dimensions 128)))
97
   (loop for n from 0 below 128
98
      do (setf (grid:aref in n)
99 100
	       (let ((x (sample-x-hankel hank n)))
		 (* x (- 1 (expt x 2))))))
101
   (grid:copy-to (apply-hankel hank in))))