source: project/release/4/sundials/trunk/tests/adex.scm @ 30926

Last change on this file since 30926 was 30926, checked in by Ivan Raikov, 6 years ago

sundials: added adaptive exponential integrate-and-fire example

File size: 2.2 KB
Line 
1
2;;
3;; Adaptive Exponential Integrate-and-Fire model.
4;;
5
6(use sundials srfi-4)
7
8
9(define  reltol 1.0e-4)
10(define  abstol 1.0e-4)
11
12
13(define  t0 0.0)
14(define  tf 1000.0)
15
16
17(define C        281.0 )
18(define gL       30.0 )
19(define EL      -70.6 )
20(define VT      -50.4 )
21(define Isyn     1000.0 )
22(define Delta    2.0 )
23(define theta   (+ VT 5.0))
24(define trefractory  0.25)
25
26(define tau_w  144.0 )
27(define a      4.0 )
28(define b      0.0805 )
29(define Vr      EL )
30
31
32
33(define (subthreshold t yy)
34  (let ((V (f64vector-ref yy 0))
35        (W (f64vector-ref yy 1)))
36    (f64vector (/ (+ (* (- gL) (- V EL)) 
37                     (* gL Delta (exp (/ (- V VT) Delta))) 
38                     (- W) 
39                     Isyn)
40                  C)
41               (/ (- (* a (- V EL)) W) tau_w)
42               )))
43
44(define (spike-detect t yy)
45  (let ((V (f64vector-ref yy 0)))
46    (s32vector (if (>= V theta) 0 1))))
47 
48   
49
50(define (main)
51 
52  (let ((yy (f64vector -65.0 0.0)))
53
54    ;; CVODE initialization
55    (let ((dt 1e-1)
56          (solver (cvode-create-solver
57                   t0 yy subthreshold 
58                   tstop: tf
59                   abstol: abstol
60                   reltol: reltol
61                   events: (s32vector 0)
62                   event-fn: spike-detect
63                   )))
64     
65      ;; In loop, call CVodeSolve, print results, and test for error.
66     
67      (let recur ((tnext (+ t0 dt)) (iout 1))
68
69        (let ((flag  (cvode-solve solver tnext)))
70
71          (if (negative? flag) (error 'main "CVODE solver error" flag))
72
73          (if (= 1 flag)
74              (let ((yy (cvode-yy solver)))
75                ;; discrete event: re-initialize solver to integrate over the discontinuity
76                (cvode-reinit-solver solver 
77                                     (+ tnext trefractory) 
78                                     (f64vector Vr (+ (f64vector-ref yy 1) b))
79                                     )))
80         
81          (print-results solver)
82
83          (if (< tnext tf)
84              (recur (if (= 1 flag) (+ (+ tnext trefractory) dt) (+ tnext dt)) (+ 1 iout)))
85          ))
86     
87      )))
88
89
90(define (print-results solver)
91  (let ((yy (cvode-yy solver)))
92    (printf "~A ~A ~A ~A ~A ~A~%" 
93            (cvode-t solver)
94            (f64vector-ref yy 0)
95            (f64vector-ref yy 1)
96            (cvode-get-last-order solver)
97            (cvode-get-num-steps solver)
98            (cvode-get-last-step solver)
99            )))
100
101     
102(main)
Note: See TracBrowser for help on using the repository browser.