import numpy as np
import matplotlib.pyplot as plt

#place the paramters here:

...

#

T = 50
dt = 0.01
time = np.arange(0, T + dt, dt)

I_ext = np.zeros(len(time))
I_ext[1000:2000] = 10

def alpha_n(V):
    return 0.01 * (V + 55) / (1 - np.exp(-0.1 * (V + 55)))

def beta_n(V):
    return 0.125 * np.exp(-0.0125 * (V + 65))

def alpha_m(V):
    return 0.1 * (V + 40) / (1 - np.exp(-0.1 * (V + 40)))

def beta_m(V):
    return 4.0 * np.exp(-0.0556 * (V + 65))

def alpha_h(V):
    return 0.07 * np.exp(-0.05 * (V + 65))

def beta_h(V):
    return 1 / (1 + np.exp(-0.1 * (V + 35)))

V = -65.0
n = alpha_n(V) / (alpha_n(V) + beta_n(V))
m = alpha_m(V) / (alpha_m(V) + beta_m(V))
h = alpha_h(V) / (alpha_h(V) + beta_h(V))

V_values = []
n_values = []
m_values = []
h_values = []

for t in time:
    I_Na = (m**3) * g_Na * h * (V - E_Na)
    I_K = (n**4) * g_K * (V - E_K)
    I_L = g_L * (V - E_L)

    dV = (I_ext[int(t / dt)] - I_Na - I_K - I_L) / C_m
    V += dV * dt

    dn = (alpha_n(V) * (1 - n) - beta_n(V) * n) * dt
    dm = (alpha_m(V) * (1 - m) - beta_m(V) * m) * dt
    dh = (alpha_h(V) * (1 - h) - beta_h(V) * h) * dt
    n += dn
    m += dm
    h += dh

    V_values.append(V)
    n_values.append(n)
    m_values.append(m)
    h_values.append(h)

V_values = np.array(V_values)
n_values = np.array(n_values)
m_values = np.array(m_values)
h_values = np.array(h_values)

plt.figure(figsize=(12, 8))

plt.subplot(2, 1, 1)
plt.plot(time, V_values, label='Membrane Potential (V)')
plt.title('Hodgkin-Huxley Neuron Model')
plt.ylabel('Membrane Potential (mV)')
plt.legend()

plt.subplot(2, 1, 2)
plt.plot(time, I_ext, label='External Current (I_ext)', color='r')
plt.xlabel('Time (ms)')
plt.ylabel('Current (uA/cm^2)')
plt.legend()

plt.tight_layout()
plt.show()