Commit 078abffa authored by Liam Healy's avatar Liam Healy

Interpolation fixes

Placed missing #'mpointer calls, added :input declarations in
defmfuns.  Fixed #'spline-example to work correctly with marrays, and
added a test.
parent 86da4ec7
;; Definition of GSLL system
;; Liam Healy
;; Time-stamp: <2009-01-22 21:37:49EST gsll-tests.asd>
;; Time-stamp: <2009-01-26 22:45:42EST gsll-tests.asd>
;; $Id$
(asdf:defsystem "gsll-tests"
......@@ -74,6 +74,7 @@
(:file "hypergeometric")
(:file "hypergeometric-randist")
(:file "index-max")
(:file "interpolation")
(:file "inverse-matrix-product")
(:file "laguerre")
(:file "lambert")
......
......@@ -5,8 +5,9 @@
(defmfun evaluate-interpolation (interpolation xa ya x acceleration)
"gsl_interp_eval"
((interpolation :pointer) (xa :pointer) (ya :pointer) (x :double)
(acceleration :pointer))
(((mpointer interpolation) :pointer) (xa :pointer) (ya :pointer) (x :double)
((mpointer acceleration) :pointer))
:inputs (xa ya)
:c-return :double
:documentation ; FDL
"Find the interpolated value of y for a given
......@@ -15,8 +16,9 @@
(defmfun derivative-evaluate-interpolation (interpolation xa ya x acceleration)
"gsl_interp_eval_deriv"
((interpolation :pointer) (xa :pointer) (ya :pointer) (x :double)
(acceleration :pointer))
(((mpointer interpolation) :pointer) (xa :pointer) (ya :pointer) (x :double)
((mpointer acceleration) :pointer))
:inputs (xa ya)
:c-return :double
:documentation ; FDL
"Find the derivative of an interpolated function for a given point
......@@ -25,8 +27,9 @@
(defmfun second-derivative-evaluate-interpolation (interpolation xa ya x acceleration)
"gsl_interp_eval_deriv2"
((interpolation :pointer) (xa :pointer) (ya :pointer) (x :double)
(acceleration :pointer))
(((mpointer interpolation) :pointer) (xa :pointer) (ya :pointer) (x :double)
((mpointer acceleration) :pointer))
:inputs (xa ya)
:c-return :double
:documentation ; FDL
"Find the second derivative of an interpolated function for a given point
......@@ -35,8 +38,9 @@
(defmfun integral-evaluate-interpolation (interpolation xa ya x low high acceleration)
"gsl_interp_eval_integ"
((interpolation :pointer) (xa :pointer) (ya :pointer) (x :double)
(low :double) (high :double) (acceleration :pointer))
(((mpointer interpolation) :pointer) (xa :pointer) (ya :pointer) (x :double)
(low :double) (high :double) ((mpointer acceleration) :pointer))
:inputs (xa ya)
:c-return :double
:documentation ; FDL
"Find the numerical integral of an interpolated function over the
......@@ -46,7 +50,7 @@
;;; Spline
(defmfun evaluate-spline (spline x acceleration)
"gsl_spline_eval"
((spline :pointer) (x :double) (acceleration :pointer))
(((mpointer spline) :pointer) (x :double) ((mpointer acceleration) :pointer))
:c-return :double
:documentation
"Find the interpolated value of y for a given
......@@ -55,7 +59,7 @@
(defmfun derivative-evaluate-spline (spline x acceleration)
"gsl_spline_eval_deriv"
((spline :pointer) (x :double) (acceleration :pointer))
(((mpointer spline) :pointer) (x :double) ((mpointer acceleration) :pointer))
:c-return :double
:documentation ; FDL
"Find the derivative of an interpolated function for a given point
......@@ -64,7 +68,7 @@
(defmfun second-derivative-evaluate-spline (spline x acceleration)
"gsl_spline_eval_deriv2"
((spline :pointer) (x :double) (acceleration :pointer))
(((mpointer spline) :pointer) (x :double) ((mpointer acceleration) :pointer))
:c-return :double
:documentation ; FDL
"Find the second derivative of an interpolated function for a given point
......@@ -73,8 +77,8 @@
(defmfun integral-evaluate-spline (spline x low high acceleration)
"gsl_spline_eval_integ"
((spline :pointer) (x :double) (low :double) (high :double)
(acceleration :pointer))
(((mpointer spline) :pointer) (x :double) (low :double) (high :double)
((mpointer acceleration) :pointer))
:c-return :double
:documentation ; FDL
"Find the numerical integral of an interpolated function over the
......
;; Interpolation allocation, initialization, and freeing.
;; Liam Healy, Sun Nov 4 2007 - 17:24
;; Time-stamp: <2009-01-25 10:14:30EST interpolation.lisp>
;; Time-stamp: <2009-01-26 22:51:09EST interpolation.lisp>
;; $Id$
(in-package :gsl)
......@@ -37,6 +37,7 @@
:initialize-suffix "init"
:initialize-args
(((c-pointer xa) :pointer) ((c-pointer ya) :pointer) ((dim0 xa) sizet))
:inputs (xa ya)
:arglists-function
(lambda (set)
`((type &optional xa-or-size (ya nil ,set))
......
;; Example spline
;; Liam Healy, Sat Nov 10 2007 - 21:18
;; Time-stamp: <2008-12-26 12:53:23EST spline-example.lisp>
;; Time-stamp: <2009-01-26 22:43:41EST spline-example.lisp>
;; $Id$
(in-package :gsl)
;;; Possible future improvement: direct double-float vectors.
(defun spline-example-arrays ()
(let ((xarr (make-array '(10) :element-type 'double-float))
(yarr (make-array '(10) :element-type 'double-float)))
(loop for i from 0 below 10
do (setf (aref xarr i) (+ i (* 0.5d0 (sin (coerce i 'double-float))))
(aref yarr i) (+ i (cos (expt (coerce i 'double-float) 2)))))
(values xarr yarr)))
(defun spline-example ()
(defun spline-example (&optional (step 0.01d0))
"The first example in Sec. 26.7 of the GSL manual."
(multiple-value-bind (xarr yarr)
(spline-example-arrays)
(let* ((acc (make-acceleration))
(cxarr (make-marray 'double-float :dimensions xarr))
(cyarr (make-marray 'double-float :dimensions yarr))
(spline (make-spline *cubic-spline-interpolation* cxarr cyarr)))
(loop for xi from (aref xarr 0) below (aref xarr 9) by 0.01d0
collect (list xi (evaluate-spline spline xi acc))))))
(let* ((acc (make-acceleration))
(xarr
(make-marray
'double-float
:initial-contents
(loop for i from 0.0d0 below 10.0d0
collect (+ i (* 0.5d0 (sin i))))))
(yarr
(make-marray
'double-float
:initial-contents
(loop for i from 0.0d0 below 10.0d0
collect (+ i (cos (expt i 2))))))
(spline (make-spline *cubic-spline-interpolation* xarr yarr)))
(loop for xi from (maref xarr 0) below (maref xarr 9) by step
collect (list xi (evaluate-spline spline xi acc)))))
(save-test interpolation (spline-example 0.1d0))
;; Regression test INTERPOLATION for GSLL, automatically generated
(in-package :gsl)
(LISP-UNIT:DEFINE-TEST INTERPOLATION
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST
(LIST (LIST 0.0d0 1.0d0)
(LIST 0.1d0 1.069722939658712d0)
(LIST 0.2d0 1.1384991026983018d0)
(LIST 0.30000000000000004d0
1.2053817124996473d0)
(LIST 0.4d0 1.2694239924436264d0)
(LIST 0.5d0 1.329679165911117d0)
(LIST 0.6d0 1.3852004562829967d0)
(LIST 0.7d0 1.4350410869401438d0)
(LIST 0.7999999999999999d0 1.4782542812634358d0)
(LIST 0.8999999999999999d0 1.5138932626337507d0)
(LIST 0.9999999999999999d0 1.5410112544319663d0)
(LIST 1.0999999999999999d0 1.5586614800389607d0)
(LIST 1.2d0 1.5658971628356115d0)
(LIST 1.3d0 1.5617715262027967d0)
(LIST 1.4000000000000001d0 1.5453377935213939d0)
(LIST 1.5000000000000002d0 1.5160839349736626d0)
(LIST 1.6000000000000003d0 1.4767879497386507d0)
(LIST 1.7000000000000004d0 1.4317331076494486d0)
(LIST 1.8000000000000005d0 1.3852104614697154d0)
(LIST 1.9000000000000006d0 1.3415110639631103d0)
(LIST 2.0000000000000004d0 1.3049259678932923d0)
(LIST 2.1000000000000005d0 1.2797462260239207d0)
(LIST 2.2000000000000006d0 1.2702628911186542d0)
(LIST 2.3000000000000007d0 1.2807670159411522d0)
(LIST 2.400000000000001d0 1.3155496532550734d0)
(LIST 2.500000000000001d0 1.3787386279170017d0)
(LIST 2.600000000000001d0 1.4697408675356218d0)
(LIST 2.700000000000001d0 1.5826333311720413d0)
(LIST 2.800000000000001d0 1.7112073731538484d0)
(LIST 2.9000000000000012d0 1.849254347808631d0)
(LIST 3.0000000000000013d0 1.9905656094639774d0)
(LIST 3.1000000000000014d0 2.1290059431234454d0)
(LIST 3.2000000000000015d0 2.2643876349412406d0)
(LIST 3.3000000000000016d0 2.406758011312636d0)
(LIST 3.4000000000000017d0 2.567175372846213d0)
(LIST 3.5000000000000018d0 2.7566980201505533d0)
(LIST 3.600000000000002d0 2.9863842538342387d0)
(LIST 3.700000000000002d0 3.2653444558817344d0)
(LIST 3.800000000000002d0 3.587530003184208d0)
(LIST 3.900000000000002d0 3.939787612326211d0)
(LIST 4.000000000000002d0 4.308923272380094d0)
(LIST 4.100000000000001d0 4.681742972418209d0)
(LIST 4.200000000000001d0 5.045052701512908d0)
(LIST 4.300000000000001d0 5.385658448736544d0)
(LIST 4.4d0 5.690366203161469d0)
(LIST 4.5d0 5.945981953860034d0)
(LIST 4.6d0 6.14101013300902d0)
(LIST 4.699999999999999d0 6.27672680272461d0)
(LIST 4.799999999999999d0 6.360219169852809d0)
(LIST 4.899999999999999d0 6.398603766104035d0)
(LIST 4.999999999999998d0 6.398997123188704d0)
(LIST 5.099999999999998d0 6.3685157728172355d0)
(LIST 5.1999999999999975d0 6.314276246700047d0)
(LIST 5.299999999999997d0 6.243395076547555d0)
(LIST 5.399999999999997d0 6.16298879407018d0)
(LIST 5.4999999999999964d0 6.0801739309783365d0)
(LIST 5.599999999999996d0 6.002067018982444d0)
(LIST 5.699999999999996d0 5.935784589792921d0)
(LIST 5.799999999999995d0 5.888443175120185d0)
(LIST 5.899999999999995d0 5.867053973966079d0)
(LIST 5.999999999999995d0 5.874461774747944d0)
(LIST 6.099999999999994d0 5.908057260357294d0)
(LIST 6.199999999999994d0 5.964862372253499d0)
(LIST 6.299999999999994d0 6.041899051895926d0)
(LIST 6.399999999999993d0 6.136189240743944d0)
(LIST 6.499999999999993d0 6.244754880256921d0)
(LIST 6.5999999999999925d0 6.364617911894225d0)
(LIST 6.699999999999992d0 6.492800277115226d0)
(LIST 6.799999999999992d0 6.626323917379292d0)
(LIST 6.8999999999999915d0 6.762210774145791d0)
(LIST 6.999999999999991d0 6.897482788874091d0)
(LIST 7.099999999999991d0 7.02916190302356d0)
(LIST 7.19999999999999d0 7.154270058053569d0)
(LIST 7.29999999999999d0 7.269829195423485d0)
(LIST 7.39999999999999d0 7.373218737499652d0)
(LIST 7.499999999999989d0 7.4653205559805d0)
(LIST 7.599999999999989d0 7.549000269792069d0)
(LIST 7.699999999999989d0 7.627146115165377d0)
(LIST 7.799999999999988d0 7.702646328331444d0)
(LIST 7.899999999999988d0 7.778389145521288d0)
(LIST 7.999999999999988d0 7.857262802965928d0)
(LIST 8.099999999999987d0 7.942155536896381d0)
(LIST 8.199999999999987d0 8.035955583543666d0)
(LIST 8.299999999999986d0 8.141551179138803d0)
(LIST 8.399999999999986d0 8.261830559912811d0)
(LIST 8.499999999999986d0 8.39968181791184d0)
(LIST 8.599999999999985d0 8.556875437368307d0)
(LIST 8.699999999999985d0 8.731369216108577d0)
(LIST 8.799999999999985d0 8.920308621429049d0)
(LIST 8.899999999999984d0 9.120839120626123d0)
(LIST 8.999999999999984d0 9.3301061809962d0)
(LIST 9.099999999999984d0 9.545255269835682d0)
(LIST 9.199999999999983d0 9.763431854440967d0)))
(MULTIPLE-VALUE-LIST (SPLINE-EXAMPLE 0.1d0))))
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment