function.h 5.06 KB
``````#ifndef RTS_FUNCTION_H
#define RTS_FUNCTION_H

#include <string>
#include <cmath>

namespace rts{

//template class for a one-dimensional function
template <class Tx, class Ty>
class function
{
protected:

std::vector<Tx> X;
std::vector<Ty> Y;
Ty range[2];
Ty bounding[2];

//comparison function for searching lambda
static bool findCeiling(Tx ax, Tx bx){
return (ax > bx);
}

//process a string to extract a function (generally delimited by tabs, commas, spaces, etc.)
void process_string(std::string s){
std::stringstream ss(s);

Tx x;
Ty y;
std::string line;

while(!ss.eof()){

std::getline(ss, line);
if(line[0] == '#') continue;

std::stringstream lstream(line);

if(ss.eof()) break;
insert(x, y);			//insert the read value into the function
}
}

void update_range(Ty y){
//if the function is empty, just set the range to the given value
if(X.size() == 0){
range[0] = range[1] = y;
}

//otherwise compare the values and set them appropriately
if(y < range[0]) range[0] = y;
if(y > range[1]) range[1] = y;
}

public:
function(){
range[0] = range[1] = NAN;
bounding[0] = bounding[1] = 0;
}

//linear interpolation
Ty linear(Tx x) const
{
if(X.size() == 0)	return bounding[0];	//return zero if the function is empty
//return bounding values if x is outside the sample domain
if(x < X[0]) return bounding[0];
if(x > X.back()) return bounding[1];

unsigned int N = X.size();			//number of sample points

//declare an iterator
typedef typename std::vector< Tx >::iterator f_iter;	//declare an iterator
f_iter it;

//find the first X-coordinate that is greater than x
unsigned int i;
for(i = 0; i<N; i++){
if(X[i] > x)
break;
}
//i currently holds the ceiling

//if the wavelength is past the end of the list, return the last sample point
if(i == N) return Y[N];
//if the wavelength is before the beginning of the list, return the front
else if(i == 0) return Y[0];
//otherwise interpolate
else{
Tx xMax = X[i];
Tx xMin = X[i-1];

Tx a = (x - xMin) / (xMax - xMin);
Ty riMin = Y[i - 1];
Ty riMax = Y[i];
Ty interp = riMax * a + riMin * (1 - a);
return interp;
}
}

///add a data point to a function
void insert(Tx x, Ty y)
{
unsigned int N = X.size();		//number of sample points

update_range(y);			//update the range of the function

if(N == 0 || X[N-1] < x){
X.push_back(x);
Y.push_back(y);
return;
}

//declare an iterator and search for the x value
typename std::vector< Tx >::iterator it;
it = search(X.begin(), X.end(), &x, &x + 1, &function<Tx, Ty>::findCeiling);

//if the function value is past the end of the vector, add it to the back
if(*it == N){
X.push_back(x);
Y.push_back(y);
}
//otherwise add the value at the iterator position
else{
X.insert(it, x);
Y.insert(Y.begin() + *it, y);
}

}

Tx getX(unsigned int i) const{
return X[i];
}

Ty getY(unsigned int i) const{
return Y[i];
}

///get the number of data points in the function
unsigned int getN() const{
return X.size();
}

//look up an indexed component
Ty operator[](int i) const{
if(i <= X.size()){
std::cout<<"ERROR: accessing non-existing sample point in 'function'"<<std::endl;
exit(1);
}
return Y[i];
}

///linear interpolation
Tx operator()(Tx x) const{
return linear(x);
}

//add a constant to the function
function<Tx, Ty> operator+(Ty r) const{

function<Tx, Ty> result;

//if there are points in the function
if(X.size() > 0){
//add r to every point in f
for(unsigned int i=0; i<X.size(); i++){
result.X.push_back(X[i]);
result.Y.push_back(Y[i] + r);
}
}
else{
result = r;
}

return result;
}

function<Tx, Ty> & operator= (const Ty & rhs){
X.clear();
Y.clear();
if(rhs != 0)			//if the RHS is zero, just clear, otherwise add one value of RHS
insert(0, rhs);

return *this;
}

//output a string description of the function (used for console output and debugging)
std::string str(){
unsigned int N = X.size();

std::stringstream ss;
ss<<"# of samples: "<<N<<std::endl;

if(N == 0) return ss.str();

ss<<"domain: ["<<X[0]<<", "<<X[X.size()-1]<<"]"<<std::endl;
ss<<"range: ["<<range[0]<<", "<<range[1]<<"]"<<std::endl;

return ss.str();

}
std::string str_vals(){
std::stringstream ss;

unsigned int N = X.size();		//number of sample points
//output each function value
for(unsigned int i = 0; i<N; i++){
if(i != 0) ss<<std::endl;
ss<<X[i]<<", "<<Y[i];
}

return ss.str();

}

std::ifstream t(filename.c_str());
std::string str((std::istreambuf_iterator<char>(t)),
std::istreambuf_iterator<char>());

process_string(str);
}

};

}	//end namespace rts

#endif
``````