# PLB module
# all rxn directions align with kinetic Saucerman model (Km unchanged,
# aside from dimensional scaling)

# 5 AUG PLB separated into PLB and Inhib1 submodules
# (parameter error was large with both submodules together)
# Use [Ip_S1] where [Ip] = [Ip_s1] + [Ip_S2]

# 25 Aug
#     return (k_kinetic, N_cT, K_C, W) kinetic parameters, constraints, and vector of volumes in each
# compartment (pL) (1 if gating variable, or in element corresponding to
# kappa)

import numpy as np 

def kinetic_parameters(M, include_all_reactions, dims, V):
    # Set the kinetic rate constants.
    # all reactions are reversible. no closed loops.
    # cAMP binds to R subunit one at a time
    
    num_cols = dims['num_cols']
    num_rows = dims['num_rows']

    # CONVERT TO fM 
    bigNum = 1e3
    fastKineticConstant = bigNum
    smallReverse = fastKineticConstant/(pow(bigNum,2))
    
    cPP2A = 0.224 # uM (Bhalla Iyengar a on PMR)    but not sure what cell type this is

    # reaction IDs
#         PLBph1, PLBph2, 
#         PLBd1 (reverse), PLBd2 (reverse)
#         inh

    k_pka_plb = 54 #: per_sec 54
    Km_pka_plb = 21 #: uM 21
    k_pp1_plb = 8.5 # per_sec 8.5
    Km_pp1_plb = 7 # uM 7
#     k_pka_i1 = 60 #: per_sec 60
#     Km_pka_i1 = 1 #: uM 1
#     Vmax_pp2a_i1 = 14 #: uM_per_sec 14
#     Km_pp2a_i1 = 1 #: uM 1
    Ki_inhib1 = 1e-3 #: uM 1e-3

    vKm = [Km_pka_plb, Km_pp1_plb 
#         Km_pka_i1, Km_pp2a_i1
        ]
    vkcat = [k_pka_plb, 
        k_pp1_plb 
#         k_pka_i1, 
#         Vmax_pp2a_i1/cPP2A
        ]

    N = len(vKm) 
    k1m = [] #np.zeros(N)
    k1p = [] #np.zeros(N)
    k2m = [] #np.zeros(N)
    k2p = [] #np.zeros(N)

    for i in range(N): 
        k2p.append(vkcat[i])
        k1m.append(fastKineticConstant) # 1/s
        k1p.append((k1m[i] + k2p[i]) / vKm[i])
        k2m.append(k1p[i]*k2p[i]/k1m[i]) # detailed balance
    
    km_inh = fastKineticConstant
    kp_inh = km_inh / Ki_inhib1 # type 2
        
    if include_all_reactions:
        k_kinetic = k1p + k2p + [kp_inh] + k1m + k2m + [km_inh]
    else:
        irr = range(3)
        k_kinetic = [
            k1p[irr], k2p[irr], k1m[irr], k2m[irr]
            ]
    

    # CONSTRAINTS
    N_cT = [] 
#     N_cT = [] 
# #     # Reaction i: [PKA:PKI] = [C][PKI] at SS   big error. Not isolated reaction
#     # repeat for type 1 and 2
#     if False:
#         N_cT[0][num_cols + 1] = 1
#         N_cT[0][num_cols + 5] = 1
#         N_cT[0][num_cols + 8) = -1
#     
#     # Gibbs free energy of L + R binding                            **MED_ERROR**
#     if False:
#         N_cT[2][num_cols + 3] = 1   # ARC
#         N_cT[2][num_cols) = -1  # cAMP
#         N_cT[2][num_cols + 2] = -1  # RC
#         G_0_bind = -45.1872 # kJ/mol
#         R = 8.314
#         T = 310
#         K_bind = np.exp(G_0_bind/(R*T))*10^6
#     
    
    K_C = []
    
    # volume vector
    W = list(np.append([1] * num_cols, [V['V_myo']] * num_rows))

    return (k_kinetic, [N_cT], K_C, W)