电解质溶液中微粒数量关系

; Li Ying created this module at Shaanxi Normal University in 2024. If you mention this model in a publication, we ask that you include the citations below.
;Ying, L. (2024). Particle Quantity Relationship in Electrolyte Solution. School of Chemistry and Chemical Engineering, Shaanxi Normal University.
extensions [py]

breed [mode6_h2c2o4_t  mode6_h2c2o4]
breed [mode6_hc2o4_t  mode6_hc2o4]
breed [mode6_c2o4_t  mode6_c2o4]
breed [mode6_h_t  mode6_h]
breed [mode6_oh_t  mode6_oh]
breed [mode6_h3po4_t  mode6_h3po4]
breed [mode6_h2po4_t  mode6_h2po4]
breed [mode6_hpo4_t  mode6_hpo4]
breed [mode6_po4_t  mode6_po4]
breed [mode6_na2co3_t  mode6_na2co3]
breed [mode6_na_t  mode6_na]
breed [mode6_co3_t  mode6_co3]
breed [mode6_hco3_t  mode6_hco3]
breed [mode6_h2co3_t  mode6_h2co3]
breed [mode6_so3_t  mode6_so3]
breed [mode6_hso3_t  mode6_hso3]
breed [mode6_h2so3_t  mode6_h2so3]
breed [mode6_cl_t  mode6_cl]
breed [mode6_nh3h2o_t  mode6_nh3h2o]
breed [mode6_nh4_t  mode6_nh4]
breed [mode6_hac_t  mode6_hac]
breed [mode6_ac_t  mode6_ac]
breed [mode6_hcn_t  mode6_hcn]
breed [mode6_cn_t  mode6_cn]

globals[
  a_
  b_
  c_
  d_
  e_
  f_
  g_
  h_
  k_
  l_
  m_

  H2C2O4
  HC2O4
  C2O4
  H_1
  OH_1

  H3PO4
  H2PO4
  HPO4
  PO4
  H_2
  OH_2

  Na_1
  CO3
  HCO3
  H2CO3
  H_3
  OH_3

  Na_2
  SO3
  HSO3
  H2SO3
  H_4
  OH_4

  CL
  NH3_H2O
  NH4
  H_5
  OH_5

  Na_3
  HAC
  AC
  H_6
  OH_6

  Na_4
  HCN
  CN
  H_7
  OH_7
]

to setup
  clear-all
  ;设置画布颜色
   ask patches [
    set pcolor white
  ]

  set-default-shape mode6_h2c2o4_t  "h2c2o4"
  set-default-shape mode6_hc2o4_t  "hc2o4"
  set-default-shape mode6_c2o4_t  "c2o4"
  set-default-shape mode6_h_t  "h"
  set-default-shape mode6_oh_t  "oh"
  set-default-shape mode6_h3po4_t  "h3po4"
  set-default-shape mode6_h2po4_t  "h2po4"
  set-default-shape mode6_hpo4_t  "hpo4"
  set-default-shape mode6_po4_t  "po4"
  set-default-shape mode6_na_t  "na"
  set-default-shape mode6_co3_t  "co3"
  set-default-shape mode6_hco3_t  "hco3"
  set-default-shape mode6_h2co3_t  "h2co3"
  set-default-shape mode6_so3_t  "so3"
  set-default-shape mode6_hso3_t  "hso3"
  set-default-shape mode6_h2so3_t  "h2so3"
  set-default-shape mode6_nh3h2o_t  "nh3h2o"
  set-default-shape mode6_nh4_t  "nh4"
  set-default-shape mode6_cl_t  "cl"
  set-default-shape mode6_hac_t  "hac"
  set-default-shape mode6_ac_t  "ac"
  set-default-shape mode6_hcn_t  "hcn"
  set-default-shape mode6_cn_t  "cn"

  ; 加载 Python 脚本
  py:setup
  py:python

  if simple = "H3PO4" [
    let resul3 H3PO4_caculate
    set H3PO4 item 0 H3PO4_caculate
    set H2PO4 item 1 H3PO4_caculate
    set HPO4 item 2 H3PO4_caculate
    set PO4 item 4 H3PO4_caculate
    set H_2 item 3 H3PO4_caculate
    set OH_2 item 5 H3PO4_caculate
    create-mode6_h3po4_t C0 * 1000 [
      setxy random-xcor random-ycor
      set size 2.5
    ]
  ]
  if simple = "H2C2O4" [
    let resul3 H2C2O4_caculate
    set H2C2O4 item 0 H2C2O4_caculate
    set HC2O4 item 1 H2C2O4_caculate
    set C2O4 item 2 H2C2O4_caculate
    set H_1 item 3 H2C2O4_caculate
    set OH_1 item 4 H2C2O4_caculate
    create-mode6_h2c2o4_t C0 * 1000 [
      setxy random-xcor random-ycor
      set size 2.5
    ]
    print (resul3)
  ]
  if simple = "Na2CO3" [
    let resul3 Na2CO3_caculate
    set Na_1 2 * C0
    set CO3 item 0 Na2CO3_caculate
    set HCO3 item 1 Na2CO3_caculate
    set H2CO3 item 2 Na2CO3_caculate
    set H_3 item 3 Na2CO3_caculate
    set OH_3 item 4 Na2CO3_caculate
    create-mode6_na_t C0 * 1000 [
      setxy random-xcor random-ycor
      set size 2.5
    ]
    print (resul3)
  ]
  if simple = "Na2SO3" [
    let resul3 Na2SO3_caculate
    set Na_2 2 * C0
    set SO3 item 0 Na2SO3_caculate
    set HSO3 item 1 Na2SO3_caculate
    set H2SO3 item 2 Na2SO3_caculate
    set H_4 item 3 Na2SO3_caculate
    set OH_4 item 4 Na2SO3_caculate
    create-mode6_na_t C0 * 1000 [
      setxy random-xcor random-ycor
      set size 2.5
    ]
    print (resul3)
  ]
  if simple = "NaHCO3" [
    let resul3 NaHCO3_caculate
    set Na_1 C0
    set CO3 item 2 NaHCO3_caculate
    set HCO3 item 0 NaHCO3_caculate
    set H2CO3 item 1 NaHCO3_caculate
    set H_3 item 3 NaHCO3_caculate
    set OH_3 item 4 NaHCO3_caculate
    create-mode6_na_t C0 * 1000 [
      setxy random-xcor random-ycor
      set size 2.5
    ]
    print (resul3)
  ]
  if simple = "NaHSO3" [
    let resul3 NaHSO3_caculate
    set Na_2 C0
    set SO3 item 2 NaHSO3_caculate
    set HSO3 item 0 NaHSO3_caculate
    set H2SO3 item 1 NaHSO3_caculate
    set H_4 item 3 NaHSO3_caculate
    set OH_4 item 4 NaHSO3_caculate
    create-mode6_na_t C0 * 1000 [
      setxy random-xcor random-ycor
      set size 2.5
    ]
    print (resul3)
  ]
  if mixed = "NH3.H2O + NH4Cl" [
    let resul3 NH3H2O_NH4Cl_caculate
    set CL C0
    set NH3_H2O item 0 NH3H2O_NH4Cl_caculate
    set  NH4 item 1 NH3H2O_NH4Cl_caculate
    set  H_5 item 2 NH3H2O_NH4Cl_caculate
    set  OH_5 item 3 NH3H2O_NH4Cl_caculate
    create-mode6_nh3h2o_t C0 * 1000 [
      setxy random-xcor random-ycor
      set size 2.5
    ]
    print (resul3)
  ]
  if mixed = "CH3COOH + CH3COONa" [
    let resul3 CH3COOH_CH3COONa_caculate
    set Na_3 C0
    set HAC item 0 CH3COOH_CH3COONa_caculate
    set AC item 1 CH3COOH_CH3COONa_caculate
    set H_6 item 2 CH3COOH_CH3COONa_caculate
    set OH_6 item 3 CH3COOH_CH3COONa_caculate
    create-mode6_na_t C0 * 1000 [
      setxy random-xcor random-ycor
      set size 2.5
    ]
    print (resul3)
  ]
  if mixed = "HCN + NaCN" [
    let resul3 HCN_NaCN_caculate
    set Na_4 C0
    set HCN item 0 HCN_NaCN_caculate
    set CN item 1 HCN_NaCN_caculate
    set H_7 item 2 HCN_NaCN_caculate
    set OH_7 item 3 HCN_NaCN_caculate
    create-mode6_na_t C0 * 1000 [
      setxy random-xcor random-ycor
      set size 2.5
    ]
    print (resul3)
  ]

  if simple = "H2C2O4" [
    let temp_h2c2o4  (H2C2O4 / (H2C2O4 + HC2O4 + C2O4 + H_1 + OH_1))
    set temp_h2c2o4 (temp_h2c2o4 * C0 * 1000)
    let temp_hc2o4  (HC2O4 / (H2C2O4 + HC2O4 + C2O4 + H_1 + OH_1))
    set temp_hc2o4 (temp_hc2o4 * C0 * 1000)
    let temp_c2o4  (C2O4 / (H2C2O4 + HC2O4 + C2O4 + H_1 + OH_1))
    set temp_c2o4 (temp_c2o4 * C0 * 1000)
    let temp_h  (H_1 / (H2C2O4 + HC2O4 + C2O4 + H_1 + OH_1))
    set temp_h (temp_h * C0 * 1000)
    let temp_oh  (OH_1 / (H2C2O4 + HC2O4 + C2O4 + H_1 + OH_1))
    set temp_oh (temp_oh * C0 * 1000)

    let temp_i 1
    let cent (list temp_h2c2o4 temp_hc2o4 temp_c2o4 temp_h temp_oh)
    let cent-list sort cent
    foreach cent-list [
      value  ->
      if value = temp_h2c2o4[
      ]
      if value = temp_hc2o4[
        ask n-of ((round value) + temp_i ) mode6_h2c2o4_t[
          set breed mode6_hc2o4_t
        ]
      ]
      if value = temp_c2o4[
        ask n-of ((round value) + temp_i ) mode6_h2c2o4_t[
          set breed mode6_c2o4_t
        ]
      ]
      if value = temp_h[
        ask n-of ((round value) + temp_i ) mode6_h2c2o4_t[
          set breed mode6_h_t
        ]
      ]
      if value = temp_oh[
        ask n-of ((round value) + temp_i ) mode6_h2c2o4_t[
          set breed mode6_oh_t
        ]
      ]
      set temp_i (temp_i + 1)
    ]
  ]
  if simple = "H3PO4" [
    let temp_h3po4  (H3PO4 / (H3PO4 + H2PO4 + HPO4 + PO4 + H_2 + OH_2))
    set temp_h3po4 (temp_h3po4 * C0 * 1000)
    let temp_h2po4  (H2PO4 / (H3PO4 + H2PO4 + HPO4 + PO4 + H_2 + OH_2))
    set temp_h2po4 (temp_h2po4 * C0 * 1000)
    let temp_hpo4  (HPO4 / (H3PO4 + H2PO4 + HPO4 + PO4 + H_2 + OH_2))
    set temp_hpo4 (temp_hpo4 * C0 * 1000)
    let temp_po4  (PO4 / (H3PO4 + H2PO4 + HPO4 + PO4 + H_2 + OH_2))
    set temp_po4 (temp_po4 * C0 * 1000)
    let temp_h  (H_2 / (H3PO4 + H2PO4 + HPO4 + PO4 + H_2 + OH_2))
    set temp_h (temp_h * C0 * 1000)
    let temp_oh  (OH_2 / (H3PO4 + H2PO4 + HPO4 + PO4 + H_2 + OH_2))
    set temp_oh (temp_oh * C0 * 1000)

    let temp_i 1
    let cent (list temp_h3po4 temp_h2po4 temp_hpo4 temp_po4 temp_h temp_oh)
    let cent-list sort cent
    foreach cent-list [
      value  ->
      if value = temp_h3po4[
      ]
      if value = temp_h2po4[
        ask n-of ((round value) + temp_i ) mode6_h3po4_t[
          set breed mode6_h2po4_t
        ]
      ]
      if value = temp_hpo4[
        ask n-of ((round value) + temp_i ) mode6_h3po4_t[
          set breed mode6_hpo4_t
        ]
      ]
      if value = temp_po4[
        ask n-of ((round value) + temp_i ) mode6_h3po4_t[
          set breed mode6_po4_t
        ]
      ]
      if value = temp_h[
        ask n-of ((round value) + temp_i ) mode6_h3po4_t[
          set breed mode6_h_t
        ]
      ]
      if value = temp_oh[
        ask n-of ((round value) + temp_i ) mode6_h3po4_t[
          set breed mode6_oh_t
        ]
      ]
      set temp_i (temp_i + 1)
    ]
  ]
  if simple = "Na2CO3" or simple = "NaHCO3" [
    let temp_na  (Na_1 / (Na_1 + CO3 + HCO3 + H2CO3 + H_3 + OH_3))
    set temp_na (temp_na * C0 * 1000)
    let temp_co3 (CO3 / (Na_1 + CO3 + HCO3 + H2CO3 + H_3 + OH_3))
    set temp_co3 (temp_co3 * C0 * 1000)
    let temp_hco3  (HCO3 / (Na_1 + CO3 + HCO3 + H2CO3 + H_3 + OH_3))
    set temp_hco3 (temp_hco3 * C0 * 1000)
    let temp_h2co3  (H2CO3 / (Na_1 + CO3 + HCO3 + H2CO3 + H_3 + OH_3))
    set temp_h2co3 (temp_h2co3 * C0 * 1000)
    let temp_h  (H_3 / (Na_1 + CO3 + HCO3 + H2CO3 + H_3 + OH_3))
    set temp_h (temp_h * C0 * 1000)
    let temp_oh  (OH_3 / (Na_1 + CO3 + HCO3 + H2CO3 + H_3 + OH_3))
    set temp_oh (temp_oh * C0 * 1000)

    let temp_i 1
    let cent (list temp_na temp_co3 temp_hco3 temp_h2co3 temp_h temp_oh)
    let cent-list sort cent
    foreach cent-list [
      value  ->
      if value = temp_na[
      ]
      if value = temp_co3[
        ask n-of ((round value) + temp_i ) mode6_na_t[
          set breed mode6_co3_t
        ]
      ]
      if value = temp_hco3[
        ask n-of ((round value) + temp_i ) mode6_na_t[
          set breed mode6_hco3_t
        ]
      ]
      if value = temp_h2co3[
        ask n-of ((round value) + temp_i ) mode6_na_t[
          set breed mode6_h2co3_t
        ]
      ]
      if value = temp_h[
        ask n-of ((round value) + temp_i ) mode6_na_t[
          set breed mode6_h_t
        ]
      ]
      if value = temp_oh[
        ask n-of ((round value) + temp_i ) mode6_na_t[
          set breed mode6_oh_t
        ]
      ]
      set temp_i (temp_i + 1)
    ]
  ]
  if simple = "Na2SO3" or simple = "NaHSO3" [
    let temp_na  (Na_2 / (Na_2 + SO3 + HSO3 + H2SO3 + H_4 + OH_4))
    set temp_na (temp_na * C0 * 1000)
    let temp_so3 (SO3 / (Na_2 + SO3 + HSO3 + H2SO3 + H_4 + OH_4))
    set temp_so3 (temp_so3 * C0 * 1000)
    let temp_hso3 (HSO3 / (Na_2 + SO3 + HSO3 + H2SO3 + H_4 + OH_4))
    set temp_hso3 (temp_hso3 * C0 * 1000)
    let temp_h2so3 (H2SO3 / (Na_2 + SO3 + HSO3 + H2SO3 + H_4 + OH_4))
    set temp_h2so3 (temp_h2so3 * C0 * 1000)
    let temp_h  (H_4 / (Na_2 + SO3 + HSO3 + H2SO3 + H_4 + OH_4))
    set temp_h (temp_h * C0 * 1000)
    let temp_oh (OH_4 / (Na_2 + SO3 + HSO3 + H2SO3 + H_4 + OH_4))
    set temp_oh (temp_oh * C0 * 1000)

    let temp_i 1
    let cent (list temp_na temp_so3 temp_hso3 temp_h2so3 temp_h temp_oh)
    let cent-list sort cent
    foreach cent-list [
      value  ->
      if value = temp_na[
      ]
      if value = temp_so3[
        ask n-of ((round value) + temp_i ) mode6_na_t[
          set breed mode6_so3_t
        ]
      ]
      if value = temp_hso3[
        ask n-of ((round value) + temp_i ) mode6_na_t[
          set breed mode6_hso3_t
        ]
      ]
      if value = temp_h2so3[
        ask n-of ((round value) + temp_i ) mode6_na_t[
          set breed mode6_h2so3_t
        ]
      ]
      if value = temp_h[
        ask n-of ((round value) + temp_i ) mode6_na_t[
          set breed mode6_h_t
        ]
      ]
      if value = temp_oh[
        ask n-of ((round value) + temp_i ) mode6_na_t[
          set breed mode6_oh_t
        ]
      ]
      set temp_i (temp_i + 1)
    ]
  ]
  if mixed = "NH3.H2O + NH4Cl" [
    let temp_cl  (CL / (CL + NH3_H2O + NH4 + H_5 + OH_5))
    set temp_cl (temp_cl * C0 * 1000)
    let temp_nh3h2o (NH3_H2O / (CL + NH3_H2O + NH4 + H_5 + OH_5))
    set temp_nh3h2o (temp_nh3h2o * C0 * 1000)
    let temp_nh4 (NH4 / (CL + NH3_H2O + NH4 + H_5 + OH_5))
    set temp_nh4 (temp_nh4 * C0 * 1000)
    let temp_h  (H_5 / (CL + NH3_H2O + NH4 + H_5 + OH_5))
    set temp_h (temp_h * C0 * 1000)
    let temp_oh (OH_5 / (CL + NH3_H2O + NH4 + H_5 + OH_5))
    set temp_oh (temp_oh * C0 * 1000)

    let temp_i 1
    let cent (list temp_cl temp_nh3h2o temp_nh4 temp_h temp_oh)
    let cent-list sort cent
    foreach cent-list [
      value  ->
      if value = temp_nh3h2o[
      ]
      if value = temp_cl[
        ask n-of ((round value) + temp_i ) mode6_nh3h2o_t[
          set breed mode6_cl_t
        ]
      ]
      if value = temp_nh4[
        ask n-of ((round value) + temp_i ) mode6_nh3h2o_t[
          set breed mode6_nh4_t
        ]
      ]
      if value = temp_h[
        ask n-of ((round value) + temp_i ) mode6_nh3h2o_t[
          set breed mode6_h_t
        ]
      ]
      if value = temp_oh[
        ask n-of ((round value) + temp_i ) mode6_nh3h2o_t[
          set breed mode6_oh_t
        ]
      ]
      set temp_i (temp_i + 1)
    ]
  ]
  if mixed = "CH3COOH + CH3COONa" [
    let temp_na  (Na_3 / (Na_3 + HAC + AC + H_6 + OH_6))
    set temp_na (temp_na * C0 * 1000)
    let temp_hac (HAC / (Na_3 + HAC + AC + H_6 + OH_6))
    set temp_hac (temp_hac * C0 * 1000)
    let temp_ac (AC / (Na_3 + HAC + AC + H_6 + OH_6))
    set temp_ac (temp_ac * C0 * 1000)
    let temp_h  (H_6 / (Na_3 + HAC + AC + H_6 + OH_6))
    set temp_h (temp_h * C0 * 1000)
    let temp_oh (OH_6 / (Na_3 + HAC + AC + H_6 + OH_6))
    set temp_oh (temp_oh * C0 * 1000)

    let temp_i 1
    let cent (list temp_na temp_hac temp_ac temp_h temp_oh)
    let cent-list sort cent
    foreach cent-list [
      value  ->
      if value = temp_na[
      ]
      if value = temp_hac[
        ask n-of ((round value) + temp_i ) mode6_na_t[
          set breed mode6_hac_t
        ]
      ]
      if value = temp_ac[
        ask n-of ((round value) + temp_i ) mode6_na_t[
          set breed mode6_ac_t
        ]
      ]
      if value = temp_h[
        ask n-of ((round value) + temp_i ) mode6_na_t[
          set breed mode6_h_t
        ]
      ]
      if value = temp_oh[
        ask n-of ((round value) + temp_i ) mode6_na_t[
          set breed mode6_oh_t
        ]
      ]
      set temp_i (temp_i + 1)
    ]
  ]
  if mixed = "HCN + NaCN" [
    let temp_na  (Na_4 / (Na_4 + HCN + CN + H_7 + OH_7))
    set temp_na (temp_na * C0 * 1000)
    let temp_hcn (HCN / (Na_4 + HCN + CN + H_7 + OH_7))
    set temp_hcn (temp_hcn * C0 * 1000)
    let temp_cn (CN / (Na_4 + HCN + CN + H_7 + OH_7))
    set temp_cn (temp_cn * C0 * 1000)
    let temp_h  (H_7 / (Na_4 + HCN + CN + H_7 + OH_7))
    set temp_h (temp_h * C0 * 1000)
    let temp_oh (OH_7 / (Na_4 + HCN + CN + H_7 + OH_7))
    set temp_oh (temp_oh * C0 * 1000)

    let temp_i 1
    let cent (list temp_na temp_hcn temp_cn temp_h temp_oh)
    let cent-list sort cent
    foreach cent-list [
      value  ->
      if value = temp_na[
      ]
      if value = temp_hcn[
        ask n-of ((round value) + temp_i ) mode6_na_t[
          set breed mode6_hcn_t
        ]
      ]
      if value = temp_cn[
        ask n-of ((round value) + temp_i ) mode6_na_t[
          set breed mode6_cn_t
        ]
      ]
      if value = temp_h[
        ask n-of ((round value) + temp_i ) mode6_na_t[
          set breed mode6_h_t
        ]
      ]
      if value = temp_oh[
        ask n-of ((round value) + temp_i ) mode6_na_t[
          set breed mode6_oh_t
        ]
      ]
      set temp_i (temp_i + 1)
    ]
  ]
end 

to go
  ask turtles
      [ fd 1                                        ;; move turtles around randomly
        rt random 360
        lt random 360 ]    ;; around the world
end 

to-report H3PO4_caculate
    set a_ C0
    set b_ a_
  while[true][
    set g_ 7.1 * 10 ^ (-3)
    set h_ 6.3 * 10 ^ (-8)
    set l_ 4.8 * 10 ^ (-13)
    set k_ 10 ^ (-14)
    set c_  calculate 1 0 (- ((g_ * b_ )+ k_)) (- 2 * g_ * h_ * b_) (- 6 * g_ * h_ * l_ * b_)
    let c_t first_positive_value c_
    set d_ (k_ / c_t)
    set e_ ((g_ * b_) / c_t)
    set f_ ((h_ * e_) / c_t)
    set m_ ((l_ * f_) / c_t)
    if a_ = (b_ + e_ + f_ + m_) [
      let result (list b_ e_ f_ c_t m_ d_)
      report result
    ]
    set b_ (a_ - e_ - f_ - m_)
  ]
end 

to-report H2C2O4_caculate
  set a_ C0
  set b_ a_
  while[true][
    set g_ 5.36 * 10 ^ (-2)
    set h_ 5.35 * 10 ^ (-5)
    set k_ 10 ^ (-14)
    set c_  calculate 0 1 0 (- (k_ + (g_ * b_))) (- 2 * g_ * h_ * b_)
    let c_t first_positive_value c_
    set d_ (k_ / c_t)
    set e_ ((g_ * b_) / c_t)
    set f_ ((h_ * e_) / c_t)
    if a_ = (b_ + e_ + f_ ) [
      let result (list b_ e_ f_ c_t d_)
      report result
    ]
    set b_ (a_ - e_ - f_ )
  ]
end 

to-report Na2CO3_caculate
  set a_ C0
  set b_ a_
  while[true][
    set g_ 4.45 * 10 ^ (-7)
    set h_ 4.69 * 10 ^ (-11)
    set k_ 10 ^ (-14)
    set c_  calculate 0 (b_ / (g_ * h_)) (1 + (b_ / h_)) 0 (- k_ )
    let c_t first_positive_value c_
    set d_ (k_ / c_t)
    set e_ ((c_t * b_) / h_)
    set f_ ((c_t * e_) / g_)
    if a_ = (b_ + e_ + f_ ) [
      let result (list b_ e_ f_ c_t d_)
      report result
    ]
    set b_ (a_ - e_ - f_ )
  ]
end 

to-report Na2SO3_caculate
  set a_ C0
  set b_ a_
  while[true][
    set g_ 1.29 * 10 ^ (-2)
    set h_ 6.24 * 10 ^ (-8)
    set k_ 10 ^ (-14)
    set c_  calculate 0 (b_ / (g_ * h_)) (1 + (b_ / h_)) 0 (- k_ )
    let c_t first_positive_value c_
    set d_ (k_ / c_t)
    set e_ ((c_t * b_) / h_)
    set f_ ((c_t * e_) / g_)
    if a_ = (b_ + e_ + f_ ) [
      let result (list b_ e_ f_ c_t d_)
      report result
    ]
    set b_ (a_ - e_ - f_ )
  ]
end 

to-report NaHCO3_caculate
  set a_ C0
  set b_ a_
  while[true][
    set g_ 4.45 * 10 ^ (-7)
    set h_ 4.69 * 10 ^ (-11)
    set k_ 10 ^ (-14)
    set c_ calculate 0 0 (1 + b_ / g_) 0 ( - (h_ * b_ + k_))
    let c_t first_positive_value c_
    set d_ (k_ / c_t)
    set e_ ((c_t * b_) / g_)
    set f_ ((h_ * b_) / c_t)
    if a_ = (b_ + e_ + f_) [
      let result (list b_ e_ f_ c_t d_)
      report result
    ]
      set b_ (a_ - e_ - f_)
  ]
end 

to-report NaHSO3_caculate
  set a_ C0
  set b_ a_
  while[true][
    set g_ 1.29 * 10 ^ (-2)
    set h_ 6.24 * 10 ^ (-8)
    set k_ 10 ^ (-14)
    set c_ calculate 0 0 (1 + b_ / g_) 0 ( - (h_ * b_ + k_))
    let c_t first_positive_value c_
    set d_ (k_ / c_t)
    set e_ ((c_t * b_) / g_)
    set f_ ((h_ * b_) / c_t)
    if a_ = (b_ + e_ + f_) [
      let result (list b_ e_ f_ c_t d_)
      report result
    ]
      set b_ (a_ - e_ - f_)
  ]
end 

to-report NH3H2O_NH4Cl_caculate
  set a_ C0
  set b_ a_
  while[true][
    set g_ 1.76 * 10 ^ (-5)
    set k_ 10 ^ (-14)
    set d_ calculate 0 0 2 b_ (- (2 * k_ + g_ * b_))
    let d_t first_positive_value d_
    print d_t
    set c_ (k_ / d_t)
    set e_ ((g_ / (g_ + d_t)) * 2 * a_ )
    if  a_ = ((b_ + e_) / 2) [
      let result (list b_ e_ c_ d_t)
      report result
    ]
      set b_ ( 2 * a_ - e_)
  ]
end 

to-report CH3COOH_CH3COONa_caculate
  set a_ C0
  set b_ a_
  while[true][
    set g_ 1.75 * 10 ^ (-5)
    set k_ 10 ^ (-14)
    set c_ calculate 0 0 2 b_ (- (2 * k_ + g_ * b_))
    let c_t first_positive_value c_
    print c_t
    set d_ (k_ / c_t)
    set e_ ((g_ / (g_ + c_t)) * 2 * a_ )
    if  a_ = ((b_ + e_) / 2) [
      let result (list b_ e_ c_t d_)
      report result
    ]
      set b_ ( 2 * a_ - e_)
  ]
end 

to-report HCN_NaCN_caculate
  set a_ C0
  set b_ a_
  while[true][
    set g_ 6.17 * 10 ^ (-10)
    set k_ 10 ^ (-14)
    set c_ calculate 0 0 2 b_ (- (2 * k_ + g_ * b_))
    let c_t first_positive_value c_
    print c_t
    set d_ (k_ / c_t)
    set e_ ((g_ / (g_ + c_t)) * 2 * a_ )
    if  a_ = ((b_ + e_) / 2) [
      let result (list b_ e_ c_t d_)
      report result
    ]
      set b_ ( 2 * a_ - e_)
  ]
end 

;;计算四次函数，
;;a,b,c,d,f分别代表四次，三次，二次，一次，常数前面的系数

to-report calculate [x1 x2 x3 x4 x5]
  print x1
  print x2
  print x3
  print x4
  print x5
  py:run "import numpy as np"
  py:run "from sympy import symbols, solveset, re"
  py:run "x = symbols('x')"
  py:set "param1" x1
  py:set "param2" x2
  py:set "param3" x3
  py:set "param4" x4
  py:set "param5" x5
  py:run "equation = param1*x**4 + param2*x**3 + param3*x**2 + param4*x + param5"
  py:run "solutions = solveset(equation, x)"
  py:run "real_solutions = [float(sol.evalf(n=30)) for sol in solutions if sol.is_real]"
  report py:runresult "real_solutions"
end 

;;返回list中第一个大于0的值

to-report first_positive_value [my-list]
  print my-list
  foreach my-list [
    xy ->
    if xy >= 0[
      report xy
    ]
  ]
end