;;; PLT Scheme Science Collection
;;; chebyshev.ss; version 0.1.0
;;; Copyright (c) 2004 M. Douglas Williams
;;;
;;; This library is free software; you can redistribute it and/or
;;; modify it under the terms of the GNU Lesser General Public
;;; License as published by the Free Software Foundation; either
;;; version 2.1 of the License, or (at your option) any later version.
;;;
;;; This library 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
;;; Lesser General Public License for more details.
;;;
;;; You should have received a copy of the GNU Lesser General Public
;;; License along with this library; if not, write to the Free
;;; Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
;;; 02111-1307 USA.
;;;
;;; -------------------------------------------------------------------
;;;
;;; This code in based on the Chebyshev Approximations in the GNU
;;; Scientific Library (GSL), which is licensed the GPL.
;;;
;;; Version  Date      Description
;;; 0.1.0    08/07/04  This is the initial release of the Chebyshev
;;;                    approximations routines ported from GSL. (Doug
;;;                    Williams)
;;; 1.0.0    09/29/04  Added contracts for functions.  Marked as ready
;;;                    for Release 1.0.  (Doug Williams)

(module chebyshev mzscheme
  
  (require (lib "contract.ss"))
  
  ;; chebyshev-series: struct
  (define-struct chebyshev-series
                 (coefficients
                  order
                  lower
                  upper))
  
  ;; Define the contracts for the module.
  (provide/contract
   (struct chebyshev-series ((coefficients (vectorof real?))
                             (order natural-number/c)
                             (lower real?)
                             (upper real?)))
   (make-chebyshev-series-order
    (-> natural-number/c chebyshev-series?))
   (chebyshev-series-init
    (-> chebyshev-series? procedure? real? real? void?))
   (chebyshev-eval
    (-> chebyshev-series? real? real?))
   (chebyshev-eval-n
    (-> chebyshev-series? natural-number/c real? real?)))
  
  (require "math.ss")
  
  ;; make-chebyshev-series-order: real -> chebyshev-series
  (define (make-chebyshev-series-order order)
    (make-chebyshev-series
     (make-vector (+ order 1))
     order
     0.0 0.0))
  
  ;; chebyshev-series-init: chebyshev-series x function x real x real
  ;;   -> void
  (define (chebyshev-series-init cs func a b)
    (if (>= a b)
        (error 'chebyshev-series-init
               "null function interval [a,b]"))
    (set-chebyshev-series-lower! cs a)
    (set-chebyshev-series-upper! cs b)
    (let* ((order (chebyshev-series-order cs))
           (coefficients (chebyshev-series-coefficients cs))
           (bma (* 0.5 (- b a)))
           (bpa (* 0.5 (+ b a)))
           (fac (/ 2.0 (+ order 1.0)))
           (f (make-vector (+ order 1))))
      (do ((k 0 (+ k 1)))
          ((> k order) (void))
        (let ((y (cos (* pi (/ (+ k 0.5) (+ order 1))))))
          (vector-set! f k
                       (func (+ (* y bma) bpa)))))
      (do ((j 0 (+ j 1)))
          ((> j order) (void))
        (let ((sum 0.0))
          (do ((k 0 (+ k 1)))
              ((> k order) (void))
            (set! sum (+ sum (* (vector-ref f k)
                                (cos (/ (* pi j (+ k 0.5))
                                        (+ order 1)))))))
          (vector-set! coefficients j (* fac sum))))))
      
  ;; chebyshev-eval: chebyshev-series x real -> real
  ;;
  ;; Evaluate the Chebyshev series at the given point.
  (define (chebyshev-eval cs x)
    (let* ((coefs (chebyshev-series-coefficients cs))
           (order (chebyshev-series-order cs))
           (lower (chebyshev-series-lower cs))
           (upper (chebyshev-series-upper cs))
           (d 0.0)
           (dd 0.0)
           (y (/ (- (* 2.0 x) lower upper) (- upper lower)))
           (y2 (* 2.0 y)))
      (do ((j order (- j 1)))
        ((= j 0) (void))
        (let ((temp d))
          (set! d (+ (* y2 d) (- dd) (vector-ref coefs j)))
          (set! dd temp)))
      (+ (* y d) (- dd) (* 0.5 (vector-ref coefs 0)))))
  
  ;; chebyshev-evaln: chebyshev-series x integer x real -> real
  ;;
  ;; Evaluate the Chebyshev series at the given point.
  (define (chebyshev-eval-n cs n x)
    (let* ((coefs (chebyshev-series-coefficients cs))
           (order (chebyshev-series-order cs))
           (lower (chebyshev-series-lower cs))
           (upper (chebyshev-series-upper cs))
           (d 0.0)
           (dd 0.0)
           (y (/ (- (* 2.0 x) lower upper) (- upper lower)))
           (y2 (* 2.0 y)))
      (do ((j (min n order) (- j 1)))
        ((= j 0) (void))
        (let ((temp d))
          (set! d (+ (* y2 d) (- dd) (vector-ref coefs j)))
          (set! dd temp)))
      (+ (* y d) (- dd) (* 0.5 (vector-ref coefs 0)))))
  
)