source: project/release/4/sundials/tags/1.4/tests/iaf.scm @ 25541

Last change on this file since 25541 was 25541, checked in by Ivan Raikov, 9 years ago

sundials release 1.4

File size: 1.7 KB
Line 
1
2;;
3;; Simple integrate-and-fire example to illustrate integrating over
4;; discontinuities.
5;;
6
7(use sundials srfi-4)
8
9
10(define  reltol 1.0e-3)
11(define  abstol 1.0e-4)
12
13
14(define  t0 0.0)
15(define  tf 100.0)
16
17
18
19(define gL        0.2 )
20(define vL      -70.0 )
21(define Isyn     20.0 )
22(define C         1.0 )
23(define theta   -25.0 )
24(define vreset  -65.0 )
25(define trefractory  5.0 )
26
27
28(define (subthreshold t yy)
29  (let ((v (f64vector-ref yy 0)))
30    (f64vector (/ (+ (* (- gL) (- v vL)) Isyn) C))))
31
32(define (spike-detect t yy)
33  (let ((v (f64vector-ref yy 0)))
34    (s32vector (if (>= v theta) 0 1))))
35 
36   
37
38(define (main)
39 
40  (let ((yy (f64vector -65.0)))
41
42    ;; CVODE initialization
43    (let ((dt 1e-2)
44          (solver (cvode-create-solver
45                   t0 yy subthreshold 
46                   tstop: tf
47                   abstol: abstol
48                   reltol: reltol
49                   lmm: cvode-lmm/bdf
50                   events: (s32vector 0)
51                   event-fn: spike-detect
52                   )))
53     
54      ;; In loop, call CVodeSolve, print results, and test for error.
55     
56      (let recur ((tnext (+ t0 dt)) (iout 1))
57
58        (let ((flag  (cvode-solve solver tnext)))
59
60          (if (negative? flag) (error 'main "CVODE solver error" flag))
61
62          (if (= 1 flag)
63              ;; discrete event: re-initialize solver to integrate over the discontinuity
64              (cvode-reinit-solver solver (+ tnext trefractory) (f64vector vreset)) )
65         
66          ;;(print-results solver tnext)
67
68          (if (< tnext tf)
69              (recur (if (= 1 flag) (+ (+ tnext trefractory) dt) (+ tnext dt)) (+ 1 iout)))
70          ))
71     
72      )))
73
74
75(define (print-results solver t)
76  (let ((yy (cvode-yy solver)))
77    (printf "~A ~A ~A ~A ~A~%" 
78            t 
79            (f64vector-ref yy 0)
80            (cvode-get-last-order solver)
81            (cvode-get-num-steps solver)
82            (cvode-get-last-step solver)
83            )))
84
85     
86(main)
Note: See TracBrowser for help on using the repository browser.