ph_sensor.cc

Go to the documentation of this file.

#include "ph_sensor.h"
#include <math.h>


double compute_ph( double wA4, double wB4, double pK1, double pK2 );
double compute_ph_rhs( double pH, double wA4, double wB4, double pK1, double pK2 );


ph_sensor::ph_sensor(dsim_model_setup *setup) : dsim_model(setup)
{
}

ph_sensor::~ph_sensor()
{
}
        
void ph_sensor::initialize_data()
{               
    // Parameter initial values
    
    double pK1  = -6.34969247686806;
    double pK2  = -10.2502636844309;
    double pH   = 0;

    // Two parameters and the pH output
    pK1_dsim = create_parameter("pK1",  sd_type_double, &pK1,   1,1,    "", "Parameter 1");     pK1_dsim.mark_telemetry();
    pK2_dsim = create_parameter("pK2",  sd_type_double, &pK2,   1,1,    "", "Parameter 2");     pK2_dsim.mark_telemetry();
    pH_dsim  = create_output("pH",      sd_type_double, &pH,    1,1,    "", "pH measurement");   pH_dsim.mark_telemetry();

    // Thismeans only call initialize_timestep - no rhs()
    configure_timestep(true,false,false);
}

void ph_sensor::initialization_complete()
{
    wA4_dsim    = parent()->request_local_variable("wA4");
    wB4_dsim    = parent()->request_local_variable("wB4");
}

void ph_sensor::initialize_timestep()
{
    // Get values of dsim variables
    double wA4 = wA4_dsim.value_as_double();
    double wB4 = wB4_dsim.value_as_double();
    double pK1 = pK1_dsim.value_as_double();
    double pK2 = pK2_dsim.value_as_double();
    
    // Compute the ph output value
    double pH = compute_ph( wA4, wB4, pK1, pK2 );
    
    // Set the output variable    
    pH_dsim.set_value_as_double(pH);
}


// Secant method
double compute_ph( double wA4, double wB4, double pK1, double pK2 )
{
    // This insures a change of sign from the initial pH of 0
    const double pH_00 = -1;
    const double tol = 1e-3;
        double fOld = compute_ph_rhs( pH_00, wA4, wB4, pK1, pK2 );
    
    double pH = 0;      
        double delta = pH - pH_00;
        double fNew;
        int n = 0;
        
    // Limited to 100 interations
        while( fabs(delta) > tol )
        {
                fNew  = compute_ph_rhs( pH, wA4, wB4, pK1, pK2 );
                
                delta = -fNew/((fNew-fOld)/delta);
        
        fOld = fNew;
                
                pH   += delta;  
        
        n++;
        
        if( n > 100) break;
        }
    return pH;
}

// Right hand side for the secant method
double compute_ph_rhs( double pH, double wA4, double wB4, double pK1, double pK2 )
{
    double y   = wA4 + pow(10,pH-14) - pow(10,-pH)
               + wB4*(1 + 2*pow(10,pH-pK2))/(1 + pow(10,pK1-pH) + pow(10,pH-pK2));

    return y;
}


extern "C"
{
    dsim_model *PHSensor(void *setup);

        dsim_model *PHSensor(void *setup)
        {
                return new ph_sensor((dsim_model_setup *)setup);
        }
        
}