eqcurve_extractor.cpp
No OneTemporary
Actions

Subscribers

None

File Metadata

	#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