# fast Na module

# Return kinetic parameters, constraints, and vector of volumes in each
# compartment (pL) (1 if gating variable, or in element corresponding to
# kappa)

# Translated from Pan 2018 cardiac AP 

import numpy as np

def kinetic_parameters(M, include_type2_reactions, dims, V):
    # Set the kinetic rate constants

    num_cols = dims['num_cols']
    num_rows = dims['num_rows']

    # load gate transition parameters
    params_m = [12.516449859118975,	0.4954094890610479,	0.0799997602584941,	-2.4283951182982153]
    params_h = [3.353932648035526E-7,	-4.189185289647097,	0.7999028222966151,	1.2995102414503217]
    params_j = [1.3078733519642442E-7,	-4.038149019814894,0.42275822436453375,	1.428087934426136]

    alpha_m0_bg = params_m[0]*1e3 # unit s^-1
    beta_m0_bg = params_m[2]*1e3 # unit s^-1
    alpha_h0_bg = params_h[0]*1e3 # unit s^-1
    beta_h0_bg = params_h[2]*1e3 # unit s^-1
    alpha_j0_bg = params_j[0]*1e3 # unit s^-1
    beta_j0_bg = params_j[2]*1e3 # unit s^-1

    # constants are stored in V
    F = V['F']
    R = V['R']
    T = V['T']
    N_A = V['N_A']
    A_cap = V['A_cap']
    cNao = V['cNao'] 
    cNai = V['cNai'] 
    G_Na = 16*A_cap # Unit mA/V
    E_Na = R*T/F*np.log(cNao/cNai)
    E_Na_norm = F*E_Na/(R*T)
    G_GHK = 2*G_Na*(1-np.exp(E_Na_norm))/(cNai - cNao*np.exp(E_Na_norm))*R*T/F # Unit mA/mM
    P_Na = G_GHK/F * 1e12 # Unit pL/s . G_GHK [=] Amp/(mol/s) 
    x_Na_channel = 122720/N_A*1e15 # Pan cell value. unit fmol
    x_Na_channel = V['numChannels']/N_A*1e15 # unit fmol

    # Calculate bond graph constants from kinetic parameters
    # Note: units of kappa are fmol/s, units of K are fmol^-1
    kf_Na = [P_Na/x_Na_channel,  # R_GHK
        3*alpha_m0_bg,  # Rm000
        3*alpha_m0_bg,  # Rm001
        3*alpha_m0_bg,  # Rm010
        3*alpha_m0_bg,  # Rm011
        2*alpha_m0_bg,  # Rm100
        2*alpha_m0_bg,  # Rm101
        2*alpha_m0_bg,  # Rm110
        2*alpha_m0_bg,  # Rm111
        alpha_m0_bg,  # Rm200
        alpha_m0_bg,  # Rm201
        alpha_m0_bg,  # Rm210
        alpha_m0_bg,  # Rm211
        alpha_h0_bg,  # Rh000
        alpha_h0_bg,  # Rh001
        alpha_h0_bg,  # Rh100
        alpha_h0_bg,  # Rh101
        alpha_h0_bg,  # Rh200
        alpha_h0_bg,  # Rh201
        alpha_h0_bg,  # Rh300
        alpha_h0_bg,  # Rh301
        alpha_j0_bg,  # Rj000
        alpha_j0_bg,  # Rj010
        alpha_j0_bg,  # Rj100
        alpha_j0_bg,  # Rj110
        alpha_j0_bg,  # Rj200
        alpha_j0_bg,  # Rj210
        alpha_j0_bg,  # Rj300
        alpha_j0_bg] # Rj310

    kr_Na = [P_Na/x_Na_channel,  # R_GHK
        beta_m0_bg,  # Rm000
        beta_m0_bg,  # Rm001
        beta_m0_bg,  # Rm010
        beta_m0_bg,  # Rm011
        2*beta_m0_bg,  # Rm100
        2*beta_m0_bg,  # Rm101
        2*beta_m0_bg,  # Rm110
        2*beta_m0_bg,  # Rm111
        3*beta_m0_bg,  # Rm202
        3*beta_m0_bg,  # Rm201
        3*beta_m0_bg,  # Rm210
        3*beta_m0_bg,  # Rm211
        beta_h0_bg,  # Rh000
        beta_h0_bg,  # Rh001
        beta_h0_bg,  # Rh100
        beta_h0_bg,  # Rh101
        beta_h0_bg,  # Rh200
        beta_h0_bg,  # Rh201
        beta_h0_bg,  # Rh300
        beta_h0_bg,  # Rh301
        beta_j0_bg,  # Rj000
        beta_j0_bg,  # Rj010
        beta_j0_bg,  # Rj100
        beta_j0_bg,  # Rj110
        beta_j0_bg,  # Rj200
        beta_j0_bg,  # Rj210
        beta_j0_bg,  # Rj300
        beta_j0_bg] # Rj310

    k_kinetic = kf_Na + kr_Na

    # CONSTRAINTS
    N_cT = []
    K_C = []

    # volume vector
    # W = list(np.append([1] * num_cols, [V['V_myo']] * num_rows))
    W = [1] * num_cols + [V['V_myo'], V['V_o']] + [1] * (num_rows-2)

    return (k_kinetic, N_cT, K_C, W)