#arnold3

FIX: CO2 O2 RuBP

#reactions
NADPH_prod:
 $pool = {1.0}NADPH
  chloroplast * NADPH_production(J, NADP, Nt)
PGA_cons:
 {1.0}PGA = {1.0}RuBP
  chloroplast * PGA_consumption(PGA, Rp, NADPH, Nt, Vcmax)
PGA_prod_Vc:
 {1.0}RuBP + {1.0}CO2 + {2.0}NADPH = {2.0}PGA
  chloroplast * carboxylation(Vc, Vj, Vp)
PGA_prod_Vo:
 {1.0}RuBP + {1.0}O2 + {2.0}NADPH = {1.5}PGA
  chloroplast * oxygenation(phi, Vc, Vj, Vp)


#parameters
Gamma = 3.81695
Jmax = 112.989573402043
Kc = 27.8535785188479
Ko = 41.9469718861558
Nt = 0.5
PPFD = 1000.0
Rd = 0.00624227154326423
Rp = 3.2
TPU = 0.227712577291006
Vcmax = 1.4749455852483
alpha = 0.19
g0 = 0.468615384615385
h = 0.6
k = 118.69
CO2 = 35.0
O2 = 21.0
RuBP = 2.0
chloroplast = 1.0


#initial values
NADPH = 0.21
PGA = 2.4


#assignment rules
!F phi = 0.21 * ((g0 + h * k / CO2 * (v_c - Rd)) * O2 / Ko) / ((g0 + (h * k - 1.6) * (v_c - Rd)) / Kc)
!F Vc = Vcmax * ((CO2 - 1.6 * (1.0 / 2.0 * ((1.6 * (Vcmax - Rd) + (k * h / CO2 * Rd + g0) * (CO2 + Kc * (1.0 + O2 / Ko)) + k * h / CO2 * Vcmax * (Gamma - CO2)) / (1.6 - k * h / CO2 * (CO2 + Kc * (1.0 + O2 / Ko))) + sqrt(pow((1.6 * (Vcmax - Rd) + (k * h / CO2 * Rd + g0) * (CO2 + Kc * (1.0 + O2 / Ko)) + k * h / CO2 * Vcmax * (Gamma - CO2)) / (1.6 - k * h / CO2 * (CO2 + Kc * (1.0 + O2 / Ko))), 2.0) + 4.0 * g0 * ((Rd * (CO2 + Kc * (1.0 + O2 / Ko)) + Vcmax * (Gamma - CO2)) / (1.6 - k * h / CO2 * (CO2 + Kc * (1.0 + O2 / Ko)))))) / (g0 + k * h / CO2 * (1.0 / 2.0) * ((1.6 * (Vcmax - Rd) + (k * h / CO2 * Rd + g0) * (CO2 + Kc * (1.0 + O2 / Ko)) + k * h / CO2 * Vcmax * (Gamma - CO2)) / (1.6 - k * h / CO2 * (CO2 + Kc * (1.0 + O2 / Ko))) + sqrt(pow((1.6 * (Vcmax - Rd) + (k * h / CO2 * Rd + g0) * (CO2 + Kc * (1.0 + O2 / Ko)) + k * h / CO2 * Vcmax * (Gamma - CO2)) / (1.6 - k * h / CO2 * (CO2 + Kc * (1.0 + O2 / Ko))), 2.0) + 4.0 * g0 * ((Rd * (CO2 + Kc * (1.0 + O2 / Ko)) + Vcmax * (Gamma - CO2)) / (1.6 - k * h / CO2 * (CO2 + Kc * (1.0 + O2 / Ko))))))))) / (CO2 - 1.6 * (1.0 / 2.0 * ((1.6 * (Vcmax - Rd) + (k * h / CO2 * Rd + g0) * (CO2 + Kc * (1.0 + O2 / Ko)) + k * h / CO2 * Vcmax * (Gamma - CO2)) / (1.6 - k * h / CO2 * (CO2 + Kc * (1.0 + O2 / Ko))) + sqrt(pow((1.6 * (Vcmax - Rd) + (k * h / CO2 * Rd + g0) * (CO2 + Kc * (1.0 + O2 / Ko)) + k * h / CO2 * Vcmax * (Gamma - CO2)) / (1.6 - k * h / CO2 * (CO2 + Kc * (1.0 + O2 / Ko))), 2.0) + 4.0 * g0 * ((Rd * (CO2 + Kc * (1.0 + O2 / Ko)) + Vcmax * (Gamma - CO2)) / (1.6 - k * h / CO2 * (CO2 + Kc * (1.0 + O2 / Ko)))))) / (g0 + k * h / CO2 * (1.0 / 2.0) * ((1.6 * (Vcmax - Rd) + (k * h / CO2 * Rd + g0) * (CO2 + Kc * (1.0 + O2 / Ko)) + k * h / CO2 * Vcmax * (Gamma - CO2)) / (1.6 - k * h / CO2 * (CO2 + Kc * (1.0 + O2 / Ko))) + sqrt(pow((1.6 * (Vcmax - Rd) + (k * h / CO2 * Rd + g0) * (CO2 + Kc * (1.0 + O2 / Ko)) + k * h / CO2 * Vcmax * (Gamma - CO2)) / (1.6 - k * h / CO2 * (CO2 + Kc * (1.0 + O2 / Ko))), 2.0) + 4.0 * g0 * ((Rd * (CO2 + Kc * (1.0 + O2 / Ko)) + Vcmax * (Gamma - CO2)) / (1.6 - k * h / CO2 * (CO2 + Kc * (1.0 + O2 / Ko)))))))) + Kc * (1.0 + O2 / Ko)))
!F Vj = J / 4.0 * ((CO2 - 1.6 * (1.0 / 2.0 * ((1.6 * (J / 4.0 - Rd) + (k * h / CO2 * Rd + g0) * (CO2 + 2.0 * Gamma) + k * h / CO2 * J / 4.0 * (Gamma - CO2)) / (1.6 - k * h / CO2 * (CO2 + 2.0 * Gamma)) + sqrt(pow((1.6 * (J / 4.0 - Rd) + (k * h / CO2 * Rd + g0) * (CO2 + 2.0 * Gamma) + k * h / CO2 * J / 4.0 * (Gamma - CO2)) / (1.6 - k * h / CO2 * (CO2 + 2.0 * Gamma)), 2.0) + 4.0 * g0 * ((Rd * (CO2 + 2.0 * Gamma) + J / 4.0 * (Gamma - CO2)) / (1.6 - k * h / CO2 * (CO2 + 2.0 * Gamma))))) / (g0 + k * h / CO2 * (1.0 / 2.0) * ((1.6 * (J / 4.0 - Rd) + (k * h / CO2 * Rd + g0) * (CO2 + 2.0 * Gamma) + k * h / CO2 * J / 4.0 * (Gamma - CO2)) / (1.6 - k * h / CO2 * (CO2 + 2.0 * Gamma)) + sqrt(pow((1.6 * (J / 4.0 - Rd) + (k * h / CO2 * Rd + g0) * (CO2 + 2.0 * Gamma) + k * h / CO2 * J / 4.0 * (Gamma - CO2)) / (1.6 - k * h / CO2 * (CO2 + 2.0 * Gamma)), 2.0) + 4.0 * g0 * ((Rd * (CO2 + 2.0 * Gamma) + J / 4.0 * (Gamma - CO2)) / (1.6 - k * h / CO2 * (CO2 + 2.0 * Gamma)))))))) / (CO2 - 1.6 * (1.0 / 2.0 * ((1.6 * (J / 4.0 - Rd) + (k * h / CO2 * Rd + g0) * (CO2 + 2.0 * Gamma) + k * h / CO2 * J / 4.0 * (Gamma - CO2)) / (1.6 - k * h / CO2 * (CO2 + 2.0 * Gamma)) + sqrt(pow((1.6 * (J / 4.0 - Rd) + (k * h / CO2 * Rd + g0) * (CO2 + 2.0 * Gamma) + k * h / CO2 * J / 4.0 * (Gamma - CO2)) / (1.6 - k * h / CO2 * (CO2 + 2.0 * Gamma)), 2.0) + 4.0 * g0 * ((Rd * (CO2 + 2.0 * Gamma) + J / 4.0 * (Gamma - CO2)) / (1.6 - k * h / CO2 * (CO2 + 2.0 * Gamma))))) / (g0 + k * h / CO2 * (1.0 / 2.0) * ((1.6 * (J / 4.0 - Rd) + (k * h / CO2 * Rd + g0) * (CO2 + 2.0 * Gamma) + k * h / CO2 * J / 4.0 * (Gamma - CO2)) / (1.6 - k * h / CO2 * (CO2 + 2.0 * Gamma)) + sqrt(pow((1.6 * (J / 4.0 - Rd) + (k * h / CO2 * Rd + g0) * (CO2 + 2.0 * Gamma) + k * h / CO2 * J / 4.0 * (Gamma - CO2)) / (1.6 - k * h / CO2 * (CO2 + 2.0 * Gamma)), 2.0) + 4.0 * g0 * ((Rd * (CO2 + 2.0 * Gamma) + J / 4.0 * (Gamma - CO2)) / (1.6 - k * h / CO2 * (CO2 + 2.0 * Gamma))))))) + 2.0 * Gamma))
!F NADP = Nt - NADPH
!F J = alpha * PPFD / sqrt(1.0 + pow(alpha * PPFD / Jmax, 2.0)) / (65.0 * 0.5)
!F v_c = (abs((abs(Vc + Vj) - abs(Vc - Vj)) / 2.0 + Vp) - abs((abs(Vc + Vj) - abs(Vc - Vj)) / 2.0 - Vp)) / 2.0
!F Vp = 3.0 * TPU * ((CO2 - 1.6 * (1.0 / 2.0 * ((1.6 * (3.0 * TPU - Rd) + (k * h / CO2 * Rd + g0) * (CO2 - Gamma) - k * h * 3.0 * TPU) / (1.6 - k * h / CO2 * (CO2 - Gamma)) + sqrt(pow((1.6 * (3.0 * TPU - Rd) + (k * h / CO2 * Rd + g0) * (CO2 - Gamma) - k * h * 3.0 * TPU) / (1.6 - k * h / CO2 * (CO2 - Gamma)), 2.0) + 4.0 * g0 * ((Rd * (CO2 - Gamma) - 3.0 * TPU * CO2) / (1.6 - k * h / CO2 * (CO2 - Gamma))))) / (g0 + k * h / CO2 * (1.0 / 2.0) * ((1.6 * (3.0 * TPU - Rd) + (k * h / CO2 * Rd + g0) * (CO2 - Gamma) - k * h * 3.0 * TPU) / (1.6 - k * h / CO2 * (CO2 - Gamma)) + sqrt(pow((1.6 * (3.0 * TPU - Rd) + (k * h / CO2 * Rd + g0) * (CO2 - Gamma) - k * h * 3.0 * TPU) / (1.6 - k * h / CO2 * (CO2 - Gamma)), 2.0) + 4.0 * g0 * ((Rd * (CO2 - Gamma) - 3.0 * TPU * CO2) / (1.6 - k * h / CO2 * (CO2 - Gamma)))))))) / (CO2 - 1.6 * (1.0 / 2.0 * ((1.6 * (3.0 * TPU - Rd) + (k * h / CO2 * Rd + g0) * (CO2 - Gamma) - k * h * 3.0 * TPU) / (1.6 - k * h / CO2 * (CO2 - Gamma)) + sqrt(pow((1.6 * (3.0 * TPU - Rd) + (k * h / CO2 * Rd + g0) * (CO2 - Gamma) - k * h * 3.0 * TPU) / (1.6 - k * h / CO2 * (CO2 - Gamma)), 2.0) + 4.0 * g0 * ((Rd * (CO2 - Gamma) - 3.0 * TPU * CO2) / (1.6 - k * h / CO2 * (CO2 - Gamma))))) / (g0 + k * h / CO2 * (1.0 / 2.0) * ((1.6 * (3.0 * TPU - Rd) + (k * h / CO2 * Rd + g0) * (CO2 - Gamma) - k * h * 3.0 * TPU) / (1.6 - k * h / CO2 * (CO2 - Gamma)) + sqrt(pow((1.6 * (3.0 * TPU - Rd) + (k * h / CO2 * Rd + g0) * (CO2 - Gamma) - k * h * 3.0 * TPU) / (1.6 - k * h / CO2 * (CO2 - Gamma)), 2.0) + 4.0 * g0 * ((Rd * (CO2 - Gamma) - 3.0 * TPU * CO2) / (1.6 - k * h / CO2 * (CO2 - Gamma))))))) - Gamma))