# PLC module, translated from Bhalla and Iyengar 1999

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

# Based on 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']

    fkc = 1e6
    smr = 1e-3

    # constants used in subunits
    a = [400e3, 0.2e3, 400e3, 0.2e3, 20e3]                  # [=] 1/mM.s
    d = [0.13e-3, 1.049e-3, 943.4e-6, 144.5e-6, 82.34e-6]   # [=] mM
    b = [d[i]*a[i] for i in range(len(a))]                  # [=] 1/s

    k1p = a[3]
    k2p = a[4]
    k3p = b[3]
    k4p = b[4]
    k5p = a[1]
    k6p = a[4]
    k7m = b[1]
    k8p = b[4]
    k9p = a[2]
    k10p = b[0]
    k11p = a[2]
    k12p = b[0]

    k1m =  b[3]
    k2m = b[4]
    k3m = a[3]
    k4m = a[4]
    k5m = b[1]
    k6m = b[4]
    k7p = a[1]
    k8m = a[4]
    k9m = b[2]
    k10m = a[0]
    k11m = b[2]
    k12m = a[0]

    # closed loops present, but parameters are adjusted from kinetic model to meet detailed balance

    # k7p = k4m*k5p*k6p*k7m/(k4p*k5m*k6m) # error is ~zero when this commented out
    # k8m = k4m*k5p*k6p*k8p/(k4p*k5m*k6m)
    k4m = k1p*k2p*k3p*k4p/(k1m*k2m*k3m)
    k8m = k5p*k6p*k7p*k8p/(k5m*k6m*k7m)
    k10m = k9p*k5m*k10p*k1p/(k9m*k5p*k1m)
    k12m = k11p*k7p*k12p*k3m/(k11m*k7m*k3p)

    # Calculate bond graph constants from kinetic parameters
    # Note: units of kappa are fmol/s, units of K are fmol^-1
    kf = [
          k1p,
          k2p,
          k3p,
          k4p,
          k5p,
          k6p,
          k7p,
          k8p,
          k9p,
          k10p,
          k11p,
          k12p
          ]
    kr = [
        k1m,
        k2m,
        k3m,
        k4m,
        k5m,
        k6m,
        k7m,
        k8m,
        k9m,
        k10m,
        k11m,
        k12m
    ]

    k_IPR_main = [6,6]      # [=] 1/s   neglecting c_1 volume adjustment as BG already takes care of this

    k_kinetic = [k_IPR_main[0]] + kf + kf + kf + [k_IPR_main[1]] + kr + kr + kr # 3 identical and independent subunits

    # 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_SR']] + [V['V_myo']] + [1]*(num_rows-3)

    return (k_kinetic, N_cT, K_C, W)