hankel.lisp 3.9 KB
Newer Older
1 2
;; Discrete Hankel Transforms.
;; Liam Healy, Sat Dec  8 2007 - 16:50
3
;; Time-stamp: <2014-12-26 13:18:38EST hankel.lisp>
4
;;
5
;; Copyright 2007, 2008, 2009, 2011, 2012, 2014 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 committed
20
;; $Id$
21 22

(in-package :gsl)
23
(named-readtables:in-readtable :antik)
24

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

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

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

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

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

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


(save-test hankel					; tests from dht/test.c
73
 (grid:copy-to (apply-hankel (make-hankel 3 1.0d0 1.0d0) #m(1.0d0 2.0d0 3.0d0)))
74 75 76 77 78 79 80 81 82
 ;; 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))
83
       (in (grid:make-foreign-array 'double-float :dimensions 128)))
84
   (loop for n from 0 below 128
85
      do (setf (grid:aref in n)
86
	       (/ (1+ (expt (sample-x-hankel hank n) 2)))))
87
   (grid:copy-to (apply-hankel hank in)))
88 89 90
 ;; 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))
91
       (in (grid:make-foreign-array 'double-float :dimensions 128)))
92
   (loop for n from 0 below 128
93
      do (setf (grid:aref in n) (exp (- (sample-x-hankel hank n)))))
94
   (grid:copy-to (apply-hankel hank in)))
95 96
 ;; 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))
97
       (in (grid:make-foreign-array 'double-float :dimensions 128)))
98
   (loop for n from 0 below 128
99
      do (setf (grid:aref in n)
100 101
	       (let ((x (sample-x-hankel hank n)))
		 (* x (- 1 (expt x 2))))))
102
   (grid:copy-to (apply-hankel hank in))))