/*
* Parse the formula
*/
void Parse (char *formula, double *C, double *H, double *O, double *N)
{
int i;
bool k;
double *ptr = NULL;
*C = 0.0;
*H = 0.0;
*O = 0.0;
*N = 0.0;
for (i=0;i<strlen (formula);i++)
{
if (formula[i] == 'C')
{
if (!k && ptr)
(*ptr) ++;
ptr = C;
k = false;
}
else if (formula[i] == 'H')
{
if (!k && ptr)
(*ptr) ++;
ptr = H;
k = false;
}
else if (formula[i] == 'O')
{
if (!k && ptr)
(*ptr) ++;
ptr = O;
k = false;
}
else if (formula[i] == 'N')
{
if (!k && ptr)
(*ptr) ++;
ptr = N;
k = false;
}
else if (formula[i] >= '0' && formula[i] <= '9' && ptr)
{
if (k)
*ptr *= 10;
k = true;
*ptr += formula[i] - '0';
}
}
}
/*
* Calculate a^b
*/
long pow (long a, int b)
{
long c = 1;
while (b--)
c *= a;
return c;
}
/*
* test
*/
double DoIt (struct xpl_s *sample, double x)
{
double CO2, CO, H2, H2O, N2, f;
CO2 = x*(6.0532*log (x) - 1.4403) - 9847.5;
CO = x*(2.8096*log (x) + 2.7605) - 5592.6;
H2 = x*(4.4180*log (x) - 11.9030) - 3953.5;
H2O = x*(9.5308*log (x) - 38.8640) - 4599.3;
N2 = x*(2.9479*log (x) + 1.3671) - 5412.1;
f = sample->products.CO*CO + sample->products.CO2*CO2 + sample->products.H2*H2 + sample->products.H2O*H2O + sample->products.N2*N2 - sample->body.E*sample->body.M;
return f;
}
/*
* Dichotomial algorithm
*/
double Dicho (struct xpl_s *sample, double x1, double x2, double epsilon)
{
double f1, f2, f3, x3;
f1 = DoIt (sample, x1);
f2 = DoIt (sample, x2);
x3 = (x1 + x2)*0.5f;
f3 = DoIt (sample, x3);
if (fabs (f3) < epsilon)
return x3;
if (f1*f3 < 0)
return Dicho (sample, x1, x3, epsilon);
else if (f2*f3 < 0)
return Dicho (sample, x3, x2, epsilon);
return 0.0;
}
/*
* Equilibrium
*/
void Equilibrium (double *pa, double *pb, double *pc, double *pd, double k)
{
double q;
double a, b, c, d, x1, x2, x;
q = *pa**pb;
if (*pc && *pd)
q /= (*pc**pd);
else
q = 9999;
a = 1-k;
b = *pc + *pd - *pa - *pb;
c = *pa**pb - k**pc**pd;
d = b*b - 4*a*c;
if (d >= 0)
{
x1 = (-b + sqrt (d))/(2*a);
x2 = (-b - sqrt (d))/(2*a);
}
else
return;
if (q > k)
x = x1 < x2 ? x2 : x1;
else
x = x1 > x2 ? x2 : x1;
*pa -= x;
*pb -= x;
*pc += x;
*pd += x;
}
/*
* SetupEquilibria
*/
void SetupEquilibria (struct xpl_s *sample)
{
double k, t;
do
{
t = sample->body.T;
k = 5.5935*log (t) - 37.348;
Equilibrium (&sample->products.CO, &sample->products.H2O, &sample->products.CO2, &sample->products.H2, k);
sample->body.E = CalculateEnergy (sample);
sample->body.T = CalculateTemperature (sample);
} while (fabs (t - sample->body.T) > 10.0);
} |
.
)
Just kidding...
(I'm not trying to say your program sucks though
What a pain in the ass.