LIF_neuron_sinusoidal_input

import os
import ipdb # ipdb.set_trace()
import numpy as np
import matplotlib.pyplot as plt

InputCurrent_maxvalue = 8e-9
C = 1e-9 #capacitance
Vth = 0.3
Vreset = 0
tstep = 1e-3
Ileak = 1e-9


times = np.arange(0,5,tstep)
#InputCurrent = 0.5*InputCurrent_maxvalue*(0.5+0.5*np.sin(2*np.pi*times))
InputCurrent = InputCurrent_maxvalue*(np.sin(2*np.pi*times))
InputCurrent = np.where(InputCurrent > 0, InputCurrent, 0)
vmembrane = np.zeros(len(times))
vmembrane[0] = 0

vout = np.zeros(len(times))
vout[0] = 0
for i, t in enumerate(times):
    # I = C_M dV/dt + Il
    # dV/dt = I/C
    dV = InputCurrent[i]/C * tstep

    # Ileak = C * dVl/dT
    # dV_l/dT = Ileak/C
    dV -= Ileak/C * tstep

    vmembrane[i] = vmembrane[i-1] + dV
    vout[i] = 0
    if vmembrane[i] > Vth:
        vmembrane[i] = 0
        vout[i] = 1
    elif vmembrane[i] < 0:
        vmembrane[i] = 0


fig, axs = plt.subplots(3)
fig.suptitle('Integrate and Fire (IF) Spiking Neuron')

axs[0].plot(times,InputCurrent,label="Input Current")
axs[0].set(xlabel="Time (s)")
axs[0].set(ylabel="Input Current (A)")
axs[0].legend()

axs[1].plot(times,vmembrane,label="Voltage Membrane")
axs[1].plot(times,Vth*np.ones(len(times)),'--y',label="threshold")
axs[1].set(xlabel="Time (s)")
axs[1].set(ylabel="Membrane Voltage (V)")
axs[1].legend()

axs[2].plot(times,vout,'r',label="Output Voltage (V)" )
axs[2].set(xlabel="Time (s)")
axs[2].set(ylabel="Output Voltage (V)")
axs[2].legend()

plt.plot
plt.show()