RateLookup.cpp

Go to the documentation of this file.

/**********************************************************************
** This program is part of 'MOOSE', the
** Messaging Object Oriented Simulation Environment.
**   copyright (C) 2003-2007 Upinder S. Bhalla, Niraj Dudani and NCBS
** It is made available under the terms of the
** GNU Lesser General Public License version 2.1
** See the file COPYING.LIB for the full notice.
**********************************************************************/

#include <vector>
using namespace std;

#include "RateLookup.h"

LookupTable::LookupTable(
    double min, double max, unsigned int nDivs, unsigned int nSpecies )
{
    min_ = min;
    max_ = max;
    // Number of points is 1 more than number of divisions.
    // Then add one more since we may interpolate at the last point in the table.
    nPts_ = nDivs + 1 + 1;
    dx_ = ( max - min ) / nDivs;
    // Every row has 2 entries for each type of gate
    nColumns_ = 2 * nSpecies;

    //~ interpolate_.resize( nSpecies );
    table_.resize( nPts_ * nColumns_ );
}

void LookupTable::addColumns(
    int species,
    const vector< double >& C1,
    const vector< double >& C2 )
    //~ const vector< double >& C2,
    //~ bool interpolate )
{
    vector< double >::const_iterator ic1 = C1.begin();
    vector< double >::const_iterator ic2 = C2.begin();
    vector< double >::iterator iTable = table_.begin() + 2 * species;
    // Loop until last but one point
    for ( unsigned int igrid = 0; igrid < nPts_ - 1 ; ++igrid ) {
        *( iTable )     = *ic1;
        *( iTable + 1 ) = *ic2;

        iTable += nColumns_;
        ++ic1, ++ic2;
    }
    // Then duplicate the last point
    *( iTable )     = C1.back();
    *( iTable + 1 ) = C2.back();

    //~ interpolate_[ species ] = interpolate;
}

void LookupTable::column( unsigned int species, LookupColumn& column )
{
    column.column = 2 * species;
    //~ column.interpolate = interpolate_[ species ];
}

void LookupTable::row( double x, LookupRow& row )
{
    if ( x < min_ )
        x = min_;
    else if ( x > max_ )
        x = max_;

    double div = ( x - min_ ) / dx_;
    unsigned int integer = ( unsigned int )( div );

    row.fraction = div - integer;
    row.row = &( table_.front() ) + integer * nColumns_;
}

void LookupTable::lookup(
    const LookupColumn& column,
    const LookupRow& row,
    double& C1,
    double& C2 )
{
    double a, b;
    double *ap, *bp;

    ap = row.row + column.column;

    //~ if ( ! column.interpolate ) {
        //~ C1 = *ap;
        //~ C2 = *( ap + 1 );
        //~
        //~ return;
    //~ }

    bp = ap + nColumns_;

    a = *ap;
    b = *bp;
    C1 = a + ( b - a ) * row.fraction;

    a = *( ap + 1 );
    b = *( bp + 1 );
    C2 = a + ( b - a ) * row.fraction;
}