Page MenuHomec4science
Files F90797095

eqcurve_extractor.cpp
No OneTemporary
Actions

File Metadata

Created
Mon, Nov 4, 20:47

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

Event Timeline

c4science · Help