""" Pyramidal CAN """
""" ------------------------------------------------------------------------ """
py_CAN_inp_eqs = '''
    dv/dt = (- I_CAN - I_M - I_leak - I_K - I_Na - I_Ca - I_SynE - I_SynExt - I_SynI - I_SynHipp + r_drive*I_exc + rin*I_tonic + r*I_stim) / ((1.*ufarad*cm**-2) * (size)) + noise : volt
    I_CAN = ((gCAN) * (size)) * mCAN**2 * (v + 20.*mV) : amp
        dmCAN/dt = (mCANInf - mCAN) / mCANTau : 1
            mCANInf = alpha2 / (alpha2 + (0.0002*ms**-1)) : 1
            mCANTau = 1. / (alpha2 + (0.0002*ms**-1)) / (3.0**((36. - 22.) / 10.)) : second
            alpha2 = (0.0002*ms**-1) * (Ca_i / (5e-4*mole*metre**-3))**2 : Hz
    I_M = ((gM) * (size)) * p * (v - Ek) : amp
        dp/dt = (pInf - p) / pTau : 1
            pInf = 1. / (1. + exp(- (v + 35.*mV) / (10.*mV))) : 1
            pTau = (1000.*ms) / (3.3 * exp((v + 35.*mV) / (20.*mV)) + exp(- (v + 35.*mV) / (20.*mV))) : second
    I_leak = ((1e-5*siemens*cm**-2) * (size)) * (v - (-70.*mV)) : amp
    I_K = ((5.*msiemens*cm**-2) * (size)) * (n**4) * (v - Ek) : amp
        dn/dt = alphan * (1. - n) - betan * n : 1
            alphan = 0.032 * (mV**-1) * (5.*mV) / exprel(-(v + 40.*mV) / (5.*mV)) / ms : Hz
            betan = 0.5 * exp(- (v + 45.*mV) / (40.*mV)) / ms : Hz
    I_Na = ((50*msiemens*cm**-2) * (size)) * (m**3) * h * (v - 50.*mV) : amp
        dm/dt = alpham * (1 - m) - betam * m : 1
            alpham = 0.32 * (mV**-1) * (4.*mV) / exprel(-(v + 42.*mV) / (4.*mV)) / ms : Hz
            betam = 0.28 * (mV**-1) * (5.*mV) / exprel( (v + 15.*mV) / (5.*mV)) / ms : Hz
        dh/dt = alphah * (1 - h) - betah * h : 1
            alphah = 0.128 * exp(- (v + 38.*mV) / (18.*mV)) / ms : Hz
            betah = 4. / (1. + exp(- (v + 15.*mV) / (5.*mV))) / ms : Hz
    I_Ca = ((1e-4 * siemens*cm**-2) * (size)) * (mCaL**2) * hCaL * (v - 120.*mV) : amp
        dmCaL/dt = (alphamCaL * (1. - mCaL)) - (betamCaL * mCaL) : 1
            alphamCaL = 0.055 * (mV**-1) * (3.8*mV) / exprel(-(v + 27.*mV) / (3.8*mV)) / ms : Hz
            betamCaL = 0.94 * exp(-(v + 75.*mV) / (17.*mV)) / ms : Hz
        dhCaL/dt = (alphahCaL * (1. - hCaL)) - (betahCaL * hCaL) : 1
            alphahCaL = 0.000457 * exp(-(v + 13.*mV) / (50.*mV)) / ms : Hz
            betahCaL = 0.0065 / (exp(-(v + 15.*mV) / (28.*mV)) + 1.) / ms : Hz
        dCa_i/dt = driveChannel + ((2.4e-4*mole*metre**-3) - Ca_i) / (200.*ms) : mole*meter**-3
            driveChannel = (-(1e4) * I_Ca / (cm**2)) / (2. * (96489*coulomb*mole**-1) * (1*umetre)) : mole*meter**-3*Hz

    I_SynE = + ge * (v - 0.*mV) : amp
        dge/dt = (-ge + he) * (1. / (0.3*ms)) : siemens
        dhe/dt = - he / (5.*ms) : siemens
    I_SynExt = + ge_ext * (v - 0.*mV) : amp
        dge_ext/dt = (- ge_ext + he_ext) * (1. / (0.3*ms)) : siemens
        dhe_ext/dt = -he_ext / (5.*ms) : siemens
    I_SynHipp = + ge_hipp * (v - 0.*mV) : amp
        dge_hipp/dt = (- ge_hipp + he_hipp) * (1. / (0.3*ms)) : siemens
        dhe_hipp/dt = - he_hipp / (5.*ms) : siemens
    # I_SynI = + gi * (v - (-50.*mV)) * int(Cl>0.5) + gi * (v - (-80.*mV)) * int(Cl<=0.5): amp
    I_SynI = + gi * (v - (-80.*mV)) : amp
        dgi/dt = (- gi + hi) * (1. / (1.*ms)) : siemens
        dhi/dt = - hi / (10.*ms) : siemens

    # dCl/dt = - Cl / tau_Cl : 1

    dglu/dt = (1. - glu) / (3.*second) : 1


    noise = sigma_noise_exc * (2. * (0.1e-3*siemens) / (1.*ufarad))**.5 * randn() / sqrt(tstep) : volt/second (constant over dt)

    x_soma : metre
    y_soma : metre
    z_soma : metre
    G_sin = 1.*int(z_soma<15*mm)*int(z_soma>0*mm) : 1 # this is the mask/scaling for which neurons get the sinusoidal input
    # I_exc : amp (linked) # this is the input theta rhythm from the MS
    I_exc = inputs_drive(t) : amp
    r_drive : 1
    I_tonic = inputs_tonic(t) : amp
    # I_tonic : amp
    rin : 1
    I_stim = inputs_stim(t) : amp
    r : 1
    size : metre**2 (shared)
'''

""" Inhibitory Neuron Types """
""" ------------------------------------------------------------------------ """
inh_inp_eqs = '''
    dv/dt = ( - I_leak - I_K - I_Na - I_SynE - I_SynExt - I_SynHipp - I_SynI + r_drive*I_exc + rin*I_tonic + r*I_stim) / ((1.*ufarad*cm**-2) * (size)) + noise : volt
    I_leak = ((0.1e-3*siemens*cm**-2) * (size)) * (v - (-65.*mV)) : amp
    I_K = ((9e-3*siemens*cm**-2) * (size)) * (n**4) * (v - (-90.*mV)) : amp
        dn/dt = (alphan * (1 - n) - betan * n) / 0.2: 1
            alphan = 0.1 / exprel(-0.1*(mV**-1)*(v + 34.*mV)) / ms : Hz
            betan = 0.125 * exp( - (v + 44.*mV) / (80.*mV)) / ms : Hz
    I_Na = ((35e-3*siemens*cm**-2) * (size)) * (m**3) * h * (v - (55.*mV)) : amp
        dm/dt = (alpham * (1 - m) - betam * m) / 0.2 : 1
            alpham = 1. / exprel(-(v + 35.*mV) / (10.*mV)) / ms : Hz
            betam = 4. * exp(- (v + 60.*mV) / (18.*mV)) / ms : Hz
        dh/dt = (alphah * (1 - h) - betah * h) / 0.2 : 1
            alphah = 0.07 * exp(- (v + 58.*mV) / (20.*mV)) / ms : Hz
            betah = 1. / (exp((- 0.1 * (mV**-1)) * (v + 28.*mV)) + 1.) / ms : Hz
    I_SynE = + ge * (v - 0.*mV) : amp
        dge/dt = (-ge+he) * (1. / (0.3*ms)) : siemens
        dhe/dt = -he/(5.*ms) : siemens
    I_SynExt = + ge_ext * (v - 0.*mV) : amp
        dge_ext/dt = (-ge_ext+he_ext) * (1. / (0.3*ms)) : siemens
        dhe_ext/dt = -he_ext/(5.*ms) : siemens
    I_SynHipp = + ge_hipp * (v - 0.*mV) : amp
        dge_hipp/dt = (-ge_hipp+he_hipp) * (1. / (0.3*ms)) : siemens
        dhe_hipp/dt = -he_hipp/(5.*ms) : siemens
    I_SynI = + gi * (v - (-80.*mV)) : amp
        dgi/dt = (-gi+hi) * (1. / (1.*ms)) : siemens
        dhi/dt = -hi/(10.*ms) : siemens


    noise = sigma_noise_inh * (2. * (0.1e-3*siemens) / (1*ufarad))**.5 * randn() / sqrt(tstep) : volt/second (constant over dt)


    x_soma : metre
    y_soma : metre
    z_soma : metre
    G_sin = 1.*int(z_soma<15*mm)*int(z_soma>0*mm) : 1  # this is the mask/scaling for which neurons get the sinusoidal input
    # I_exc : amp (linked)    # same as in the pyCAN group, excitatory input from MS
    I_exc = inputs_drive(t) : amp
    r_drive : 1
    I_tonic = inputs_tonic(t) : amp
    # I_tonic : amp
    rin : 1
    I_stim = inputs_stim(t) : amp
    r : 1
    size : metre**2 (shared)
'''

# Spike and reset
reset_eqs = '''
    glu = glu - 0.
    Cl = Cl + 0.2
'''

# Spike and reset
reset_eqs_noCl = '''
    glu = glu - 0.
'''


# Kuramoto oscillators
kuramoto_eqs_stim = '''
    dTheta/dt = ((omega + (kN * PIF) - G_in*(X+stim_MS(t))*sin(Theta + offset)) * second**-1) : 1
    PIF = .5 * (sin(ThetaPreInput - Theta)) : 1

    ThetaPreInput : 1
    omega : 1 (constant)
    kN : 1 (shared)         # k/N ratio, affects sync.
    G_in : 1 (shared)       # input gain, affects the phase reset aggressiveness
    offset : 1 (shared)     # range [0, 2*pi], controls phase reset curve
    X : 1 (linked)          # this is linked to the firing rates
'''

# synapses
syn_kuramoto_eqs = '''
    ThetaPreInput_post = Theta_pre
'''

# Order parameter group calculation equations
pop_avg_eqs = '''
    coherence = sqrt(x**2 + y**2) : 1
    phase = arctan(y/x) + int(x<0 and y>0)*pi - int(x<0 and y<0)*pi: 1

    # Rhythms
    rhythm_default = coherence * cos(phase) : 1
    rhythm_positive = coherence * (cos(phase)+1)/2 : 1
    rhythm_abs = abs(rhythm_default) : 1
    rhythm_rect = rhythm_positive : 1
    rhythm_zero = 0.*rhythm : amp   # for debugging

    # Output selection
    rhythm = G_out*rhythm_rect : amp

    x : 1
    y : 1
    G_out : amp                 # output rhythm gain
'''

syn_avg_eqs = '''
    x_post = cos(Theta_pre)/N_incoming : 1 (summed)
    y_post = sin(Theta_pre)/N_incoming : 1 (summed)
'''


# Vm avgs
eq_record_neurons = '''
    sum_v : volt
'''

eq_record_synapses = '''
    sum_v_post = v_pre/N_incoming : volt (summed)
'''


# LFPs estimation
eq_record_LFP_neurons = '''
    sum_I_SynE : amp
    sum_I_SynI : amp
'''

eq_record_LFP_synapses = '''
    sum_I_SynE_post = I_SynE_pre : amp (summed)
    sum_I_SynI_post = I_SynI_pre : amp (summed)
'''