source: project/release/4/sundials/trunk/tests/ml.scm @ 25082

;;
;; Morris-Lecar model
;;

(use mathh sundials srfi-4)


(define TEND  500.0)

                           
;; Model parameters

(define Istim  20)
(define vk     -84)
(define vl     -60)
(define vca    120)
(define gk     2.5)
(define gl     0.5)
(define gca    1)
(define c      2)
(define v1     -1.2)
(define v2     18)
(define v3     12)
(define v4     17)
(define phi    0.0667)


;; State functions

(define (minf v) (* 0.5 (+ 1 (tanh (/ (- v v1) v2)))))
(define (winf v) (* 0.5 (+ 1 (tanh (/ (- v v3) v4)))))
(define (lamw v) (* phi (cosh (/ (- v v3) (* 2 v4)))))
                           
;; Model equations

(define (rhs t yy)
  (let ((v (f64vector-ref yy 0))
        (w (f64vector-ref yy 1)))

  (let ((ica (* gca (* (minf v)  (- vca v))))
        (ik  (* gk  (* w (- vk v )))))
    
    (let ((dv (/ (+ Istim (* gl (- vl v)) ica ik) c))
          (dw (* (lamw v) (- (winf v) w))))
    
      (f64vector dv dw)

      ))
  ))

(define (rhs/unsafe t yy yp _)
  (let ((v (pointer-f64-ref yy))
        (w (pointer-f64-ref (pointer+-f64 yy 1))))

  (let ((ica (* gca (* (minf v)  (- vca v))))
        (ik  (* gk  (* w (- vk v )))))
    
    (let ((dv (/ (+ Istim (* gl (- vl v)) ica ik) c))
          (dw (* (lamw v) (- (winf v) w))))
    
      (pointer-f64-set! yp dv)
      (pointer-f64-set! (pointer+-f64 yp 1) dw)

      ))
  ))



(define (main)
  
  (let ((yy (f64vector -60.899 0.0149)) ;; v w

        ;; Integration limits 
        (t0  0.0)
        (tf  TEND)
        (dt  1e-2))
    
    ;; CVODE initialization 
    (let ((solver (cvode-create-solver
                   t0 tf yy rhs
                   abstol: 1e-4
                   reltol: 1e-4)))

      ;; In loop, call CVodeSolve, print results, and test for error. 

      (let recur ((tnext (+ t0 dt)) (iout 1))

        (let ((flag  (cvode-solve solver tnext)))
          (if (negative? flag) (error 'main "CVODE solver error" flag))

          (if (< tnext tf)
              (recur (+ tnext dt) (+ 1 iout)))
          ))
      
      (let ((yy (cvode-yy solver)))
        (let ((v (f64vector-ref yy 0))
              (w (f64vector-ref yy 1)))
          (assert (< (abs (- v -30.0536210523504)) 1e-12))
          (assert (< (abs (- w 0.0402640532124851 )) 1e-12))))

      (cvode-destroy-solver solver)
      
      )))

(define (main/unsafe)
  
  (let ((yy (f64vector -60.899 0.0149)) ;; v w

        ;; Integration limits 
        (t0  0.0)
        (tf  TEND)
        (dt  1e-2))
    
    ;; CVODE initialization 
    (let ((solver (cvode-create-solver/unsafe
                   t0 tf yy rhs/unsafe
                   abstol: 1e-4
                   reltol: 1e-4)))

      ;; In loop, call CVodeSolve, print results, and test for error. 

      (let recur ((tnext (+ t0 dt)) (iout 1))

        (let ((flag  (cvode-solve solver tnext)))
          (if (negative? flag) (error 'main "CVODE solver error" flag))

          (if (< tnext tf)
              (recur (+ tnext dt) (+ 1 iout)))
          ))
      
      (let ((yy (cvode-yy solver)))
        (let ((v (f64vector-ref yy 0))
              (w (f64vector-ref yy 1)))
          (assert (< (abs (- v -30.0536210523504)) 1e-12))
          (assert (< (abs (- w 0.0402640532124851 )) 1e-12))))

      (cvode-destroy-solver solver)
      
      )))


(define (print-results solver t)
  (let ((yy (cvode-yy solver)))
    (printf "~A ~A ~A ~A ~A ~A~%" 
            t 
            (f64vector-ref yy 0)
            (f64vector-ref yy 1)
            (cvode-get-last-order solver)
            (cvode-get-num-steps solver)
            (cvode-get-last-step solver)
            )))
      
      
(main)
(main/unsafe)