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

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

sundials: small tweaks to example model parameters

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 500.0)
15
16
17(define C       200.0 )
18(define gL      10.0 )
19(define EL      -58.0 )
20(define VT      -50.0 )
21(define Delta   2.0 )
22(define theta   0.0 )
23(define trefractory  0.25)
24
25(define a       2.0 )
26(define tau_w  120.0 )
27(define b       100.0 )
28(define Vr     -46.0 )
29
30(define Isyn     210.0 )
31
32(define (subthreshold t yy)
33  (let ((V (f64vector-ref yy 0))
34        (W (f64vector-ref yy 1)))
35    (f64vector (/ (+ (* (- gL ) (- V EL))
36                     (* gL Delta (exp (/ (- V VT) Delta)))
37                     (- W) Isyn) 
38                  C)
39               (/ (- (* a (- V EL)) W) tau_w)
40               )))
41
42(define (spike-detect t yy)
43  (let ((V (f64vector-ref yy 0)))
44    (s32vector (if (>= V theta) 0 1))))
45 
46   
47
48(define (main)
49 
50  (let ((yy (f64vector -70.0 0.0)))
51
52    ;; CVODE initialization
53    (let ((dt 1e-1)
54          (solver (cvode-create-solver
55                   t0 yy subthreshold 
56                   tstop: tf
57                   abstol: abstol
58                   reltol: reltol
59                   events: (s32vector 0)
60                   event-fn: spike-detect
61                   )))
62     
63      ;; In loop, call CVodeSolve, print results, and test for error.
64     
65      (let recur ((tnext (+ t0 dt)) (iout 1))
66
67        (let ((flag  (cvode-solve solver tnext)))
68
69          (if (negative? flag) (error 'main "CVODE solver error" flag))
70
71          (if (= 1 flag)
72              (let ((yy (cvode-yy solver)))
73                ;; discrete event: re-initialize solver to integrate over the discontinuity
74                (cvode-reinit-solver solver 
75                                     (+ tnext trefractory) 
76                                     (f64vector Vr (+ (f64vector-ref yy 1) b))
77                                     )))
78         
79          (print-results solver)
80
81          (if (< tnext tf)
82              (recur (if (= 1 flag) (+ (+ tnext trefractory) dt) (+ tnext dt)) (+ 1 iout)))
83          ))
84     
85      )))
86
87
88(define (print-results solver)
89  (let ((yy (cvode-yy solver)))
90    (printf "~A ~A ~A ~A ~A ~A~%" 
91            (cvode-t solver)
92            (f64vector-ref yy 0)
93            (f64vector-ref yy 1)
94            (cvode-get-last-order solver)
95            (cvode-get-num-steps solver)
96            (cvode-get-last-step solver)
97            )))
98
99     
100(main)
Note: See TracBrowser for help on using the repository browser.