eqcurve_extractor.cpp
9fc5d580b72bfixpoint

eqcurve_extractor.cpp
View Options

			#include "eqcurve_extractor.hpp"

			#include <iostream>

			namespace specmicp {
			namespace reactmicp {
			namespace eqcurve {


			index_t EquilibriumCurveExtractor::find_point(scalar_t x)
			{
			index_t index;
			const index_t size = m_eqcurve.rows();
			index_t lower_index = 0;
			index_t upper_index = size+1;

			while (upper_index - lower_index > 1)
			{
			index = (upper_index + lower_index) /2;
			if ((x >= xs_offset(index)) == is_increasing())
			lower_index = index;
			else
			upper_index = index;
			}
			if (x == xs_offset(1)) index = 0;
			else if (x == xs_offset(size)) index = size-1;
			else index = std::max((specmicp::index_t) 0, lower_index-1);

			return index;
			}


			scalar_t EquilibriumCurveExtractor::interpolate(index_t index, scalar_t x, index_t col)
			{
			// limit
			if (is_increasing())
			{
			if (x <= m_eqcurve(first(), 0))
			return m_eqcurve(first(), col);
			else if (x >= m_eqcurve(last(), 0))
			return m_eqcurve(last(), col);
			}
			else
			{
			if (x >= m_eqcurve(first(), 0))
			return m_eqcurve(first(), col);
			else if (x <= m_eqcurve(last(), 0))
			return m_eqcurve(last(), col);
			}
			scalar_t y;

			const scalar_t diff = slope(index, col);
			y = m_eqcurve(index, col) + diff*(m_eqcurve(index+1,0)-x);

			return y;
			}

			scalar_t EquilibriumCurveExtractor::slope(index_t index, index_t col)
			{
			if (index == last()) return (m_eqcurve(index-1,col)-m_eqcurve(index,col))/(m_eqcurve(index-1,0)-m_eqcurve(index,0));
			return (m_eqcurve(index+1,col)-m_eqcurve(index,col))/(m_eqcurve(index+1,0)-m_eqcurve(index,0));
			}

			} // end namespace eqcurve
			} // end namespace reactmicp
			} // end namespace specmicp