From 5d67faded31d32fdc4845a2d3572c213937d11bf Mon Sep 17 00:00:00 2001 From: Harry Moffat Date: Mon, 23 Mar 2009 23:25:30 +0000 Subject: [PATCH] Implemented fixes, so that a python demo problem would work as advertised. --- Cantera/src/numerics/Func1.cpp | 821 +++++++++-------- Cantera/src/numerics/Func1.h | 1557 ++++++++++++++++++++------------ 2 files changed, 1452 insertions(+), 926 deletions(-) diff --git a/Cantera/src/numerics/Func1.cpp b/Cantera/src/numerics/Func1.cpp index 6fcc8d53f..1f94b4b78 100644 --- a/Cantera/src/numerics/Func1.cpp +++ b/Cantera/src/numerics/Func1.cpp @@ -6,400 +6,515 @@ using namespace std; namespace Cantera { - static Func1* checkDupl(Func1& f) { - if (f.parent() != 0) - return &f.duplicate(); - else - return &f; - } - Func1& Sin1::derivative() const { - Func1* c = new Cos1(m_c); - Func1* r = &newTimesConstFunction(*c, m_c); + Func1::Func1() : + m_c(0.0), + m_f1(0), + m_f2(0), + m_parent(0) + { + } + + Func1::Func1(const Func1 &right) : + m_c(right.m_c), + m_f1(right.m_f1), + m_f2(right.m_f2), + m_parent(right.m_parent) + { + } + + Func1::~Func1() { + } + + Func1& Func1::operator=(const Func1 &right) { + if (&right == this) return *this; + m_c = right.m_c; + m_f1 = right.m_f1; + m_f2 = right.m_f2; + m_parent = right.m_parent; + return *this; + } + + Func1& Func1::duplicate() const { + Func1 *nfunc = new Func1(*this); + return *nfunc; + } + + int Func1::ID() const { + return 0; + } + + // Calls method eval to evaluate the function + doublereal Func1::operator()(doublereal t) const { + return eval(t); + } + + // Evaluate the function. + doublereal Func1::eval(doublereal t) const { + return 0.0; + } + + Func1& Func1::derivative() const { + cout << "derivative error... ERR: ID = " << ID() << endl; + cout << write("x") << endl; + return *(new Func1); + } + + bool Func1::isIdentical(Func1& other) const { + if ((ID() != other.ID()) || (m_c != other.m_c)) + return false; + if (m_f1) { + if (!other.m_f1) return false; + if (!m_f1->isIdentical(*other.m_f1)) return false; + } + if (m_f2) { + if (!other.m_f2) return false; + if (!m_f2->isIdentical(*other.m_f2)) return false; + } + return true; + } + + //! accessor function for the returned constant + doublereal Func1::c() const { + return m_c; + } + + // Function to set the storred constant + void Func1::setC(doublereal c) { + m_c = c; + } + + //! accessor function for m_f1 + Func1& Func1::func1() const { + return *m_f1; + } + + Func1& Func1::func2() const { + return *m_f2; + } + + int Func1::order() const { + return 3; + } + + + Func1& Func1::func1_dup() const { + return m_f1->duplicate(); + } + + + Func1& Func1::func2_dup() const { + return m_f2->duplicate(); + } + + + Func1* Func1::parent() const { + return m_parent; + } + + + void Func1::setParent(Func1* p) { + m_parent = p; + } + + /*****************************************************************************/ + + string Sin1::write(string arg) const { + string c = ""; + if (m_c != 1.0) c = fp2str(m_c); + return "\\sin(" + c + arg + ")"; + } + + Func1& Sin1::derivative() const { + Func1* c = new Cos1(m_c); + Func1* r = &newTimesConstFunction(*c, m_c); #ifdef DEBUG_FUNC - cout << "Sin1::derivative: \n"; - cout << "function = \'" + write("x") + "\'\n"; - cout << "derivative = \'" + r->write("x") + "\'\n"; + cout << "Sin1::derivative: \n"; + cout << "function = \'" + write("x") + "\'\n"; + cout << "derivative = \'" + r->write("x") + "\'\n"; #endif - return *r; + return *r; + } + /*****************************************************************************/ + + Func1& Cos1::derivative() const { + Func1* s = new Sin1(m_c); + Func1* r = &newTimesConstFunction(*s, -m_c); +#ifdef DEBUG_FUNC + cout << "Cos1::derivative: \n"; + cout << "function = \'" + write("x") + "\'\n"; + cout << "derivative = \'" + r->write("x") + "\'\n"; +#endif + return *r; + } + + std::string Cos1::write(std::string arg) const { + string c = ""; + if (m_c != 1.0) c = fp2str(m_c); + return "\\cos("+c+arg+")"; + } + + /**************************************************************************/ + + Func1& Exp1::derivative() const { + Func1* f = new Exp1(m_c); + if (m_c != 1.0) + return newTimesConstFunction(*f, m_c); + else + return *f; + } + + std::string Exp1::write(std::string arg) const { + string c = ""; + if (m_c != 1.0) c = fp2str(m_c); + return "\\exp("+c+arg+")"; + } + + /******************************************************************************/ + + Func1& Pow1::derivative() const { + Func1* r; + if (m_c == 0.0) { + r = new Const1(0.0); + } + else if (m_c == 1.0) { + r = new Const1(1.0); + } + else { + Func1* f = new Pow1(m_c - 1.0); + r = &newTimesConstFunction(*f, m_c); + } +#ifdef DEBUG_FUNC + cout << "Pow1::derivative: \n"; + cout << "function = \'" + write("x") + "\'\n"; + cout << "derivative = \'" + r->write("x") + "\'\n"; +#endif + return *r; + } + + string Func1::write(std::string arg) const { + return "("+arg+")"; + } + + + + + string Pow1::write(string arg) const { + //cout << "Pow1" << endl; + string c = ""; + if (m_c == 0.5) { + return "\\sqrt{" + arg + "}"; } - - Func1& Cos1::derivative() const { - Func1* s = new Sin1(m_c); - Func1* r = &newTimesConstFunction(*s, -m_c); -#ifdef DEBUG_FUNC - cout << "Cos1::derivative: \n"; - cout << "function = \'" + write("x") + "\'\n"; - cout << "derivative = \'" + r->write("x") + "\'\n"; -#endif - return *r; + if (m_c == -0.5) { + return "\\frac{1}{\\sqrt{" + arg + "}}"; + } + if (m_c != 1.0) { + c = fp2str(m_c); + return "\\left("+arg+"\\right)^{"+c+"}"; } - - Func1& Exp1::derivative() const { - Func1* f = new Exp1(m_c); - if (m_c != 1.0) - return newTimesConstFunction(*f, m_c); - else - return *f; + else { + return arg; } + } - Func1& Pow1::derivative() const { - Func1* r; - if (m_c == 0.0) { - r = new Const1(0.0); - } - else if (m_c == 1.0) { - r = new Const1(1.0); - } - else { - Func1* f = new Pow1(m_c - 1.0); - r = &newTimesConstFunction(*f, m_c); - } -#ifdef DEBUG_FUNC - cout << "Pow1::derivative: \n"; - cout << "function = \'" + write("x") + "\'\n"; - cout << "derivative = \'" + r->write("x") + "\'\n"; -#endif - return *r; - } - string Func1::write(std::string arg) const { - return "("+arg+")"; - } + string Const1::write(string arg) const { + //cout << "Const1" << endl; + string c = ""; + c = fp2str(m_c); + return c; + } - string Sin1::write(string arg) const { - string c = ""; - if (m_c != 1.0) c = fp2str(m_c); - return "\\sin("+c+arg+")"; - } + string Ratio1::write(string arg) const { + //cout << "Ratio1" << endl; + return "\\frac{" + m_f1->write(arg) + "}{" + + m_f2->write(arg) + "}"; + } - string Cos1::write(string arg) const { - string c = ""; - if (m_c != 1.0) c = fp2str(m_c); - return "\\cos("+c+arg+")"; - } + string Product1::write(string arg) const { + //cout << "Product1" << endl; + string s = m_f1->write(arg); + if (m_f1->order() < order()) s = "\\left(" + s + "\\right)"; + string s2 = m_f2->write(arg); + if (m_f2->order() < order()) s2 = "\\left(" + s2 + "\\right)"; + return s + " " + s2; + } - string Pow1::write(string arg) const { - //cout << "Pow1" << endl; - string c = ""; - if (m_c == 0.5) { - return "\\sqrt{" + arg + "}"; - } - if (m_c == -0.5) { - return "\\frac{1}{\\sqrt{" + arg + "}}"; - } - if (m_c != 1.0) { - c = fp2str(m_c); - return "\\left("+arg+"\\right)^{"+c+"}"; - } - else { - return arg; - } - } + string Sum1::write(string arg) const { + //cout << "Sum1" << endl; + string s1 = m_f1->write(arg); + string s2 = m_f2->write(arg); + if (s2[0] == '-') return s1 + " - " + s2.substr(1,s2.size()); + else return s1 + " + " + s2; + } - string Exp1::write(string arg) const { - string c = ""; - if (m_c != 1.0) c = fp2str(m_c); - return "\\exp("+c+arg+")"; - } + string Diff1::write(string arg) const { + //cout << "Diff1" << endl; + string s1 = m_f1->write(arg); + string s2 = m_f2->write(arg); + if (s2[0] == '-') return s1 + " + " + s2.substr(1,s2.size()); + else return s1 + " - " + s2; + } - string Const1::write(string arg) const { - //cout << "Const1" << endl; - string c = ""; - c = fp2str(m_c); - return c; - } + string Composite1::write(string arg) const { + //cout << "Composite1" << endl; + string g = m_f2->write(arg); + return m_f1->write(g); + } - string Ratio1::write(string arg) const { - //cout << "Ratio1" << endl; - return "\\frac{" + m_f1->write(arg) + "}{" - + m_f2->write(arg) + "}"; - } + string TimesConstant1::write(string arg) const { + //cout << "TimesConstant1" << endl; + string s = m_f1->write(arg); + if (m_f1->order() < order()) s = "\\left(" + s + "\\right)"; + if (m_c == 1.0) return s; + if (m_c == -1.0) return "-"+s; + char n = s[0]; + if (n >= '0' && n <= '9') + s = "\\left(" + s + "\\right)"; + return fp2str(m_c) + s; + } - string Product1::write(string arg) const { - //cout << "Product1" << endl; - string s = m_f1->write(arg); - if (m_f1->order() < order()) s = "\\left(" + s + "\\right)"; - string s2 = m_f2->write(arg); - if (m_f2->order() < order()) s2 = "\\left(" + s2 + "\\right)"; - return s + " " + s2; - } + string PlusConstant1::write(string arg) const { + //cout << "PlusConstant1" << endl; + if (m_c == 0.0) return m_f1->write(arg); + return m_f1->write(arg) + " + " + fp2str(m_c); + } - string Sum1::write(string arg) const { - //cout << "Sum1" << endl; - string s1 = m_f1->write(arg); - string s2 = m_f2->write(arg); - if (s2[0] == '-') return s1 + " - " + s2.substr(1,s2.size()); - else return s1 + " + " + s2; - } + doublereal Func1::isProportional(TimesConstant1& other) { + if (isIdentical(other.func1())) return other.c(); + return 0.0; + } + doublereal Func1::isProportional(Func1& other) { + if (isIdentical(other)) return 1.0; + else return 0.0; + } - string Diff1::write(string arg) const { - //cout << "Diff1" << endl; - string s1 = m_f1->write(arg); - string s2 = m_f2->write(arg); - if (s2[0] == '-') return s1 + " + " + s2.substr(1,s2.size()); - else return s1 + " - " + s2; - } + static bool isConstant(Func1& f) { + if (f.ID() == ConstFuncType) + return true; + else + return false; + } - string Composite1::write(string arg) const { - //cout << "Composite1" << endl; - string g = m_f2->write(arg); - return m_f1->write(g); - } + static bool isZero(Func1& f) { + if (f.ID() == ConstFuncType && f.c() == 0.0) + return true; + else + return false; + } - string TimesConstant1::write(string arg) const { - //cout << "TimesConstant1" << endl; - string s = m_f1->write(arg); - if (m_f1->order() < order()) s = "\\left(" + s + "\\right)"; - if (m_c == 1.0) return s; - if (m_c == -1.0) return "-"+s; - char n = s[0]; - if (n >= '0' && n <= '9') - s = "\\left(" + s + "\\right)"; - return fp2str(m_c) + s; - } + static bool isOne(Func1& f) { + if (f.ID() == ConstFuncType && f.c() == 1.0) + return true; + else + return false; + } - string PlusConstant1::write(string arg) const { - //cout << "PlusConstant1" << endl; - if (m_c == 0.0) return m_f1->write(arg); - return m_f1->write(arg) + " + " + fp2str(m_c); - } + static bool isTimesConst(Func1& f) { + if (f.ID() == TimesConstantFuncType) + return true; + else + return false; + } - doublereal Func1::isProportional(TimesConstant1& other) { - if (isIdentical(other.func1())) return other.c(); - return 0.0; - } - doublereal Func1::isProportional(Func1& other) { - if (isIdentical(other)) return 1.0; - else return 0.0; - } + static bool isExp(Func1& f) { + if (f.ID() == ExpFuncType) + return true; + else + return false; + } - static bool isConstant(Func1& f) { - if (f.ID() == ConstFuncType) - return true; - else - return false; - } + static bool isPow(Func1& f) { + if (f.ID() == PowFuncType) + return true; + else + return false; + } - static bool isZero(Func1& f) { - if (f.ID() == ConstFuncType && f.c() == 0.0) - return true; - else - return false; + Func1& newSumFunction(Func1& f1, Func1& f2) { + if (f1.isIdentical(f2)) + return newTimesConstFunction(f1, 2.0); + if (isZero(f1)) { + delete &f1; + return f2; } - - static bool isOne(Func1& f) { - if (f.ID() == ConstFuncType && f.c() == 1.0) - return true; - else - return false; + if (isZero(f2)) { + delete &f2; + return f1; } - - static bool isTimesConst(Func1& f) { - if (f.ID() == TimesConstantFuncType) - return true; - else - return false; + doublereal c = f1.isProportional(f2); + if (c != 0) { + if (c == -1.0) + return *(new Const1(0.0)); + else { + return newTimesConstFunction(f1, c + 1.0); + } } + return *(new Sum1(f1, f2)); + } - static bool isExp(Func1& f) { - if (f.ID() == ExpFuncType) - return true; - else - return false; + Func1& newDiffFunction(Func1& f1, Func1& f2) { + if (isZero(f2)) { + delete &f2; return f1; } - - static bool isPow(Func1& f) { - if (f.ID() == PowFuncType) - return true; - else - return false; + if (f1.isIdentical(f2)) { + delete &f1; delete &f2; + return *(new Const1(0.0)); } - - Func1& newSumFunction(Func1& f1, Func1& f2) { - if (f1.isIdentical(f2)) - return newTimesConstFunction(f1, 2.0); - if (isZero(f1)) { - delete &f1; - return f2; - } - if (isZero(f2)) { - delete &f2; - return f1; - } - doublereal c = f1.isProportional(f2); - if (c != 0) { - if (c == -1.0) - return *(new Const1(0.0)); - else { - return newTimesConstFunction(f1, c + 1.0); - } - } - return *(new Sum1(f1, f2)); + doublereal c = f1.isProportional(f2); + if (c != 0.0) { + if (c == 1.0) return *(new Const1(0.0)); + else return newTimesConstFunction(f1, 1.0 - c); } + return *(new Diff1(f1, f2)); + } - Func1& newDiffFunction(Func1& f1, Func1& f2) { - if (isZero(f2)) { - delete &f2; return f1; - } - if (f1.isIdentical(f2)) { - delete &f1; delete &f2; - return *(new Const1(0.0)); - } - doublereal c = f1.isProportional(f2); - if (c != 0.0) { - if (c == 1.0) return *(new Const1(0.0)); - else return newTimesConstFunction(f1, 1.0 - c); - } - return *(new Diff1(f1, f2)); + Func1& newProdFunction(Func1& f1, Func1& f2) { + if (isOne(f1)) { + delete &f1; return f2; + } + if (isOne(f2)) { + delete &f2; return f1; + } + if (isZero(f1) || isZero(f2)) { + delete &f1; delete &f2; + return *(new Const1(0.0)); + } + if (isConstant(f1) && isConstant(f2)) { + doublereal c1c2 = f1.c() * f2.c(); + delete &f1; delete &f2; + return *(new Const1(c1c2)); + } + if (isConstant(f1)) { + doublereal c = f1.c(); + delete &f1; + return newTimesConstFunction(f2, c); + } + if (isConstant(f2)) { + doublereal c = f2.c(); + delete &f2; + return newTimesConstFunction(f1, c); } - - Func1& newProdFunction(Func1& f1, Func1& f2) { - if (isOne(f1)) { - delete &f1; return f2; - } - if (isOne(f2)) { - delete &f2; return f1; - } - if (isZero(f1) || isZero(f2)) { - delete &f1; delete &f2; - return *(new Const1(0.0)); - } - if (isConstant(f1) && isConstant(f2)) { - doublereal c1c2 = f1.c() * f2.c(); - delete &f1; delete &f2; - return *(new Const1(c1c2)); - } - if (isConstant(f1)) { - doublereal c = f1.c(); - delete &f1; - return newTimesConstFunction(f2, c); - } - if (isConstant(f2)) { - doublereal c = f2.c(); - delete &f2; - return newTimesConstFunction(f1, c); - } - if (isPow(f1) && isPow(f2)) { - Func1& p = *(new Pow1(f1.c() + f2.c())); - delete &f1; delete &f2; - return p; - } - - if (isExp(f1) && isExp(f2)) { - Func1& p = *(new Exp1(f1.c() + f2.c())); - delete &f1; delete &f2; - return p; - } - - bool tc1 = isTimesConst(f1); - bool tc2 = isTimesConst(f2); - - if (tc1 || tc2) { - doublereal c1 = 1.0, c2 = 1.0; - Func1 *ff1 = 0, *ff2 = 0; - if (tc1) { - c1 = f1.c(); - ff1 = &f1.func1_dup(); - delete &f1; - } - else ff1 = &f1; - if (tc2) { - c2 = f2.c(); - ff2 = &f2.func1_dup(); - delete &f2; - } - else ff2 = &f2; - Func1& p = newProdFunction(*ff1, *ff2); - - if (c1*c2 != 1.0) { - return newTimesConstFunction(p, c1*c2); - } - else - return p; - } - else - return *(new Product1(f1, f2)); + if (isPow(f1) && isPow(f2)) { + Func1& p = *(new Pow1(f1.c() + f2.c())); + delete &f1; delete &f2; + return p; } - Func1& newRatioFunction(Func1& f1, Func1& f2) { - if (isOne(f2)) return f1; - if (isZero(f1)) return *(new Const1(0.0)); - if (f1.isIdentical(f2)) { - delete &f1; delete &f2; - return *(new Const1(1.0)); - } - if (f1.ID() == PowFuncType && f2.ID() == PowFuncType) { - return *(new Pow1(f1.c() - f2.c())); - } - if (f1.ID() == ExpFuncType && f2.ID() == ExpFuncType) { - return *(new Exp1(f1.c() - f2.c())); - } - return *(new Ratio1(f1, f2)); + if (isExp(f1) && isExp(f2)) { + Func1& p = *(new Exp1(f1.c() + f2.c())); + delete &f1; delete &f2; + return p; } + + bool tc1 = isTimesConst(f1); + bool tc2 = isTimesConst(f2); + + if (tc1 || tc2) { + doublereal c1 = 1.0, c2 = 1.0; + Func1 *ff1 = 0, *ff2 = 0; + if (tc1) { + c1 = f1.c(); + ff1 = &f1.func1_dup(); + delete &f1; + } + else ff1 = &f1; + if (tc2) { + c2 = f2.c(); + ff2 = &f2.func1_dup(); + delete &f2; + } + else ff2 = &f2; + Func1& p = newProdFunction(*ff1, *ff2); + + if (c1*c2 != 1.0) { + return newTimesConstFunction(p, c1*c2); + } + else + return p; + } + else + return *(new Product1(f1, f2)); + } + + Func1& newRatioFunction(Func1& f1, Func1& f2) { + if (isOne(f2)) return f1; + if (isZero(f1)) return *(new Const1(0.0)); + if (f1.isIdentical(f2)) { + delete &f1; delete &f2; + return *(new Const1(1.0)); + } + if (f1.ID() == PowFuncType && f2.ID() == PowFuncType) { + return *(new Pow1(f1.c() - f2.c())); + } + if (f1.ID() == ExpFuncType && f2.ID() == ExpFuncType) { + return *(new Exp1(f1.c() - f2.c())); + } + return *(new Ratio1(f1, f2)); + } - Func1& newCompositeFunction(Func1& f1, Func1& f2) { - //#ifdef DEBUG_FUNC - //cout << "creating new composite function." << endl; - //cout << "f1 = " << f1.write("x") << " " << f1.ID() << endl; - //cout << "f2 = " << f2.write("x") << " " << f2.ID() << endl; - //#endif - if (isZero(f1)) { - delete &f1; delete &f2; - return *(new Const1(0.0)); - } - if (isConstant(f1)) { - delete &f2; - return f1; - } - if (isPow(f1) && f1.c() == 1.0) { - delete &f1; - return f2; - } - if (isPow(f1) && f1.c() == 0.0) { - delete &f1; - delete &f2; - return *(new Const1(1.0)); - } - if (isPow(f1) && isPow(f2)) { - doublereal c1c2 = f1.c() * f2.c(); - delete &f1; - delete &f2; - return *(new Pow1(c1c2)); - } - return *(new Composite1(f1, f2)); + Func1& newCompositeFunction(Func1& f1, Func1& f2) { + //#ifdef DEBUG_FUNC + //cout << "creating new composite function." << endl; + //cout << "f1 = " << f1.write("x") << " " << f1.ID() << endl; + //cout << "f2 = " << f2.write("x") << " " << f2.ID() << endl; + //#endif + if (isZero(f1)) { + delete &f1; delete &f2; + return *(new Const1(0.0)); } + if (isConstant(f1)) { + delete &f2; + return f1; + } + if (isPow(f1) && f1.c() == 1.0) { + delete &f1; + return f2; + } + if (isPow(f1) && f1.c() == 0.0) { + delete &f1; + delete &f2; + return *(new Const1(1.0)); + } + if (isPow(f1) && isPow(f2)) { + doublereal c1c2 = f1.c() * f2.c(); + delete &f1; + delete &f2; + return *(new Pow1(c1c2)); + } + return *(new Composite1(f1, f2)); + } - Func1& newTimesConstFunction(Func1& f, doublereal c) { - if (c == 0.0) { - delete &f; - return *(new Const1(0.0)); - } - if (c == 1.0) { - return f; - } - if (f.ID() == TimesConstantFuncType) { - f.setC(f.c() * c); - return f; - } - return *(new TimesConstant1(f, c)); + Func1& newTimesConstFunction(Func1& f, doublereal c) { + if (c == 0.0) { + delete &f; + return *(new Const1(0.0)); } + if (c == 1.0) { + return f; + } + if (f.ID() == TimesConstantFuncType) { + f.setC(f.c() * c); + return f; + } + return *(new TimesConstant1(f, c)); + } - Func1& newPlusConstFunction(Func1& f, doublereal c) { - if (c == 0.0) { - return f; - } - if (isConstant(f)) { - doublereal cc = f.c() + c; - delete &f; - return *(new Const1(cc)); - } - if (f.ID() == PlusConstantFuncType) { - f.setC(f.c() + c); - return f; - } - return *(new PlusConstant1(f, c)); + Func1& newPlusConstFunction(Func1& f, doublereal c) { + if (c == 0.0) { + return f; } + if (isConstant(f)) { + doublereal cc = f.c() + c; + delete &f; + return *(new Const1(cc)); + } + if (f.ID() == PlusConstantFuncType) { + f.setC(f.c() + c); + return f; + } + return *(new PlusConstant1(f, c)); + } } diff --git a/Cantera/src/numerics/Func1.h b/Cantera/src/numerics/Func1.h index 246845d3e..fa3e6b23f 100644 --- a/Cantera/src/numerics/Func1.h +++ b/Cantera/src/numerics/Func1.h @@ -27,623 +27,1034 @@ using namespace std; namespace Cantera { - const int FourierFuncType = 1; - const int PolyFuncType = 2; - const int ArrheniusFuncType = 3; - const int GaussianFuncType = 4; - const int SumFuncType = 20; - const int DiffFuncType = 25; - const int ProdFuncType = 30; - const int RatioFuncType = 40; - const int PeriodicFuncType = 50; - const int CompositeFuncType = 60; - const int TimesConstantFuncType = 70; - const int PlusConstantFuncType = 80; - const int SinFuncType = 100; - const int CosFuncType = 102; - const int ExpFuncType = 104; - const int PowFuncType = 106; - const int ConstFuncType = 110; + const int FourierFuncType = 1; + const int PolyFuncType = 2; + const int ArrheniusFuncType = 3; + const int GaussianFuncType = 4; + const int SumFuncType = 20; + const int DiffFuncType = 25; + const int ProdFuncType = 30; + const int RatioFuncType = 40; + const int PeriodicFuncType = 50; + const int CompositeFuncType = 60; + const int TimesConstantFuncType = 70; + const int PlusConstantFuncType = 80; + const int SinFuncType = 100; + const int CosFuncType = 102; + const int ExpFuncType = 104; + const int PowFuncType = 106; + const int ConstFuncType = 110; - class Sin1; - class Cos1; - class Exp1; - class Pow1; - class TimesConstant1; + class Sin1; + class Cos1; + class Exp1; + class Pow1; + class TimesConstant1; - /** - * Base class for 'functor' classes that evaluate a function of - * one variable. + /** + * Base class for 'functor' classes that evaluate a function of + * one variable. + */ + class Func1 { + public: + Func1(); + + virtual ~Func1(); + + + Func1(const Func1 &right); + + Func1& operator=(const Func1 &right); + + //! Duplicate the current function. + /*! + * This duplicates the current function, returning a + * reference to the new malloced function. */ - class Func1 { - public: - Func1() : m_c(0.0), m_f1(0), m_f2(0), m_parent(0) {} - virtual ~Func1() {} - virtual int ID() const { return 0; } + virtual Func1& duplicate() const; - virtual Func1& duplicate() { cout << "DUPL ERR: ID = " << ID() << endl; - return *(new Func1);} + virtual int ID() const; - /// Calls method eval to evaluate the function - doublereal operator()(doublereal t) const { return eval(t); } + //! Calls method eval to evaluate the function + doublereal operator()(doublereal t) const; - /// Evaluate the function. - virtual doublereal eval(doublereal t) const { return 0.0; } + /// Evaluate the function. + virtual doublereal eval(doublereal t) const; - virtual Func1& derivative() const { - cout << "derivative error... ERR: ID = " << ID() << endl; - cout << write("x") << endl; - return *(new Func1); - } - - bool isIdentical(Func1& other) const { - if ((ID() != other.ID()) || (m_c != other.m_c)) - return false; - if (m_f1) { - if (!other.m_f1) return false; - if (!m_f1->isIdentical(*other.m_f1)) return false; - } - if (m_f2) { - if (!other.m_f2) return false; - if (!m_f2->isIdentical(*other.m_f2)) return false; - } - return true; - } - - virtual doublereal isProportional(TimesConstant1& other); - virtual doublereal isProportional(Func1& other); - - virtual std::string write(std::string arg) const; - - doublereal c() const { return m_c; } - void setC(doublereal c) { m_c = c; } - Func1& func1() { return *m_f1; } - Func1& func2() { return *m_f2; } - virtual int order() const { return 3; } - Func1& func1_dup() const { return m_f1->duplicate(); } - Func1& func2_dup() const { return m_f2->duplicate(); } - Func1* parent() { return m_parent; } - void setParent(Func1* p) { m_parent = p; } - - protected: - doublereal m_c; - Func1 *m_f1, *m_f2; - Func1* m_parent; - - private: - }; - - - Func1& newSumFunction(Func1& f1, Func1& f2); - Func1& newDiffFunction(Func1& f1, Func1& f2); - Func1& newProdFunction(Func1& f1, Func1& f2); - Func1& newRatioFunction(Func1& f1, Func1& f2); - Func1& newCompositeFunction(Func1& f1, Func1& f2); - Func1& newTimesConstFunction(Func1& f1, doublereal c); - Func1& newPlusConstFunction(Func1& f1, doublereal c); - - /// sin - class Sin1 : public Func1 { - public: - Sin1(doublereal omega = 1.0) { - m_c = omega; - } - virtual ~Sin1() {} - virtual std::string write(std::string arg) const; - virtual Func1& duplicate() { return *(new Sin1(m_c)); } - virtual int ID() const { return SinFuncType; } - virtual doublereal eval(doublereal t) const { - return sin(m_c*t); - } - virtual Func1& derivative() const; - - protected: - - }; - - /// cos - class Cos1 : public Func1 { - public: - Cos1(doublereal omega = 1.0) { - m_c = omega; - } - virtual ~Cos1() {} - virtual std::string write(std::string arg) const; - virtual Func1& duplicate() { return *(new Cos1(m_c)); } - virtual int ID() const { return CosFuncType; } - virtual doublereal eval(doublereal t) const { - return cos(m_c * t); - } - virtual Func1& derivative() const; - - protected: - }; - - /// exp - class Exp1 : public Func1 { - public: - Exp1(doublereal A = 1.0) {m_c = A;} - virtual ~Exp1() {} - virtual std::string write(std::string arg) const; - virtual int ID() const { return ExpFuncType; } - virtual Func1& duplicate() { return *(new Exp1(m_c)); } - virtual doublereal eval(doublereal t) const { - return exp(m_c*t); - } - - virtual Func1& derivative() const; - - protected: - - }; - - /// pow - class Pow1 : public Func1 { - public: - Pow1(doublereal n) {m_c = n;} - virtual ~Pow1() {} - virtual std::string write(std::string arg) const; - virtual int ID() const { return PowFuncType; } - virtual Func1& duplicate() { return *(new Pow1(m_c)); } - virtual doublereal eval(doublereal t) const { - return pow(t, m_c); - } - virtual Func1& derivative() const; - - protected: - - }; - - /** - * Constant. - */ - class Const1 : public Func1 { - public: - Const1(doublereal A) { - m_c = A; - } - virtual ~Const1() {} - virtual std::string write(std::string arg) const; - virtual int ID() const { return ConstFuncType; } - virtual doublereal eval(doublereal t) const { - return m_c; - } - virtual Func1& duplicate() { return *(new Const1(m_c)); } - virtual Func1& derivative() const { - Func1* z = new Const1(0.0); - return *z; - } - - protected: - }; - - - - /** - * Sum of two functions. + //! Creates a derivative to the current function + /*! + * This will malloc a derivative function and + * return a reference to the function. */ - class Sum1 : public Func1 { - public: - Sum1(Func1& f1, Func1& f2) { - m_f1 = &f1; - m_f2 = &f2; - m_f1->setParent(this); - m_f2->setParent(this); - } - virtual ~Sum1() { - delete m_f1; - delete m_f2; - } - virtual int ID() const { return SumFuncType; } - virtual doublereal eval(doublereal t) const { - return m_f1->eval(t) + m_f2->eval(t); - } - virtual Func1& duplicate() { - Func1& f1d = m_f1->duplicate(); - Func1& f2d = m_f2->duplicate(); - Func1& dup = newSumFunction(f1d, f2d); - return dup; - } - virtual Func1& derivative() const { - Func1& d1 = m_f1->derivative(); - Func1& d2 = m_f2->derivative(); - Func1& d = newSumFunction(d1, d2); - return d; - } - virtual int order() const { return 0; } - virtual std::string write(std::string arg) const; - protected: - }; + virtual Func1& derivative() const; - - /** - * Difference of two functions. + //! Routine to determine if two functions are the same. + /*! + * Two functions are the same if they are teh same function. + * This means that the ID and storred constant is the same. + * This means that the m_f1 and m_f2 are identical if they + * are non-null. */ - class Diff1 : public Func1 { - public: - Diff1(Func1& f1, Func1& f2) { - m_f1 = &f1; - m_f2 = &f2; - } - virtual ~Diff1() { - delete m_f1; - delete m_f2; - } - virtual int ID() const { return DiffFuncType; } - virtual doublereal eval(doublereal t) const { - return m_f1->eval(t) - m_f2->eval(t); - } - virtual Func1& duplicate() { - Func1& f1d = m_f1->duplicate(); - Func1& f2d = m_f2->duplicate(); - Func1& dup = newDiffFunction(f1d, f2d); - return dup; - } - virtual Func1& derivative() const { - Func1& d = newDiffFunction(m_f1->derivative(), m_f2->derivative()); - return d; - } - virtual int order() const { return 0; } - virtual std::string write(std::string arg) const; + bool isIdentical(Func1& other) const; - protected: - }; + virtual doublereal isProportional(TimesConstant1& other); + virtual doublereal isProportional(Func1& other); + + virtual std::string write(std::string arg) const; - /** - * Product of two functions. - */ - class Product1 : public Func1 { - public: - Product1(Func1& f1, Func1& f2) { - m_f1 = &f1; - m_f2 = &f2; - } + //! accessor function for the storred constant + doublereal c() const; - virtual ~Product1() { - delete m_f1; - delete m_f2; - } - virtual int ID() const { return ProdFuncType; } - virtual Func1& duplicate() { - Func1& f1d = m_f1->duplicate(); - Func1& f2d = m_f2->duplicate(); - Func1& dup = newProdFunction(f1d, f2d); - return dup; - } - virtual std::string write(std::string arg) const; - virtual doublereal eval(doublereal t) const { - return m_f1->eval(t) * m_f2->eval(t); - } - virtual Func1& derivative() const { - Func1& a1 = newProdFunction(m_f1->duplicate(), m_f2->derivative()); - Func1& a2 = newProdFunction(m_f2->duplicate(), m_f1->derivative()); - Func1& s = newSumFunction(a1, a2); - return s; - } - virtual int order() const { return 1; } + //! Function to set the storred constant + void setC(doublereal c); + + //! accessor function for m_f1 + Func1& func1() const; + + //! accessor function for m_f2 + Func1& func2() const; + + //! Return the order of the function, if it makes sense + virtual int order() const; + + + Func1& func1_dup() const; + + + Func1& func2_dup() const; + + Func1* parent() const; + + + void setParent(Func1* p); + + protected: + doublereal m_c; + Func1 *m_f1; + Func1 *m_f2; + Func1 *m_parent; + }; + + + Func1& newSumFunction(Func1& f1, Func1& f2); + Func1& newDiffFunction(Func1& f1, Func1& f2); + Func1& newProdFunction(Func1& f1, Func1& f2); + Func1& newRatioFunction(Func1& f1, Func1& f2); + Func1& newCompositeFunction(Func1& f1, Func1& f2); + Func1& newTimesConstFunction(Func1& f1, doublereal c); + Func1& newPlusConstFunction(Func1& f1, doublereal c); + + //! implements the sin() function + /*! + * The argument to sin() is in radians + */ + class Sin1 : public Func1 { + public: + + Sin1(doublereal omega = 1.0) : + Func1() + { + m_c = omega; + } + + virtual ~Sin1() {} + + Sin1(const Sin1 &b) : + Func1(b) + { + } + + Sin1& operator=(const Sin1 &right) { + if (&right == this) return *this; + Func1::operator=(right); + return *this; + } + + virtual Func1& duplicate() const { + Sin1 *nfunc = new Sin1(*this); + return (Func1&) *nfunc; + } + + virtual std::string write(std::string arg) const; + + virtual int ID() const { return SinFuncType; } + + virtual doublereal eval(doublereal t) const { + return sin(m_c*t); + } + + virtual Func1& derivative() const; + }; + + /// cos + class Cos1 : public Func1 { + public: + Cos1(doublereal omega = 1.0) : + Func1() + { + m_c = omega; + } + virtual ~Cos1() {} + Cos1(const Cos1 &b) : + Func1(b) + { + } + + Cos1& operator=(const Cos1 &right) { + if (&right == this) return *this; + Func1::operator=(right); + return *this; + } + + virtual Func1& duplicate() const { + Cos1 *nfunc = new Cos1(*this); + return (Func1&) *nfunc; + } + virtual std::string write(std::string arg) const; + virtual int ID() const { return CosFuncType; } + virtual doublereal eval(doublereal t) const { + return cos(m_c * t); + } + virtual Func1& derivative() const; + + protected: + }; + + /// exp + class Exp1 : public Func1 { + public: + Exp1(doublereal A = 1.0) : + Func1() + { + m_c = A; + } + virtual ~Exp1() { + } + Exp1(const Exp1 &b) : + Func1(b) + { + } + Exp1& operator=(const Exp1 &right) { + if (&right == this) return *this; + Func1::operator=(right); + return *this; + } + virtual std::string write(std::string arg) const; + virtual int ID() const { return ExpFuncType; } + virtual Func1& duplicate() const { return *(new Exp1(m_c)); } + virtual doublereal eval(doublereal t) const { + return exp(m_c*t); + } + + virtual Func1& derivative() const; + + protected: + + }; + + /// pow + class Pow1 : public Func1 { + public: + Pow1(doublereal n) : + Func1() + { + m_c = n; + } + virtual ~Pow1() {} + Pow1(const Pow1 &b) : + Func1(b) + { + } + Pow1& operator=(const Pow1 &right) { + if (&right == this) return *this; + Func1::operator=(right); + return *this; + } + virtual std::string write(std::string arg) const; + virtual int ID() const { return PowFuncType; } + virtual Func1& duplicate() const { return *(new Pow1(m_c)); } + virtual doublereal eval(doublereal t) const { + return pow(t, m_c); + } + virtual Func1& derivative() const; + + protected: + + }; + + /** + * Constant. + */ + class Const1 : public Func1 { + public: + Const1(doublereal A) : + Func1() + { + m_c = A; + } + virtual ~Const1() { + } + + Const1(const Const1 &b) : + Func1(b) + { + } + + Const1& operator=(const Const1 &right) { + if (&right == this) return *this; + Func1::operator=(right); + return *this; + } + + virtual std::string write(std::string arg) const; + virtual int ID() const { return ConstFuncType; } + virtual doublereal eval(doublereal t) const { + return m_c; + } + virtual Func1& duplicate() const { return *(new Const1(m_c)); } + + virtual Func1& derivative() const { + Func1* z = new Const1(0.0); + return *z; + } + }; + + + + /** + * Sum of two functions. + */ + class Sum1 : public Func1 { + public: + Sum1(Func1& f1, Func1& f2) : + Func1() + { + m_f1 = &f1; + m_f2 = &f2; + m_f1->setParent(this); + m_f2->setParent(this); + } + + virtual ~Sum1() { + delete m_f1; + delete m_f2; + } + + Sum1(const Sum1 &b) : + Func1(b) + { + *this = Sum1::operator=(b); + } + + Sum1& operator=(const Sum1 &right) { + if (&right == this) return *this; + Func1::operator=(right); + m_f1 = &(m_f1->duplicate()); + m_f2 = &(m_f2->duplicate()); + m_f1->setParent(this); + m_f2->setParent(this); + m_parent = 0; + return *this; + } + + virtual int ID() const { + return SumFuncType; + } + + virtual doublereal eval(doublereal t) const { + return m_f1->eval(t) + m_f2->eval(t); + } + + virtual Func1& duplicate() const { + Func1& f1d = m_f1->duplicate(); + Func1& f2d = m_f2->duplicate(); + Func1& dup = newSumFunction(f1d, f2d); + return dup; + } + + virtual Func1& derivative() const { + Func1& d1 = m_f1->derivative(); + Func1& d2 = m_f2->derivative(); + Func1& d = newSumFunction(d1, d2); + return d; + } + virtual int order() const { return 0; } + + virtual std::string write(std::string arg) const; + }; + + + /** + * Difference of two functions. + */ + class Diff1 : public Func1 { + public: + Diff1(Func1& f1, Func1& f2) { + m_f1 = &f1; + m_f2 = &f2; + m_f1->setParent(this); + m_f2->setParent(this); + } + + virtual ~Diff1() { + delete m_f1; + delete m_f2; + } + + Diff1(const Diff1 &b) : + Func1(b) + { + *this = Diff1::operator=(b); + } + + Diff1& operator=(const Diff1 &right) { + if (&right == this) return *this; + Func1::operator=(right); + m_f1 = &(m_f1->duplicate()); + m_f2 = &(m_f2->duplicate()); + m_f1->setParent(this); + m_f2->setParent(this); + m_parent = 0; + return *this; + } + + virtual int ID() const { + return DiffFuncType; + } + + virtual doublereal eval(doublereal t) const { + return m_f1->eval(t) - m_f2->eval(t); + } + + virtual Func1& duplicate() const { + Func1& f1d = m_f1->duplicate(); + Func1& f2d = m_f2->duplicate(); + Func1& dup = newDiffFunction(f1d, f2d); + return dup; + } + virtual Func1& derivative() const { + Func1& d = newDiffFunction(m_f1->derivative(), m_f2->derivative()); + return d; + } + virtual int order() const { return 0; } + + virtual std::string write(std::string arg) const; + + }; + + + /** + * Product of two functions. + */ + class Product1 : public Func1 { + public: + Product1(Func1& f1, Func1& f2) : + Func1() + { + m_f1 = &f1; + m_f2 = &f2; + m_f1->setParent(this); + m_f2->setParent(this); + } + + virtual ~Product1() { + delete m_f1; + delete m_f2; + } + + Product1(const Product1 &b) : + Func1(b) + { + *this = Product1::operator=(b); + } + + Product1& operator=(const Product1 &right) { + if (&right == this) return *this; + Func1::operator=(right); + m_f1 = &(m_f1->duplicate()); + m_f2 = &(m_f2->duplicate()); + m_f1->setParent(this); + m_f2->setParent(this); + m_parent = 0; + return *this; + } + + virtual int ID() const { + return ProdFuncType; + } + + virtual Func1& duplicate() const { + Func1& f1d = m_f1->duplicate(); + Func1& f2d = m_f2->duplicate(); + Func1& dup = newProdFunction(f1d, f2d); + return dup; + } + + virtual std::string write(std::string arg) const; + + virtual doublereal eval(doublereal t) const { + return m_f1->eval(t) * m_f2->eval(t); + } + + virtual Func1& derivative() const { + Func1& a1 = newProdFunction(m_f1->duplicate(), m_f2->derivative()); + Func1& a2 = newProdFunction(m_f2->duplicate(), m_f1->derivative()); + Func1& s = newSumFunction(a1, a2); + return s; + } + virtual int order() const { return 1; } - protected: - }; + protected: + }; - /** - * Product of two functions. - */ - class TimesConstant1 : public Func1 { - public: - TimesConstant1(Func1& f1, doublereal A) { - m_f1 = &f1; - m_c = A; - } + /** + * Product of two functions. + */ + class TimesConstant1 : public Func1 { + public: + TimesConstant1(Func1& f1, doublereal A) : + Func1() + { + m_f1 = &f1; + m_c = A; + m_f1->setParent(this); + } - virtual ~TimesConstant1() { - delete m_f1; - } - virtual int ID() const { return TimesConstantFuncType; } - virtual Func1& duplicate() { - Func1& f1 = m_f1->duplicate(); - Func1* dup = new TimesConstant1(f1, m_c); - return *dup; - } - virtual doublereal isProportional(TimesConstant1& other) { - if (func1().isIdentical(other.func1())) - return (other.c()/c()); - else - return 0.0; - } - virtual doublereal isProportional(Func1& other) { - if (func1().isIdentical(other)) return 1.0/c(); - else return 0.0; - } - virtual doublereal eval(doublereal t) const { - return m_f1->eval(t) * m_c; - } - virtual Func1& derivative() const { - Func1& f1d = m_f1->derivative(); - Func1* d = &newTimesConstFunction(f1d, m_c); - return *d; - } - virtual std::string write(std::string arg) const; - virtual int order() const { return 0; } - protected: - }; + virtual ~TimesConstant1() { + delete m_f1; + } - /** - * A function plus a constant. - */ - class PlusConstant1 : public Func1 { - public: - PlusConstant1(Func1& f1, doublereal A) { - m_f1 = &f1; - m_c = A; - } + TimesConstant1(const TimesConstant1 &b) : + Func1(b) + { + *this = TimesConstant1::operator=(b); + } - virtual ~PlusConstant1() { - delete m_f1; - } - virtual int ID() const { return PlusConstantFuncType; } - virtual Func1& duplicate() { - Func1& f1 = m_f1->duplicate(); - Func1* dup = new PlusConstant1(f1, m_c); - return *dup; - } + TimesConstant1& operator=(const TimesConstant1 &right) { + if (&right == this) return *this; + Func1::operator=(right); + m_f1 = &(m_f1->duplicate()); + m_f1->setParent(this); + m_parent = 0; + return *this; + } + virtual int ID() const { return TimesConstantFuncType; } - virtual doublereal eval(doublereal t) const { - return m_f1->eval(t) + m_c; - } - virtual Func1& derivative() const { - Func1& f1d = m_f1->derivative(); - return f1d; - } - virtual std::string write(std::string arg) const; - virtual int order() const { return 0; } + virtual Func1& duplicate() const { + Func1& f1 = m_f1->duplicate(); + Func1* dup = new TimesConstant1(f1, m_c); + return *dup; + } + + virtual doublereal isProportional(TimesConstant1& other) { + if (func1().isIdentical(other.func1())) + return (other.c()/c()); + else + return 0.0; + } + + virtual doublereal isProportional(Func1& other) { + if (func1().isIdentical(other)) return 1.0/c(); + else return 0.0; + } + + virtual doublereal eval(doublereal t) const { + return m_f1->eval(t) * m_c; + } + + virtual Func1& derivative() const { + Func1& f1d = m_f1->derivative(); + Func1* d = &newTimesConstFunction(f1d, m_c); + return *d; + } + + virtual std::string write(std::string arg) const; + + virtual int order() const { return 0; } + protected: + }; + + /** + * A function plus a constant. + */ + class PlusConstant1 : public Func1 { + public: + PlusConstant1(Func1& f1, doublereal A) : + Func1() + { + m_f1 = &f1; + m_c = A; + m_f1->setParent(this); + } + + virtual ~PlusConstant1() { + delete m_f1; + } + + PlusConstant1(const PlusConstant1 &b) : + Func1(b) + { + *this = PlusConstant1::operator=(b); + } + + PlusConstant1& operator=(const PlusConstant1 &right) { + if (&right == this) return *this; + Func1::operator=(right); + m_f1 = &(m_f1->duplicate()); + m_f1->setParent(this); + m_parent = 0; + return *this; + } + + virtual int ID() const { return PlusConstantFuncType; } + + virtual Func1& duplicate() const { + Func1& f1 = m_f1->duplicate(); + Func1* dup = new PlusConstant1(f1, m_c); + return *dup; + } + + virtual doublereal eval(doublereal t) const { + return m_f1->eval(t) + m_c; + } + virtual Func1& derivative() const { + Func1& f1d = m_f1->derivative(); + return f1d; + } + virtual std::string write(std::string arg) const; + + virtual int order() const { return 0; } - protected: - }; + protected: + }; - /** - * Ratio of two functions. - */ - class Ratio1 : public Func1 { - public: - Ratio1(Func1& f1, Func1& f2) { - m_f1 = &f1; - m_f2 = &f2; - } - virtual ~Ratio1() { - delete m_f1; - delete m_f2; - } - virtual int ID() const { return RatioFuncType; } - virtual doublereal eval(doublereal t) const { - return m_f1->eval(t) / m_f2->eval(t); - } - virtual Func1& duplicate() { - Func1& f1d = m_f1->duplicate(); - Func1& f2d = m_f2->duplicate(); - Func1& dup = newRatioFunction(f1d, f2d); - return dup; - } - virtual Func1& derivative() const { - Func1& a1 = newProdFunction(m_f1->derivative(), m_f2->duplicate()); - Func1& a2 = newProdFunction(m_f1->duplicate(), m_f2->derivative()); - Func1& s = newDiffFunction(a1, a2); - Func1& p = newProdFunction(m_f2->duplicate(), m_f2->duplicate()); - Func1& r = newRatioFunction(s, p); - return r; - } - virtual std::string write(std::string arg) const; - virtual int order() const { return 1; } + /** + * Ratio of two functions. + */ + class Ratio1 : public Func1 { + public: + Ratio1(Func1& f1, Func1& f2) : + Func1() + { + m_f1 = &f1; + m_f2 = &f2; + m_f1->setParent(this); + m_f2->setParent(this); + } - protected: - }; + virtual ~Ratio1() { + delete m_f1; + delete m_f2; + } - /** - * Composite function. - */ - class Composite1 : public Func1 { - public: - Composite1(Func1& f1, Func1& f2) { - m_f1 = &f1; - m_f2 = &f2; - } - virtual ~Composite1() { - delete m_f1; - delete m_f2; - } - virtual int ID() const { return CompositeFuncType; } - virtual doublereal eval(doublereal t) const { - return m_f1->eval( m_f2->eval(t) ); - } - virtual Func1& duplicate() { - Func1& f1d = m_f1->duplicate(); - Func1& f2d = m_f2->duplicate(); - Func1& dup = newCompositeFunction(f1d, f2d); - return dup; - } - virtual Func1& derivative() const { - Func1* d1 = &m_f1->derivative(); + Ratio1(const Ratio1 &b) : + Func1(b) + { + *this = Ratio1::operator=(b); + } - Func1* d3 = &newCompositeFunction(*d1, m_f2->duplicate()); - Func1* d2 = &m_f2->derivative(); - Func1* p = &newProdFunction(*d3, *d2); + Ratio1& operator=(const Ratio1 &right) { + if (&right == this) return *this; + Func1::operator=(right); + m_f1 = &(m_f1->duplicate()); + m_f2 = &(m_f2->duplicate()); + m_f1->setParent(this); + m_f2->setParent(this); + m_parent = 0; + return *this; + } + + virtual int ID() const { return RatioFuncType; } + + virtual doublereal eval(doublereal t) const { + return m_f1->eval(t) / m_f2->eval(t); + } + + virtual Func1& duplicate() const { + Func1& f1d = m_f1->duplicate(); + Func1& f2d = m_f2->duplicate(); + Func1& dup = newRatioFunction(f1d, f2d); + return dup; + } + + virtual Func1& derivative() const { + Func1& a1 = newProdFunction(m_f1->derivative(), m_f2->duplicate()); + Func1& a2 = newProdFunction(m_f1->duplicate(), m_f2->derivative()); + Func1& s = newDiffFunction(a1, a2); + Func1& p = newProdFunction(m_f2->duplicate(), m_f2->duplicate()); + Func1& r = newRatioFunction(s, p); + return r; + } + + virtual std::string write(std::string arg) const; + + virtual int order() const { return 1; } + + protected: + }; + + /** + * Composite function. + */ + class Composite1 : public Func1 { + public: + Composite1(Func1& f1, Func1& f2) : + Func1() + { + m_f1 = &f1; + m_f2 = &f2; + m_f1->setParent(this); + m_f2->setParent(this); + } + + virtual ~Composite1() { + delete m_f1; + delete m_f2; + } + + Composite1(const Composite1 &b) : + Func1(b) + { + *this = Composite1::operator=(b); + } + + Composite1& operator=(const Composite1 &right) { + if (&right == this) return *this; + Func1::operator=(right); + m_f1 = &(m_f1->duplicate()); + m_f2 = &(m_f2->duplicate()); + m_f1->setParent(this); + m_f2->setParent(this); + m_parent = 0; + return *this; + } + + virtual int ID() const { return CompositeFuncType; } + + virtual doublereal eval(doublereal t) const { + return m_f1->eval( m_f2->eval(t) ); + } + + virtual Func1& duplicate() const { + Func1& f1d = m_f1->duplicate(); + Func1& f2d = m_f2->duplicate(); + Func1& dup = newCompositeFunction(f1d, f2d); + return dup; + } + + virtual Func1& derivative() const { + Func1* d1 = &m_f1->derivative(); + + Func1* d3 = &newCompositeFunction(*d1, m_f2->duplicate()); + Func1* d2 = &m_f2->derivative(); + Func1* p = &newProdFunction(*d3, *d2); #ifdef DEBUG_FUNC - cout << "Composite1::derivative: \n"; - cout << "f1 = " << m_f1->write("x") << endl; - cout << "f2 = " << m_f2->write("x") << endl; - cout << "d1 = " << d1 << " " << d1->write("x") << endl; - cout << "d3 = " << d3->write("x") << endl; - cout << "d2 = " << d2->write("x") << endl; - cout << "function = \'" + write("x") + "\'\n"; - cout << "derivative = \'" + p->write("x") + "\'\n"; + cout << "Composite1::derivative: \n"; + cout << "f1 = " << m_f1->write("x") << endl; + cout << "f2 = " << m_f2->write("x") << endl; + cout << "d1 = " << d1 << " " << d1->write("x") << endl; + cout << "d3 = " << d3->write("x") << endl; + cout << "d2 = " << d2->write("x") << endl; + cout << "function = \'" + write("x") + "\'\n"; + cout << "derivative = \'" + p->write("x") + "\'\n"; #endif - return *p; - } - virtual std::string write(std::string arg) const; - virtual int order() const { return 2; } - protected: - }; + return *p; + } - // - // The functors below are the old-style ones. They still work, - // but can't do derivatives. - // + virtual std::string write(std::string arg) const; - /** - * A Gaussian. - * \f[ - * f(t) = A e^{-[(t - t_0)/\tau]^2} - * \f] - * where \f[ \tau = \frac{fwhm}{2\sqrt{\ln 2}} \f] - * @param A peak value - * @param t0 offset - * @param fwhm full width at half max - */ - class Gaussian : public Func1 { - public: - Gaussian(double A, double t0, double fwhm) { - m_A = A; - m_t0 = t0; - m_tau = fwhm/(2.0*std::sqrt(std::log(2.0))); - } - virtual ~Gaussian() {} - virtual doublereal eval(doublereal t) const { - doublereal x = (t - m_t0)/m_tau; - return m_A*std::exp(-x*x); - } - protected: - doublereal m_A, m_t0, m_tau; - private: - }; + virtual int order() const { return 2; } + + protected: + }; + + // + // The functors below are the old-style ones. They still work, + // but can't do derivatives. + // + + /** + * A Gaussian. + * \f[ + * f(t) = A e^{-[(t - t_0)/\tau]^2} + * \f] + * where \f[ \tau = \frac{fwhm}{2\sqrt{\ln 2}} \f] + * @param A peak value + * @param t0 offset + * @param fwhm full width at half max + */ + class Gaussian : public Func1 { + public: + Gaussian(double A, double t0, double fwhm) : + Func1() + { + m_A = A; + m_t0 = t0; + m_tau = fwhm/(2.0*std::sqrt(std::log(2.0))); + } + + virtual ~Gaussian() {} + + Gaussian(const Gaussian &b) : + Func1(b) + { + *this = Gaussian::operator=(b); + } + + Gaussian& operator=(const Gaussian &right) { + if (&right == this) return *this; + Func1::operator=(right); + m_A = right.m_A; + m_t0 = right.m_t0; + m_tau = right.m_tau; + m_parent = 0; + return *this; + } + + virtual Func1& duplicate() const { + Gaussian *np = new Gaussian(*this); + return *((Func1 *)np); + } + + virtual doublereal eval(doublereal t) const { + doublereal x = (t - m_t0)/m_tau; + return m_A * std::exp(-x*x); + } + + protected: + doublereal m_A, m_t0, m_tau; + private: + }; - /** - * Polynomial of degree n. - */ - class Poly1 : public Func1 { - public: - Poly1(int n, doublereal* c) { - m_n = n+1; - m_c.resize(n+1); - std::copy(c, c+m_n, m_c.begin()); - } - virtual ~Poly1() {} + /** + * Polynomial of degree n. + */ + class Poly1 : public Func1 { + public: + Poly1(int n, doublereal* c) : + Func1() + { + m_n = n+1; + m_cpoly.resize(n+1); + std::copy(c, c+m_n, m_cpoly.begin()); + } - virtual doublereal eval(doublereal t) const { - int n; - doublereal r = m_c[m_n-1]; - for (n = 1; n < m_n; n++) { - r *= t; - r += m_c[m_n - n - 1]; - } - return r; - } + virtual ~Poly1() { + } - protected: - int m_n; - vector_fp m_c; - }; + Poly1(const Poly1 &b) : + Func1(b) + { + *this = Poly1::operator=(b); + } + + Poly1& operator=(const Poly1 &right) { + if (&right == this) return *this; + Func1::operator=(right); + m_cpoly = right.m_cpoly; + m_n = right.m_n; + m_parent = 0; + return *this; + } - /** - * Fourier cosine/sine series. - * - * \f[ - * f(t) = \frac{A_0}{2} + - * \sum_{n=1}^N A_n \cos (n \omega t) + B_n \sin (n \omega t) - * \f] - */ - class Fourier1 : public Func1 { - public: - Fourier1(int n, doublereal omega, doublereal a0, - doublereal* a, doublereal* b) { - m_n = n; - m_omega = omega; - m_a0_2 = 0.5*a0; - m_ccos.resize(n); - m_csin.resize(n); - std::copy(a, a+n, m_ccos.begin()); - std::copy(b, b+n, m_csin.begin()); - } - virtual ~Fourier1() {} + virtual Func1& duplicate() const { + Poly1 *np = new Poly1(*this); + return *((Func1 *)np); + } - virtual doublereal eval(doublereal t) const { - int n, nn; - doublereal sum = m_a0_2; - for (n = 0; n < m_n; n++) { - nn = n + 1; - sum += m_ccos[n]*std::cos(m_omega*nn*t) - + m_csin[n]*std::sin(m_omega*nn*t); - } - return sum; - } + virtual doublereal eval(doublereal t) const { + int n; + doublereal r = m_cpoly[m_n-1]; + for (n = 1; n < m_n; n++) { + r *= t; + r += m_cpoly[m_n - n - 1]; + } + return r; + } - protected: - int m_n; - doublereal m_omega, m_a0_2; - vector_fp m_ccos, m_csin; - }; + protected: + int m_n; + vector_fp m_cpoly; + }; - /** - * Sum of Arrhenius terms. - * \f[ - * f(T) = \sum_{n=1}^N A_n T^b_n \exp(-E_n/T) - * \f] - */ - class Arrhenius1 : public Func1 { - public: - Arrhenius1(int n, doublereal* c) { - m_n = n; - m_A.resize(n); - m_b.resize(n); - m_E.resize(n); - int loc; - for (int i = 0; i < n; i++) { - loc = 3*i; - m_A[i] = c[loc]; - m_b[i] = c[loc+1]; - m_E[i] = c[loc+2]; - } - } - virtual ~Arrhenius1() {} + /** + * Fourier cosine/sine series. + * + * \f[ + * f(t) = \frac{A_0}{2} + + * \sum_{n=1}^N A_n \cos (n \omega t) + B_n \sin (n \omega t) + * \f] + */ + class Fourier1 : public Func1 { + public: + Fourier1(int n, doublereal omega, doublereal a0, + doublereal* a, doublereal* b) : + Func1() + { + m_n = n; + m_omega = omega; + m_a0_2 = 0.5*a0; + m_ccos.resize(n); + m_csin.resize(n); + std::copy(a, a+n, m_ccos.begin()); + std::copy(b, b+n, m_csin.begin()); + } - virtual doublereal eval(doublereal t) const { - int n; - doublereal sum = 0.0; - for (n = 0; n < m_n; n++) { - sum += m_A[n]*std::pow(t,m_b[n])*std::exp(-m_E[n]/t); - } - return sum; - } + virtual ~Fourier1() {} - protected: - int m_n; - vector_fp m_A, m_b, m_E; - }; + Fourier1(const Fourier1 &b) : + Func1(b) + { + *this = Fourier1::operator=(b); + } - /** - * Periodic function. Takes any function and makes it - * periodic with period T. - */ - class Periodic1 : public Func1 { - public: - Periodic1(Func1& f, doublereal T) { - m_func = &f; - m_c = T; - } - virtual ~Periodic1() { delete m_func; } - virtual doublereal eval(doublereal t) const { - int np = int(t/m_c); - doublereal time = t - np*m_c; - return m_func->eval(time); - } - protected: - Func1* m_func; + Fourier1& operator=(const Fourier1 &right) { + if (&right == this) return *this; + Func1::operator=(right); + m_omega = right.m_omega; + m_a0_2 = right.m_a0_2; + m_ccos = right.m_ccos; + m_csin = right.m_csin; + m_n = right.m_n; + m_parent = 0; + return *this; + } - private: - }; + virtual Func1& duplicate() const { + Fourier1 *np = new Fourier1(*this); + return *((Func1 *)np); + } + + virtual doublereal eval(doublereal t) const { + int n, nn; + doublereal sum = m_a0_2; + for (n = 0; n < m_n; n++) { + nn = n + 1; + sum += m_ccos[n]*std::cos(m_omega*nn*t) + + m_csin[n]*std::sin(m_omega*nn*t); + } + return sum; + } + + protected: + int m_n; + doublereal m_omega, m_a0_2; + vector_fp m_ccos, m_csin; + }; + + + /** + * Sum of Arrhenius terms. + * \f[ + * f(T) = \sum_{n=1}^N A_n T^b_n \exp(-E_n/T) + * \f] + */ + class Arrhenius1 : public Func1 { + public: + Arrhenius1(int n, doublereal* c) : + Func1() + { + m_n = n; + m_A.resize(n); + m_b.resize(n); + m_E.resize(n); + int loc; + for (int i = 0; i < n; i++) { + loc = 3*i; + m_A[i] = c[loc]; + m_b[i] = c[loc+1]; + m_E[i] = c[loc+2]; + } + } + virtual ~Arrhenius1() { + } + + + Arrhenius1(const Arrhenius1 &b) : + Func1() + { + *this = Arrhenius1::operator=(b); + } + + Arrhenius1& operator=(const Arrhenius1 &right) { + if (&right == this) return *this; + Func1::operator=(right); + m_n = right.m_n; + m_A = right.m_A; + m_b = right.m_b; + m_E = right.m_E; + m_parent = 0; + return *this; + } + + virtual Func1& duplicate() const { + Arrhenius1 *np = new Arrhenius1(*this); + return *((Func1 *)np); + } + + virtual doublereal eval(doublereal t) const { + int n; + doublereal sum = 0.0; + for (n = 0; n < m_n; n++) { + sum += m_A[n]*std::pow(t,m_b[n])*std::exp(-m_E[n]/t); + } + return sum; + } + + protected: + int m_n; + vector_fp m_A, m_b, m_E; + }; + + /** + * Periodic function. Takes any function and makes it + * periodic with period T. + */ + class Periodic1 : public Func1 { + public: + Periodic1(Func1& f, doublereal T) : + Func1() + { + m_func = &f; + m_c = T; + } + + Periodic1(const Periodic1 &b) : + Func1() + { + *this = Periodic1::operator=(b); + } + + Periodic1& operator=(const Periodic1 &right) { + if (&right == this) return *this; + Func1::operator=(right); + m_func = &(right.m_func->duplicate()); + return *this; + } + + virtual Func1& duplicate() const { + Periodic1 *np = new Periodic1(*this); + return *((Func1 *)np); + } + + virtual ~Periodic1() { delete m_func; } + + virtual doublereal eval(doublereal t) const { + int np = int(t/m_c); + doublereal time = t - np*m_c; + return m_func->eval(time); + } + + protected: + Func1* m_func; + + private: + }; }