scirs2-integrate 0.4.0

Numerical integration module for SciRS2 (scirs2-integrate)
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92

//! Simple Chemical Equilibrium Test
//!
//! A minimal test of the chemical equilibrium functionality to verify it works.

use scirs2_integrate::ode::chemical_equilibrium::{systems, ActivityModel};

#[allow(dead_code)]
fn main() -> Result<(), Box<dyn std::error::Error>> {
    println!("Simple Chemical Equilibrium Test\n");

    // Test 1: Weak acid equilibrium setup
    println!("1. Testing weak acid equilibrium setup");
    let ka: f64 = 1.8e-5; // Acetic acid
    let initial_acid = 0.1;

    let calculator = systems::weak_acid_equilibrium(ka, initial_acid, None);

    println!("   ✅ Calculator created successfully");
    println!("   Species: {:?}", calculator.species_names);
    println!("   Reactions: {:?}", calculator.reaction_names);
    println!("   Ka: {:.2e}", calculator.equilibrium_constants[0]);

    // Simple analytical approximation for weak acid
    let h_analytical = (ka * initial_acid).sqrt();
    let ph_analytical = -h_analytical.log10();
    println!("   Analytical pH approximation: {ph_analytical:.2}");

    // Test 2: Buffer system setup
    println!("\n2. Testing buffer equilibrium setup");
    let buffer_calc = systems::buffer_equilibrium(ka, 0.1, 0.1);

    println!("   ✅ Buffer calculator created");
    println!("   Species: {:?}", buffer_calc.species_names);
    println!("   Ka: {:.2e}", buffer_calc.equilibrium_constants[0]);

    // Henderson-Hasselbalch equation for equal concentrations
    let pka = -ka.log10();
    let ph_hh = pka; // pH = pKa when [A-] = [HA]
    println!("   Henderson-Hasselbalch pH: {ph_hh:.2}");

    // Test 3: Complex formation setup
    println!("\n3. Testing complex formation setup");
    let k_formation = 1e6;
    let complex_calc = systems::complex_formation(k_formation, 0.001, 0.01);

    println!("   ✅ Complex formation calculator created");
    println!("   Species: {:?}", complex_calc.species_names);
    println!(
        "   K_formation: {:.2e}",
        complex_calc.equilibrium_constants[0]
    );

    // Estimate complex formation
    let ligand = 0.01; // L
    let alpha = 1.0 / (1.0 + k_formation * ligand);
    let fraction_free = alpha * 100.0;
    println!("   Estimated fraction of free metal: {fraction_free:.1}%");

    // Test 4: Solubility equilibrium setup
    println!("\n4. Testing solubility equilibrium setup");
    let ksp: f64 = 1.8e-10;
    let solubility_calc = systems::solubility_equilibrium(ksp, 1.0, 1.0);

    println!("   ✅ Solubility calculator created");
    println!("   Species: {:?}", solubility_calc.species_names);
    println!("   Ksp: {:.2e}", solubility_calc.equilibrium_constants[0]);

    // Analytical solubility for 1:1 salt
    let solubility_analytical = ksp.sqrt();
    println!("   Analytical solubility: {solubility_analytical:.3e} M");
    println!(
        "   Solubility in mg/L: {:.2}",
        solubility_analytical * 143.32
    ); // AgCl MW

    // Test 5: Activity coefficient models
    println!("\n5. Testing activity coefficient models");
    let mut test_calc = systems::weak_acid_equilibrium(1e-5, 0.1, None);

    // Test different activity models
    test_calc.set_activity_model(ActivityModel::Ideal);
    println!("   ✅ Ideal activity model set");

    test_calc.set_activity_model(ActivityModel::ExtendedDebyeHuckel);
    println!("   ✅ Extended Debye-Hückel model set");

    println!("   Activity coefficient models working properly");

    println!("\n✅ All chemical equilibrium tests passed!");

    Ok(())
}