source: project/release/4/sundials/trunk/tests/iaf.scm @ 25539

;;
;; Simple integrate-and-fire example to illustrate integrating over
;; discontinuities.
;;

(use sundials srfi-4)


(define  reltol 1.0e-3)
(define  abstol 1.0e-4)


(define  t0 0.0)
(define  tf 100.0)



(define gL        0.2 )
(define vL      -70.0 )
(define Isyn     20.0 )
(define C         1.0 )
(define theta   -25.0 )
(define vreset  -65.0 )
(define trefractory  5.0 )


(define (subthreshold t yy)
  (let ((v (f64vector-ref yy 0)))
    (f64vector (/ (+ (* (- gL) (- v vL)) Isyn) C))))

(define (spike-detect t yy)
  (let ((v (f64vector-ref yy 0)))
    (s32vector (if (>= v theta) 0 1))))
  
    

(define (main)
  
  (let ((yy (f64vector -65.0)))

    ;; CVODE initialization 
    (let ((dt 1e-2)
          (solver (cvode-create-solver
                   t0 yy subthreshold  
                   tstop: tf
                   abstol: abstol
                   reltol: reltol
                   lmm: cvode-lmm/bdf
                   events: (s32vector 0)
                   event-fn: spike-detect
                   )))
      
      ;; 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 (= 1 flag)
              ;; discrete event: re-initialize solver to integrate over the discontinuity
              (cvode-reinit-solver solver (+ tnext trefractory) (f64vector vreset)) )
          
          ;;(print-results solver tnext)

          (if (< tnext tf)
              (recur (if (= 1 flag) (+ (+ tnext trefractory) dt) (+ tnext dt)) (+ 1 iout)))
          ))
      
      )))


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

      
(main)