vertools/fbound.c at master · ryantkoehler/vertools · GitHub

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
/*
* fbound.c
*
* Copyright 2019 Ryan Koehler, VerdAscend Sciences, ryan@verdascend.com
*
* The programs and source code of the vertools collection are free software.
* They are distributed in the hope that they will be useful,
* WITHOUT ANY WARRANTY OF FITNESS FOR ANY PARTICULAR PURPOSE.
*
* Permission is granted for research, educational, and possibly commercial use
*   and modification as long as 1) Code and any derived works are not
*   redistributed for any fee, and 2) Proper credit is given to the authors.
*   If you wish to include this software in a product, or use it commercially,
*   please contact the authors.
*
* See https://www.verdascend.com/ for more
*
*/

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include "prim.h"
#include "fbound.h"


/*****************************************************************/
double Fraction1D(double conc,double conc2,double temp, double gD)
{
    double aD,bD,cD,qD,eD,xD,climD,cexcD;
    int ns;

    if (conc>conc2)
    {
        climD=conc2;
        cexcD=conc;
    }
    else
    {
        climD=conc;
        cexcD=conc2;
    }
    eD=exp( (gD* -1000.0) / (1.98722*(temp+273.15)) );
    aD=-eD*climD;
    bD=1.0+eD*climD+eD*cexcD;
    cD=-eD*cexcD;

    ns = NUM_SIGN(bD);
    qD=-0.5*(bD+RNUM(ns)*sqrt(pow(bD,2.0)-4.0*aD*cD));

    xD=qD/aD*100.0;
    return xD;
}
/*****************************************************************/
double Fraction2D(double conc,double conc2,double temp ,double gD)
    {
    double aD,bD,cD,qD,eD,xD,climD,cexcD;
    int ns;

    if (conc>conc2)
        {climD=conc2;
        cexcD=conc;}
    else
        {climD=conc;
        cexcD=conc2;}
    eD=exp(-(gD*1000.0)/(1.98722*(temp+273.15)));
    aD=-eD*climD;
    bD=1+eD*climD+eD*cexcD;
    cD=-eD*cexcD;

    ns = NUM_SIGN(bD);
    qD=-0.5*(bD+RNUM(ns)*sqrt(pow(bD,2.0)-4.0*aD*cD));

    xD=cD/qD*100.0;
    return xD;
    }