SO31.ss
#|  SO31.ss: FFI Interface to the SO(3,1) C library.
    Copyright (C) 2007 Will M. Farr <farr@mit.edu>

    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 2 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, write to the Free Software Foundation, Inc.,
    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
|#

(module SO31 mzscheme
  (require (planet "all.ss" ("wmfarr" "plt-linalg.plt" 1 1))
           (all-except (lib "contract.ss") ->)
           (rename (lib "contract.ss") ->/c ->)
           (lib "foreign.ss")
           (all-except (lib "42.ss" "srfi") :)
           (lib "etc.ss"))
  
  (define (SO31-matrix? obj)
    (and (matrix? obj)
         (= (matrix-rows obj) 4)
         (= (matrix-cols obj) 4)))
  
  (define (SO31-params? obj)
    (and (f64vector? obj)
         (= (f64vector-length obj) 6)))
  
  (provide SO31-matrix? SO31-params?)
  
  (provide/contract
   (Rx (->/c number? SO31-matrix?))
   (Ry (->/c number? SO31-matrix?))
   (Rz (->/c number? SO31-matrix?))
   (Bx (->/c number? SO31-matrix?))
   (By (->/c number? SO31-matrix?))
   (Bz (->/c number? SO31-matrix?))
   (matrix->params (->/c SO31-matrix? SO31-params?))
   (params->matrix (->/c SO31-params? SO31-matrix?))
   (params->inverse-params (->/c SO31-params? SO31-params?))
   (params-compose (->/c SO31-params? SO31-params? SO31-params?))
   (params->dM (->/c SO31-params? (vectorof SO31-matrix?)))
   (matrix->dM (->/c SO31-matrix? (vectorof SO31-matrix?))))
  
  (define *lib* (ffi-lib (build-path (this-expression-source-directory) "SO31c" "compiled" "native" (system-library-subpath) "SO31")))
  
  (unsafe!) (matrix-unsafe!)
  
  (define-syntax (generate-matrix-bindings stx)
    (syntax-case stx ()
      ((generate-matrix-bindings name)
       (with-syntax ((C-name (datum->syntax-object
                              (syntax name)
                              (string->symbol
                               (string-append "lorentz_"
                                              (symbol->string (syntax-object->datum (syntax name))))))))
         (syntax
          (define name
            (get-ffi-obj (quote C-name) *lib*
                         (_fun _double*
                               (m : (_matrix o 4 4)) ->
                               _void ->
                               m))))))))
  
  (generate-matrix-bindings Rx)
  (generate-matrix-bindings Ry)
  (generate-matrix-bindings Rz)
  (generate-matrix-bindings Bx)
  (generate-matrix-bindings By)
  (generate-matrix-bindings Bz)
  
  (define matrix->params
    (get-ffi-obj 'lorentz_params_from_matrix *lib*
                 (_fun _matrix
                       (p : (_f64vector o 6)) ->
                       _void ->
                       p)))
  
  (define params->matrix
    (get-ffi-obj 'lorentz_matrix_from_params *lib*
                 (_fun _f64vector
                       (m : (_matrix o 4 4)) ->
                       _void ->
                       m)))
  
  (define params->inverse-params
    (get-ffi-obj 'params_to_inverse_params *lib*
                 (_fun _f64vector
                       (pinv : (_f64vector o 6)) ->
                       _void ->
                       pinv)))
  
  (define params-compose
    (get-ffi-obj 'params_compose *lib*
                 (_fun _f64vector _f64vector
                       (pc : (_f64vector o 6)) ->
                       _void ->
                       pc)))
  
  (define (make-dM-storage)
    (let ((v (make-cvector _pointer 6)))
      (do-ec (:range i 6) (cvector-set! v i (malloc _double 16 'atomic)))
      v))
  
  (define params->dM
    (get-ffi-obj 'lorentz_deriv_matrix_at_params *lib*
                 (_fun _f64vector
                       (dMs : _cvector = (make-dM-storage)) ->
                       _void ->
                       (vector-of-length-ec 6 (:range i 6)
                         (ptr->matrix (cvector-ref dMs i) 4 4)))))
  
  (define matrix->dM
    (get-ffi-obj 'lorentz_deriv_matrix_at_matrix *lib*
                 (_fun _matrix
                       (dMs : _cvector = (make-dM-storage)) ->
                       _void ->
                       (vector-of-length-ec 6 (:range i 6)
                         (ptr->matrix (cvector-ref dMs i) 4 4))))))