化学离子平衡作业偷懒神器

下一次加入沉淀平衡

这一段放到K.dat中，有需要自己添加

1e-14
4
1 CH3COOH 1.74e-5 CH3COO-
2 H2CO3 4.2e-7 HCO3- 5.6e-11 CO3--
2 H2S 1.3e-7 HS- 7.1e-15 S--
1 HCl 1e100 Cl-
3
1 NH3_H2O 1.78e-5 NH4+
3 Fe(OH)3 3.16e-9 Fe(OH)2- 5.01e-10 Fe(OH)-- 1.35e-12 Fe---
1 NaOH 1e100 Na+

以下为main.cpp

#include<bits/stdc++.h>
using namespace std;
const double eps = 1e-14;
double Kacid[110][5], Kbase[110][5];
int Na, Nb;
int Tacid[110], Tbase[110];

double Cacid[110][5], Cbase[110][5];
double CH, COH;
double Kw; // the ionic product of water
string nameAcid[110][5];
string nameBase[110][5];

void readData(){
	ifstream in("K.dat");
	in >> Kw;
	in >> Na;
	for (int i=1;i<=Na;i++){
		in >> Tacid[i];
		for (int j=1;j<=Tacid[i];j++){
			in >> nameAcid[i][j];
			in >> Kacid[i][j];
		}
		in >> nameAcid[i][Tacid[i]+1];
	}
	in >> Nb;
	for (int i=1;i<=Nb;i++){
		in >> Tbase[i];
		for (int j=1;j<=Tbase[i];j++){
			in >> nameBase[i][j];
			in >> Kbase[i][j];
		}
		in >> nameBase[i][Tbase[i]+1];
	}
	in.close();
}

double Solve(double &C1, double &C2, double &C3, double &C4, double K){
	if (K > 1e30){
		double T = min(C1, C2);
		C3 += T;
		C4 += T;
		C1 -= T;
		C2 -= T;
		return 0;
	}
	else if (K < 1e-30){
		double T = min(C3, C4);
		C3 -= T;
		C4 -= T;
		C1 += T;
		C2 += T;
		return 0;
	}
	double A = K-1;
	double B = -K*(C1+C2)-(C3+C4);
	double C = K*C1*C2 - C3*C4;
//	cout << A << ' ' << B << ' ' << C << endl;
	assert(B*B - 4 * A * C + eps >= 0);
	double X1 = (-B + sqrt(B*B-4*A*C)) / 2 / A;
	double X2 = (-B - sqrt(B*B-4*A*C)) / 2 / A;
	if (abs(X1) < eps && abs(X2) < eps) return 0;
//	cout << A << ' ' << B << ' ' << C << ' ' << X1 << ' ' << X2 << ' ' << K << endl;
	if (X1 < min(C1,C2) + eps && -X1 < min(C3,C4) + eps) {
		C1 -= X1; C2 -= X1; C3 += X1; C4 += X1;
		return X1;
	}
	if (X2 < min(C1,C2) + eps && -X2 < min(C3,C4) + eps) {
		C1 -= X2; C2 -= X2; C3 += X2; C4 += X2;
		return X2;
	}
	assert(0);
	return 0;
}

double Solve(double &C1, double &C2, double &C3, double K){
	if (K > 1e30){
		double T = C1;
		C3 += T;
		C1 -= T;
		C2 += T;
		return 0;
	}
	else if (K < 1e-30){
		double T = min(C3, C2);
		C1 += T;
		C2 -= T;
		C3 -= T;
		return 0;
	}
	double A = 1;
	double B = C2+C3+K*C1;
	double C = C2*C3-K*C1;
	assert(B*B - 4 * A * C + eps >= 0);
	double X1 = (-B + sqrt(B*B-4*A*C)) / 2 / A;
	double X2 = (-B - sqrt(B*B-4*A*C)) / 2 / A;
//	cout << X1 << ' ' << X2 << ' ' << K << endl;
	if (abs(X1) < eps && abs(X2) < eps) return 0;
	if (X1 < C1 + eps && -X1 < min(C2,C3) + eps) {
		C1 -= X1; C2 += X1; C3 += X1;
		return X1;
	}
	if (X2 < C1 + eps && -X2 < min(C2,C3) + eps) {
		C1 -= X2; C2 += X2; C3 += X2;
		return X2;
	}
	assert(0);
	return 0;
}

double SolveWater(){
	double A = 1;
	double B = -CH - COH;
	double C = CH * COH - Kw;
 // A = 1 B = -([H] + [OH]) C = [H][OH]-Kw;
	assert(B*B - 4 * A * C >= 0);
	double X1 = (-B + sqrt(B*B-4*A*C)) / 2 / A;
	double X2 = (-B - sqrt(B*B-4*A*C)) / 2 / A;
	if (X1 < min(CH,COH) + eps) {
		CH -= X1; COH -= X1;
		return X1;
	}
	if (X2 < min(CH,COH) + eps) {
		CH -= X2; COH -= X2;
		return X2;
	}
	assert(0);
	return 0;
}

void calc(){
	/* Find the reaction that can be do max*/
//	int Ka = 1, Kal = 1, Kb = 1, Kbl = 1;
//	double val = 0;
	for (int i=1;i<=Na;i++){
		for (int j=1;j<=Tacid[i];j++){
			for (int k=1;k<=Nb;k++){
				for (int l=1;l<=Tbase[k];l++){
						double K = Kacid[i][j] * Kbase[k][l] / Kw;
//						cout << Cacid[i][j] << ' ' << Cbase[k][l] << endl;
						Solve(Cacid[i][j], Cbase[k][l], Cacid[i][j+1], Cbase[k][l+1], K);
				}
			}
		}
	}
	
	for (int i=1;i<=Na;i++){
		for (int j=1;j<=Tacid[i];j++){
			double K = Kacid[i][j];
			Solve(Cacid[i][j], Cacid[i][j+1], CH, K);
		}
	}

	for (int i=1;i<=Nb;i++){
		for (int j=1;j<=Tbase[i];j++){
			double K = Kbase[i][j];
			Solve(Cbase[i][j], Cbase[i][j+1], COH, K);
		}
	}
	
	SolveWater();

}
/* 
 HA == A- + H+ Ka
 BOH == B+ + OH- Kb
 H+ + OH- == H2O 1/Kw
 HA + BOH == AB + H2O 
 [A-][B+]   [A-][H+][B+][OH-]
 -------- = ----------------- = Ka * Kb / Kw
  [AB]         [AB][H2O]   
  Sol.
  HA + BOH == A- + B+ + H2O K
  C1   C2     C3   C4
  (C3+x * C4+x) / (C1-x * C2-x) == K
  K(C1C2 - (C1+C2)x + x^2) == C3C4 + (C3+C4)x + x^2
  --> A = (K-1) B = -K*(C1+C2)-(C3+C4) C = KC1C2 - C3C4

  --------------------------------------------------------

  HA == A- + H+ Ka
  C1    C2   C3
   ([A-]+x)([H+]+x)
  ------------------ = K
         [HA]-x
  
  --> A = 1 B = C2+C3+K*C1 C = C2*C3-K*C1

  --------------------------------------------------------

  [H] * [OH] = Kw;
  ([H]-x) * ([OH]-x) = Kw;
  A = 1 B = -([H] + [OH]) C = [H][OH]-Kw;
*/
int main(){
	readData();
	int n;
	int T = 10000;
	Cacid[1][1] = 0.1;
	Cacid[1][2] = 0.1;
	Cbase[3][2] = 0.2;
	while (T--) calc();
	cout << "[H] = " << CH << endl;
	cout << "[OH] = " << COH << endl;
	cout << "pH = " << -log10(CH) << endl;
	cout << "Acid:" << endl;
	for (int i=1;i<=Na;i++){
		for (int j=1;j<=Tacid[i]+1;j++)
			cout << "[" << nameAcid[i][j] << "] = " << Cacid[i][j] << ' ';
		cout << endl;
	}
	cout << "Base" << endl;
	for (int i=1;i<=Nb;i++){
		for (int j=1;j<=Tbase[i]+1;j++)
			cout << "[" << nameBase[i][j] << "] = " << Cbase[i][j] << ' ';
		cout << endl;
	}
}

示例：
0.1mol/L CH3COOH 与 0.1mol/L CH3COONa混，最后pH值

	Cacid[1][1] = 0.1;
	Cacid[1][2] = 0.1;
	Cbase[3][2] = 0.1;

Solution:

[H] = 1.73939e-05
[OH] = 5.74913e-10
pH = 4.7596
Acid:
[CH3COOH] = 0.0999826 [CH3COO-] = 0.100017
[H2CO3] = 0 [HCO3-] = 0 [CO3--] = 0
[H2S] = 0 [HS-] = 0 [S--] = 0
[HCl] = 0 [Cl-] = 0
Base
[NH3_H2O] = 0 [NH4+] = 0
[Fe(OH)3] = 0 [Fe(OH)2-] = 0 [Fe(OH)--] = 0 [Fe---] = 0
[NaOH] = 0 [Na+] = 0.1

示例：
0.01mol/L CH3COOH 与 0.40mol/L CH3COONa混，最后pH值

	Cacid[1][1] = 0.01;
	Cacid[1][2] = 0.4;
	Cbase[3][2] = 0.4;

Solution:

[H] = 4.34982e-07
[OH] = 2.29895e-08
pH = 6.36153
Acid:
[CH3COOH] = 0.00999959 [CH3COO-] = 0.4
[H2CO3] = 0 [HCO3-] = 0 [CO3--] = 0
[H2S] = 0 [HS-] = 0 [S--] = 0
[HCl] = 0 [Cl-] = 0
Base
[NH3_H2O] = 0 [NH4+] = 0
[Fe(OH)3] = 0 [Fe(OH)2-] = 0 [Fe(OH)--] = 0 [Fe---] = 0
[NaOH] = 0 [Na+] = 0.4