;; Demo of Spline interpolation using GSL.
;; Derived from an example in chapter 25 of the GSL book.

(libload "gsl/aux_globalvar")
(libload "gsl/aux_structure")
(libload "gsl/interpolation")

;; synthesize 10 sample points
(defparameter numpoints 10)
(defparameter x (double-matrix numpoints))
(defparameter y (double-matrix numpoints))
(for (i 0 (- numpoints 1))
     (x i (+ i (* 0.5 (sin i))))
     (y i (+ i (cos (* i i)))))

;; build spline evaluator from sample points
(defparameter acc (gsl_interp_accel_alloc ))
(defparameter spline (gsl_spline_alloc  (var_gsl_interp_cspline) numpoints))
(gsl_spline_init spline (idx-ptr x) (idx-ptr y) numpoints)

;; evaluate and plot spline
(let* ((xx (range 0.1 (- numpoints 0.5) 0.1))
       (yy (all ((x xx)) (gsl_spline_eval spline x acc))))
  ;; plot curve
  (defparameter plot-context 
    (graph-xy xx yy 
	      '(title "GSL Spline Interpolation")
	      '(xgrid t) '(ygrid t)
	      '(color-rgb 1 0 0)
	      '(object object-nil)))
  ;; plot control points
  (graph-xy (x ()) (y ())
	    plot-context
	    '(color-rgb 0 0 0.7)
	    '(linestyle 1)
	    '(object closed-circle)
	    '(object-size 5)))

;; free stuff
(gsl_spline_free  spline)
(gsl_interp_accel_free acc)