187605/general_math.f95 at master · ModelDBRepository/187605 · GitHub

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
!    -*- f95 -*-
! (c) 2016 - author information redacted by request 2026-08-07
! INRIA Rhone-Alpes
! STDP model : math functions
module general_math
    implicit none

    contains

    real*8 function heaviside(x)
        implicit none
        real*8 :: x
        if (x == 0) then
            heaviside = 0.5
        elseif (x < 0) then
            heaviside = 0
        else
            heaviside = 1.0
        endif
    end function heaviside

    real*8 function Rect(x, xlo, xhi)
        implicit none
        real*8 :: x, xlo, xhi
        if (x < xlo) then
            Rect = 0
        elseif (x > xhi) then
            Rect = 0
        else
            Rect = 1.0
        endif
    end function Rect

    real*8 function Hill(x,k,n)
        implicit none
        real*8 :: x,k
        integer :: n
        real*8 :: xn
        xn=x**n
        Hill = xn/(xn+k**n)
    end function Hill

    real*8 function lin_interpTable(Table,li,i_point, default)
        implicit none
        real*8 :: i_point
        integer :: ip, li
        real*8 :: Table(li), default
        ip=floor(i_point)
        if (ip<1) then
            lin_interpTable = default
        elseif (ip<li) then
            lin_interpTable = Table(ip)+(Table(ip+1)-Table(ip))*(i_point-ip)
        else
            lin_interpTable = default
        endif
    end function lin_interpTable

    real*8 function lin_interpTable_brds(Table,li,i_point)
        implicit none
        real*8 :: i_point
        integer :: ip, li
        real*8 :: Table(li)
        ip=floor(i_point)
        if (ip<1) then
            lin_interpTable_brds = Table(1)
        elseif (ip<li) then
            lin_interpTable_brds = Table(ip)+(Table(ip+1)-Table(ip))*(i_point-ip)
        else
            lin_interpTable_brds = Table(li)
        endif
    end function lin_interpTable_brds

    subroutine search_xa_inds_around_x(xa,n,x, klo, khi)
        implicit none
        integer, intent(in) :: n
        real*8, intent(in) :: xa(n), x
        integer, intent(out) :: khi,klo
        integer :: k
        klo=1
        khi=n
        do
            if (khi-klo <= 1) exit
            k=(khi+klo)/2
            if(xa(k) > x) then
                khi=k
            else
                klo=k
            endif
        end do
    end subroutine search_xa_inds_around_x

    real*8 function lin_interpTable_xy(x,y,n, x_new)
        implicit none
        integer :: n, klo, khi
        real*8 :: x(n), y(n), x_new
        if (x_new<=x(1)) then
            lin_interpTable_xy = y(1)
        elseif (x_new>=x(n)) then
            lin_interpTable_xy = y(n)
        else
            call search_xa_inds_around_x(x,n,x_new, klo, khi)
            if (y(khi)/=y(klo)) then
                lin_interpTable_xy = y(klo)+(y(khi)-y(klo))*(x_new-x(klo))/(x(khi)-x(klo))
            else
                lin_interpTable_xy = y(klo)
            end if
        endif
    end function lin_interpTable_xy

    real*8 function kernel_repeat(t, Freq, n, kernel, l_ker, step_ker)
        !t -grid starts from 0.0 (t=t'-t0)
        implicit none
        real*8 :: t, Freq
        integer :: n, l_ker
        real*8 :: kernel(l_ker), step_ker
        integer :: ic, idepth, k
        ic = floor(t*Freq)
        idepth = floor(l_ker*step_ker*Freq)
        kernel_repeat = 0
        if (t>l_ker*step_ker+(n-1)/Freq) then
            return
        end if
        do k=0,idepth
            kernel_repeat = kernel_repeat + lin_interpTable(kernel, l_ker, (t-(k+ic)/Freq)/step_ker, 0d0)
        end do
    end function kernel_repeat

end module general_math