source: project/release/4/sundials/trunk/tests/izhfs.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.5 KB
Line 
1;;
2;; Example of Izhikevich fast spiking neuron
3;;
4
5
6(use sundials srfi-4)
7
8
9(define  reltol 1.0e-6)
10(define  abstol 1.0e-6)
11
12
13(define  t0 0.0)
14(define  tf 800.0)
15
16(define k     1.0)
17(define Vpeak 25.0)
18(define Vt    -55.0)
19(define Vr    -40.0)
20(define Vb    -55.0)
21(define Cm    20.0)
22(define Isyn  0.0)
23
24;; State initial values
25(define  V      -65.0)
26(define  tspike   0.0)
27(define  spike     #f)
28
29(define  FS_a 0.2)
30(define  FS_b 0.025)
31(define  FS_c -45.0)
32(define  FS_U (* FS_b V))
33
34(define Iext 100.)
35
36(define (UU v)  (if (< v Vb) 0. (* FS_b (- v Vb) (- v Vb) (- v Vb))))
37
38(define (subthreshold t yy)
39  (let ((v (f64vector-ref yy 0))
40        (u (f64vector-ref yy 1)))
41    (let (
42          (dv (/ (+ (* k (- v Vr) (- v Vt)) (- u) Iext) Cm))
43          (du (* FS_a (- (UU v) u)))
44          )
45      (f64vector dv du)
46      )))
47
48
49(define (vreset t h yy)
50  (let ((v (f64vector-ref yy 0))
51        (u (f64vector-ref yy 1)))
52  (let ((v1 FS_c)
53        (u1 (+ u (* h FS_a (- (* FS_b (- FS_c (* FS_b (- v Vb) (- v Vb) (- v Vb)) )) u)))))
54
55    (f64vector v1 u1)
56    )))
57   
58
59
60(define (spike-detect t yy)
61  (let ((v (f64vector-ref yy 0)))
62    (s32vector (floor (- v Vpeak)))))
63             
64
65(define (main)
66 
67  (let ((yy (f64vector V FS_U)))
68
69    ;; CVODE initialization
70    (let ((h 1e-1)
71          (solver (cvode-create-solver
72                   t0 yy subthreshold 
73                   tstop: tf
74                   abstol: abstol
75                   reltol: reltol
76                   events: (s32vector 0)
77                   event-fn: spike-detect
78                   )))
79     
80      ;; In loop, call CVodeSolve, print results, and test for error.
81     
82      (let recur ((tnext (+ t0 h)) (iout 1))
83
84        (let ((flag  (cvode-solve solver tnext)))
85
86          (if (negative? flag) (error 'main "CVODE solver error" flag))
87
88          (if (= 1 flag)
89              (begin
90                ;; discrete event: re-initialize solver to integrate over the discontinuity
91                (print-results solver)
92                (cvode-reinit-solver solver (+ tnext h) (vreset (cvode-t solver) h (cvode-yy solver)) )
93                (recur tnext (+ iout 1)))
94
95              (if (< (cvode-t solver) tnext) 
96                  (recur tnext (+ 1 iout))
97                  (begin
98
99                    (print-results solver)
100                    (if (< tnext tf)
101                        (recur (if (= 1 flag) (+ (+ tnext h) h) (+ tnext h)) (+ 1 iout)))
102                    )))
103          ))
104      ))
105  )
106
107(define (print-results solver)
108  (let ((yy (cvode-yy solver))
109        (t (cvode-t solver)))
110    (let ((spike (= (s32vector-ref (spike-detect t yy) 0) 0)))
111    (printf "~A ~A ~A ~A ~A ~A ~A ~A~%" 
112            t (or (and spike 1) 0) (or (and spike t) 0.0)
113            (f64vector-ref yy 0)
114            (f64vector-ref yy 1)
115            (cvode-get-last-order solver)
116            (cvode-get-num-steps solver)
117            (cvode-get-last-step solver)
118            ))
119    ))
120
121     
122(main)
Note: See TracBrowser for help on using the repository browser.