/
stims.f95
160 lines (143 loc) · 6 KB
/
stims.f95
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
! -*- f95 -*-
! (c) 2016 - Ilya Prokin - isprokin@gmail.com - https://sites.google.com/site/ilyaprokin
! INRIA Rhone-Alpes
! STDP model : file defines functions describing periodical stimulation,
! tabulating functions needed for stimulation and helpers defining
! singularity time points for integration
module stims
use pars_mod
use general_math
implicit none
type stims_tabs_type
real*8, allocatable :: Glu_tab(:), Iact_tab(:)
integer :: Glu_tab_li, Iact_tab_li
real*8 :: Glu_tab_step, Iact_tab_step
end type stims_tabs_type
type(stims_tabs_type), save :: stims_tabs
contains
subroutine stims_first_tcrit_set(pars, stims_tabs, stims_first_tcrit)
use qsort_c_module
implicit none
type(pars_type), intent(in) :: pars
type(stims_tabs_type), intent(in) :: stims_tabs
real*8, intent(out) :: stims_first_tcrit(5)
! Glu
stims_first_tcrit(1)=pars%integration%t_start + pars%stimulation%tpost-pars%stimulation%tsdt-pars%stimulation%dt_stim
stims_first_tcrit(2)=stims_first_tcrit(1)+stims_tabs%Glu_tab_li*pars%stimulation%tables_step
! Iact
stims_first_tcrit(3)=pars%integration%t_start + pars%stimulation%tpost-2.0*pars%stimulation%tsdt
stims_first_tcrit(4) = stims_first_tcrit(3) + pars%stimulation%tsdt
stims_first_tcrit(5)=stims_first_tcrit(3)+stims_tabs%Iact_tab_li*pars%stimulation%tables_step
call QsortC(stims_first_tcrit)
end subroutine stims_first_tcrit_set
subroutine tcrit_per_stim_set_helper(t, tcritlast, period, ITASK, RWORK, LRW, tcrit, ncrit, icrit)
! to correctly set TCRIT for the solvers from ODEPACK
! pediod is a period of periodic stimulation = 1d0/pars%stimulation%Freq
! icrit should be declared in main program and set icrit=1 if all values in tcrit should be considered
! call in a loop before a solver. Ex. with LSODA:
! ! initializations...
! do i=2,n
! call tcrit_set_helper(t(i), tcrit, ncrit, ITASK, RWORK, LRW, iCRIT)
! call DLSODA(FEX,NEQ,ynew,t_in,t(i),ITOL,RTOL,ATOL,ITASK,ISTATE, IOPT,RWORK,LRW,IWORK,LIW,JDUM,JT)
! y(:,i)=ynew
! end do
implicit none
integer, intent(in) :: LRW, ncrit
real*8, intent(in) :: t, tcritlast, period
real*8, intent(inout) :: RWORK(LRW)
integer, intent(inout) :: ITASK
integer, intent(inout) :: icrit
real*8, intent(inout) :: tcrit(ncrit)
do
if (t>tcritlast) then
ITASK=1
exit
end if
if (icrit<=ncrit) then
if (tcrit(icrit)>=t) then
ITASK=4
RWORK(1)=tcrit(icrit)
exit
else
icrit=icrit+1
end if
else
icrit = 1
tcrit = tcrit + period
end if
end do
end subroutine tcrit_per_stim_set_helper
subroutine stims_tables_make(pars, stims_tabs)
implicit none
type(pars_type), intent(in) :: pars
type(stims_tabs_type), intent(out) :: stims_tabs
stims_tabs%Glu_tab_step = pars%stimulation%tables_step
call Glu_pulse(stims_tabs%Glu_tab_step, pars, stims_tabs%Glu_tab)
stims_tabs%Glu_tab_li = size(stims_tabs%Glu_tab)
stims_tabs%Iact_tab_step = pars%stimulation%tables_step
call Iact_pulse(stims_tabs%Iact_tab_step, pars, stims_tabs%Iact_tab)
stims_tabs%Iact_tab_li = size(stims_tabs%Iact_tab)
end subroutine stims_tables_make
subroutine stims_tables_clean(stims_tabs)
type(stims_tabs_type) :: stims_tabs
deallocate(stims_tabs%Glu_tab)
deallocate(stims_tabs%Iact_tab)
end subroutine stims_tables_clean
subroutine Glu_pulse(table_step, pars, Tabl)
implicit none
real*8, intent(in) :: table_step
type(pars_type), intent(in) :: pars
integer :: i
integer :: li
real*8, allocatable, intent(out) :: Tabl(:)
real*8 :: t
li = pars%Glu_release%tauGlu*30/table_step+1
allocate(Tabl(li))
do i = 1,li
t=(i-1)*table_step
Tabl(i) = pars%Glu_release%Glumax*exp(-t/pars%Glu_release%tauGlu)
end do
end subroutine Glu_pulse
subroutine Iact_pulse(table_step, pars, Tabl)
implicit none
real*8, intent(in) :: table_step
type(pars_type), intent(in) :: pars
integer :: i
integer :: li
real*8, allocatable, intent(out) :: Tabl(:)
real*8 :: t
li = int(pars%action%APdur/table_step)+1
allocate(Tabl(li))
do i = 1,li
t=(i-1)*table_step - pars%stimulation%tsdt
Tabl(i) = -pars%action%DPmax -pars%action%APmax * exp(-t/pars%action%tausbAP)*heaviside(t)
end do
end subroutine Iact_pulse
! functions for periodical stims
real*8 function Glu_func(t, stims_tabs, pars)
implicit none
type(pars_type) :: pars
type(stims_tabs_type) :: stims_tabs
real*8 :: t, t0, tp
t0=pars%integration%t_start + pars%stimulation%tpost-pars%stimulation%tsdt-pars%stimulation%dt_stim
tp=t-t0
if (tp>0) then
Glu_func = kernel_repeat(tp, pars%stimulation%Freq, int(pars%stimulation%num_stim), stims_tabs%Glu_tab, stims_tabs%Glu_tab_li, stims_tabs%Glu_tab_step)
else
Glu_func = 0.0
end if
end function Glu_func
real*8 function Iact_func(t, stims_tabs, pars)
implicit none
type(pars_type) :: pars
type(stims_tabs_type) :: stims_tabs
real*8 :: t,t0, tp
t0=pars%integration%t_start + pars%stimulation%tpost-2.0*pars%stimulation%tsdt
tp=t-t0
if (tp>0) then
Iact_func = kernel_repeat(tp, pars%stimulation%Freq, int(pars%stimulation%num_stim), stims_tabs%Iact_tab, stims_tabs%Iact_tab_li, stims_tabs%Iact_tab_step)
else
Iact_func = 0.0
end if
end function Iact_func
end module stims