From 752ddd68b00c9821c1f73509b0f8c737d6500bb5 Mon Sep 17 00:00:00 2001 From: Ingmar Schoegl Date: Thu, 2 Jan 2025 08:45:59 -0700 Subject: [PATCH] [unittests] Add googletests for clib_experimental --- test/SConscript | 2 + test/clib_experimental/test_clib3.cpp | 304 +++++++++++++++++++++++ test/clib_experimental/test_ctfunc3.cpp | 305 ++++++++++++++++++++++++ 3 files changed, 611 insertions(+) create mode 100644 test/clib_experimental/test_clib3.cpp create mode 100644 test/clib_experimental/test_ctfunc3.cpp diff --git a/test/SConscript b/test/SConscript index 823ccc13a9..91709ce9d5 100644 --- a/test/SConscript +++ b/test/SConscript @@ -209,6 +209,8 @@ def addPythonTest(testname, subset): # Instantiate tests addTestProgram('clib', 'clib') +if localenv['clib_experimental']: + addTestProgram('clib_experimental', 'clib-experimental') addTestProgram('equil', 'equil') addTestProgram('general', 'general') addTestProgram('kinetics', 'kinetics') diff --git a/test/clib_experimental/test_clib3.cpp b/test/clib_experimental/test_clib3.cpp new file mode 100644 index 0000000000..5cb1672799 --- /dev/null +++ b/test/clib_experimental/test_clib3.cpp @@ -0,0 +1,304 @@ +#include +#include +#include + +#include "cantera/core.h" +#include "cantera/clib/clib_defs.h" +#include "cantera/clib_experimental/ct3.h" +#include "cantera/clib_experimental/ctsol3.h" +#include "cantera/clib_experimental/ctthermo3.h" +#include "cantera/clib_experimental/ctkin3.h" +#include "cantera/clib_experimental/cttrans3.h" + +using namespace Cantera; +using ::testing::HasSubstr; + +string reportError() +{ + vector output_buf; + int buflen = ct3_getCanteraError(0, output_buf.data()); + output_buf.resize(buflen); + ct3_getCanteraError(buflen, output_buf.data()); + return string(output_buf.data()); +} + +TEST(ct3, cabinet_exceptions) +{ + sol3_newSolution("h2o2.yaml", "ohmech", "default"); + sol3_name(999, 0, 0); + + string err = reportError(); + EXPECT_THAT(err, HasSubstr("Index 999 out of range.")); + + sol3_thermo(998); + err = reportError(); + EXPECT_THAT(err, HasSubstr("Index 998 out of range.")); + + int ret = sol3_del(997); + ASSERT_EQ(ret, -1); + err = reportError(); + EXPECT_THAT(err, HasSubstr("Index 997 out of range.")); + + int ref = sol3_newSolution("h2o2.yaml", "ohmech", "default"); + sol3_del(ref); + int thermo = sol3_thermo(ref); + EXPECT_EQ(thermo, -2); + err = reportError(); + EXPECT_THAT(err, HasSubstr("has been deleted.")); + + ct3_resetStorage(); + ret = sol3_del(0); + ASSERT_EQ(ret, -1); + err = reportError(); + EXPECT_THAT(err, HasSubstr("Index 0 out of range.")); +} + +TEST(ct3, new_solution) +{ + ct3_resetStorage(); + + string name = "ohmech"; + int ref = sol3_newSolution("h2o2.yaml", name.c_str(), "default"); + ASSERT_EQ(ref, 0); + + ASSERT_EQ(sol3_cabinetSize(), 1); + ASSERT_EQ(thermo3_cabinetSize(), 1); + ASSERT_EQ(kin3_cabinetSize(), 1); + + int buflen = sol3_name(ref, 0, 0); // includes \0 + ASSERT_EQ(buflen, int(name.size() + 1)); + + int thermo = sol3_thermo(ref); + ASSERT_EQ(thermo3_parentHandle(thermo), ref); + + vector buf(buflen); + sol3_name(ref, buflen, buf.data()); + string solName(buf.data()); + ASSERT_EQ(solName, name); +} + +TEST(ct3, sol3_objects) +{ + ct3_resetStorage(); + + int ref = sol3_newSolution("gri30.yaml", "gri30", "none"); + ASSERT_EQ(ref, 0); + ASSERT_EQ(thermo3_cabinetSize(), 1); // one ThermoPhase object + + int ref2 = sol3_newSolution("h2o2.yaml", "ohmech", "default"); + ASSERT_EQ(ref2, 1); + ASSERT_EQ(thermo3_cabinetSize(), 2); // two ThermoPhase objects + + int thermo = sol3_thermo(ref); + ASSERT_EQ(thermo3_parentHandle(thermo), ref); + + int thermo2 = sol3_thermo(ref2); + ASSERT_EQ(thermo2, 1); // references stored object with index '1' + ASSERT_EQ(thermo3_nSpecies(thermo2), 10u); + ASSERT_EQ(thermo3_parentHandle(thermo2), ref2); + + int kin = sol3_kinetics(ref); + + int kin2 = sol3_kinetics(ref2); + ASSERT_EQ(kin2, 1); + ASSERT_EQ(kin3_nReactions(kin2), 29u); + ASSERT_EQ(kin3_parentHandle(kin2), ref2); + ASSERT_EQ(kin3_parentHandle(kin), ref); + + int trans = sol3_transport(ref); + ASSERT_EQ(trans3_parentHandle(trans), ref); + + int trans2 = sol3_transport(ref2); + ASSERT_EQ(trans2, 1); + int buflen = trans3_transportModel(trans2, 0, 0); + vector buf(buflen); + trans3_transportModel(trans2, buflen, buf.data()); + string trName(buf.data()); + ASSERT_EQ(trName, "mixture-averaged"); + ASSERT_EQ(trans3_parentHandle(trans2), ref2); + + sol3_del(ref2); + int nsp = thermo3_nSpecies(thermo2); + ASSERT_EQ(nsp, ERR); + string err = reportError(); + EXPECT_THAT(err, HasSubstr("has been deleted.")); + + nsp = thermo3_nSpecies(thermo2); + ASSERT_EQ(nsp, ERR); + err = reportError(); + EXPECT_THAT(err, HasSubstr("has been deleted.")); + + trans2 = sol3_setTransportModel(ref, "mixture-averaged"); + ASSERT_EQ(trans2, 2); + buflen = trans3_transportModel(trans2, 0, 0); + buf.resize(buflen); + trans3_transportModel(trans2, buflen, buf.data()); + trName = buf.data(); + ASSERT_EQ(trName, "mixture-averaged"); +} + +TEST(ct3, new_interface) +{ + ct3_resetStorage(); + + int sol = sol3_newSolution("ptcombust.yaml", "gas", "none"); + ASSERT_EQ(sol, 0); + + vector adj{sol}; + int surf = sol3_newInterface("ptcombust.yaml", "Pt_surf", 1, adj.data()); + ASSERT_EQ(surf, 1); + + int ph_surf = sol3_thermo(surf); + int buflen = sol3_name(ph_surf, 0, 0) + 1; // include \0 + vector buf(buflen); + sol3_name(ph_surf, buflen, buf.data()); + string solName(buf.data()); + ASSERT_EQ(solName, "Pt_surf"); + + int kin3_surf = sol3_kinetics(surf); + buflen = kin3_kineticsType(kin3_surf, 0, 0) + 1; // include \0 + buf.resize(buflen); + kin3_kineticsType(ph_surf, buflen, buf.data()); + string kinType(buf.data()); + ASSERT_EQ(kinType, "surface"); +} + +TEST(ct3, new_interface_auto) +{ + ct3_resetStorage(); + + vector adj; + int surf = sol3_newInterface("ptcombust.yaml", "Pt_surf", 0, adj.data()); + ASSERT_EQ(surf, 0); + + ASSERT_EQ(sol3_nAdjacent(surf), 1u); + int gas = sol3_adjacent(surf, 0); + ASSERT_EQ(gas, 1); + + int buflen = sol3_name(gas, 0, 0) + 1; // include \0 + vector buf(buflen); + sol3_name(gas, buflen, buf.data()); + string solName(buf.data()); + ASSERT_EQ(solName, "gas"); +} + +TEST(ct3, thermo) +{ + int ret; + int sol = sol3_newSolution("gri30.yaml", "gri30", "none"); + int thermo = sol3_thermo(sol); + ASSERT_GE(thermo, 0); + size_t nsp = thermo3_nSpecies(thermo); + ASSERT_EQ(nsp, 53u); + + ret = thermo3_setTemperature(thermo, 500); + ASSERT_EQ(ret, 0); + ret = thermo3_setPressure(thermo, 5 * 101325); + ASSERT_EQ(ret, 0); + ret = thermo3_setMoleFractionsByName(thermo, "CH4:1.0, O2:2.0, N2:7.52"); + ASSERT_EQ(ret, 0); + + ret = thermo3_equilibrate(thermo, "HP", "auto", 1e-9, 50000, 1000, 0); + ASSERT_EQ(ret, 0); + double T = thermo3_temperature(thermo); + ASSERT_GT(T, 2200); + ASSERT_LT(T, 2300); + + size_t ns = thermo3_nSpecies(thermo); + vector work(ns); + vector X(ns); + thermo3_getMoleFractions(thermo, ns, X.data()); + + thermo3_getPartialMolarEnthalpies(thermo, ns, work.data()); + double prod = std::inner_product(X.begin(), X.end(), work.begin(), 0.0); + ASSERT_NEAR(prod, thermo3_enthalpy_mole(thermo), 1e-6); + + thermo3_getPartialMolarEntropies(thermo, ns, work.data()); + prod = std::inner_product(X.begin(), X.end(), work.begin(), 0.0); + ASSERT_NEAR(prod, thermo3_entropy_mole(thermo), 1e-6); + + thermo3_getPartialMolarIntEnergies(thermo, ns, work.data()); + prod = std::inner_product(X.begin(), X.end(), work.begin(), 0.0); + ASSERT_NEAR(prod, thermo3_intEnergy_mole(thermo), 1e-6); + + thermo3_getPartialMolarCp(thermo, ns, work.data()); + prod = std::inner_product(X.begin(), X.end(), work.begin(), 0.0); + ASSERT_NEAR(prod, thermo3_cp_mole(thermo), 1e-6); + + thermo3_getPartialMolarVolumes(thermo, ns, work.data()); + prod = std::inner_product(X.begin(), X.end(), work.begin(), 0.0); + ASSERT_NEAR(prod, 1./thermo3_molarDensity(thermo), 1e-6); +} + +TEST(ct3, kinetics) +{ + int sol0 = sol3_newSolution("gri30.yaml", "gri30", "none"); + int thermo = sol3_thermo(sol0); + int kin = sol3_kinetics(sol0); + ASSERT_GE(kin, 0); + + size_t nr = kin3_nReactions(kin); + ASSERT_EQ(nr, 325u); + + thermo3_equilibrate(thermo, "HP", "auto", 1e-9, 50000, 1000, 0); + double T = thermo3_temperature(thermo); + thermo3_setTemperature(thermo, T - 200); + + auto sol = newSolution("gri30.yaml", "gri30", "none"); + auto phase = sol->thermo(); + auto kinetics = sol->kinetics(); + + phase->equilibrate("HP"); + ASSERT_NEAR(T, phase->temperature(), 1e-2); + phase->setTemperature(T - 200); + + vector c_ropf(nr); + kin3_getFwdRatesOfProgress(kin, 325, c_ropf.data()); + vector cpp_ropf(nr); + kinetics->getFwdRatesOfProgress(cpp_ropf.data()); + + for (size_t n = 0; n < nr; n++) { + ASSERT_NEAR(cpp_ropf[n], c_ropf[n], 1e-6); + } +} + +TEST(ct3, transport) +{ + int sol0 = sol3_newSolution("gri30.yaml", "gri30", "default"); + int thermo = sol3_thermo(sol0); + int tran = sol3_transport(sol0); + + size_t nsp = thermo3_nSpecies(thermo); + vector c_dkm(nsp); + int ret = trans3_getMixDiffCoeffs(tran, 53, c_dkm.data()); + ASSERT_EQ(ret, 0); + + vector cpp_dkm(nsp); + auto sol = newSolution("gri30.yaml", "gri30"); + auto transport = sol->transport(); + transport->getMixDiffCoeffs(cpp_dkm.data()); + + for (size_t n = 0; n < nsp; n++) { + ASSERT_NEAR(cpp_dkm[n], c_dkm[n], 1e-10); + } +} + + +int main(int argc, char** argv) +{ + printf("Running main() from test_clib3.cpp\n"); + testing::InitGoogleTest(&argc, argv); + make_deprecation_warnings_fatal(); + printStackTraceOnSegfault(); + Cantera::CanteraError::setStackTraceDepth(20); + /* commented code below is relevant to future PRs */ + // vector fileNames = {"gtest-freeflame.yaml", "gtest-freeflame.h5"}; + // for (const auto& fileName : fileNames) { + // if (std::ifstream(fileName).good()) { + // std::remove(fileName.c_str()); + // } + // } + int result = RUN_ALL_TESTS(); + appdelete(); + return result; +} diff --git a/test/clib_experimental/test_ctfunc3.cpp b/test/clib_experimental/test_ctfunc3.cpp new file mode 100644 index 0000000000..f46cd56309 --- /dev/null +++ b/test/clib_experimental/test_ctfunc3.cpp @@ -0,0 +1,305 @@ +#include +#include +#include + +#include "cantera/core.h" +#include "cantera/clib_experimental/ctfunc3.h" +#include "cantera/numerics/Func1Factory.h" + +using namespace Cantera; + + +TEST(ctfunc3, invalid) +{ + // exceptions return -1 + ASSERT_EQ(func13_newBasic("spam", 0.), -2); + ASSERT_EQ(func13_newAdvanced("eggs", 0, NULL), -2); +} + +TEST(ctfunc3, sin) +{ + double omega = 2.1; + int fcn = func13_newBasic("sin", omega); + ASSERT_GE(fcn, 0); + EXPECT_DOUBLE_EQ(func13_eval(fcn, 0.5), sin(omega * 0.5)); + + int dfcn = func13_newDerivative(fcn); + EXPECT_DOUBLE_EQ(func13_eval(dfcn, 0.5), omega * cos(omega * 0.5)); + + int buflen = func13_write(fcn, "x", 0, 0); + vector buf(buflen); + func13_write(fcn, "x", buflen, buf.data()); + string rep(buf.data()); + ASSERT_EQ(rep, "\\sin(2.1x)"); +} + +TEST(ctfunc3, cos) +{ + double omega = 2.; + int fcn = func13_newBasic("cos", omega); + ASSERT_GE(fcn, 0); + EXPECT_DOUBLE_EQ(func13_eval(fcn, 0.5), cos(omega * 0.5)); + + int dfcn = func13_newDerivative(fcn); + EXPECT_DOUBLE_EQ(func13_eval(dfcn, 0.5), -omega * sin(omega * 0.5)); +} + +TEST(ctfunc3, exp) +{ + double omega = 2.; + int fcn = func13_newBasic("exp", omega); + ASSERT_GE(fcn, 0); + EXPECT_DOUBLE_EQ(func13_eval(fcn, 0.5), exp(omega * 0.5)); + + int dfcn = func13_newDerivative(fcn); + EXPECT_DOUBLE_EQ(func13_eval(dfcn, 0.5), omega * exp(omega * 0.5)); +} + +TEST(ctfunc3, log) +{ + double omega = 2.; + int fcn = func13_newBasic("log", omega); + ASSERT_GE(fcn, 0); + EXPECT_DOUBLE_EQ(func13_eval(fcn, 1. / omega), 0.); + + int dfcn = func13_newDerivative(fcn); + EXPECT_DOUBLE_EQ(func13_eval(dfcn, .5), omega / .5); +} + +TEST(ctfunc3, pow) +{ + double exp = .5; + int fcn = func13_newBasic("pow", exp); + ASSERT_GE(fcn, 0); + EXPECT_DOUBLE_EQ(func13_eval(fcn, 0.5), pow(0.5, exp)); + + int dfcn = func13_newDerivative(fcn); + EXPECT_DOUBLE_EQ(func13_eval(dfcn, 0.5), exp * pow(0.5, exp - 1)); +} + +TEST(ctfunc3, constant) +{ + double a = .5; + int fcn = func13_newBasic("constant", a); + ASSERT_GE(fcn, 0); + EXPECT_DOUBLE_EQ(func13_eval(fcn, 0.5), a); + + int dfcn = func13_newDerivative(fcn); + EXPECT_DOUBLE_EQ(func13_eval(dfcn, .5), 0.); +} + +TEST(ctfunc3, tabulated_linear) +{ + vector params = {0., 1., 2., 1., 0., 1.}; + + int fcn = func13_newAdvanced("tabulated-linear", params.size(), params.data()); + ASSERT_GE(fcn, 0); + EXPECT_DOUBLE_EQ(func13_eval(fcn, 0.5), 0.5); + + // exceptions return -1 + EXPECT_EQ(func13_newAdvanced("tabulated-linear", 5, params.data()), -2); +} + +TEST(ctfunc3, tabulated_previous) +{ + vector params = {0., 1., 2., 1., 0., 1.}; + + int fcn = func13_newAdvanced("tabulated-previous", params.size(), params.data()); + ASSERT_GE(fcn, 0); + EXPECT_DOUBLE_EQ(func13_eval(fcn, 0.5), 1.); +} + +TEST(ctfunc3, poly) +{ + double a0 = .5; + double a1 = .25; + double a2 = .125; + vector params = {a2, a1, a0}; + int fcn = func13_newAdvanced("polynomial3", params.size(), params.data()); + ASSERT_GE(fcn, 0); + EXPECT_DOUBLE_EQ(func13_eval(fcn, .5), (a2 * .5 + a1) * .5 + a0); + + params = {1, 0, -2.2, 3.1}; + fcn = func13_newAdvanced("polynomial3", params.size(), params.data()); + int buflen = func13_write(fcn, "x", 0, 0); + vector buf(buflen); + func13_write(fcn, "x", buflen, buf.data()); + string rep(buf.data()); + ASSERT_EQ(rep, "x^3 - 2.2x + 3.1"); +} + +TEST(ctfunc3, Fourier) +{ + double a0 = .5; + double a1 = .25; + double b1 = .125; + double omega = 2.; + vector params = {a0, a1, omega, b1}; + int fcn = func13_newAdvanced("Fourier", params.size(), params.data()); + ASSERT_GE(fcn, 0); + EXPECT_DOUBLE_EQ( + func13_eval(fcn, .5), .5 * a0 + a1 * cos(omega * .5) + b1 * sin(omega * .5)); +} + +TEST(ctfunc3, Gaussian) +{ + double A = .5; + double t0 = .6; + double fwhm = .25; + vector params = {A, t0, fwhm}; + int fcn = func13_newAdvanced("Gaussian", params.size(), params.data()); + ASSERT_GE(fcn, 0); + double tau = fwhm / (2. * sqrt(log(2.))); + double x = - t0 / tau; + EXPECT_DOUBLE_EQ(func13_eval(fcn, 0.), A * exp(-x * x)); + + // exceptions return -1 + EXPECT_EQ(func13_newAdvanced("Gaussian", 2, params.data()), -2); +} + +TEST(ctfunc3, Arrhenius) +{ + double A = 38.7; + double b = 2.7; + double E = 2.619184e+07 / GasConstant; + vector params = {A, b, E}; + int fcn = func13_newAdvanced("Arrhenius", params.size(), params.data()); + ASSERT_GE(fcn, 0); + EXPECT_DOUBLE_EQ(func13_eval(fcn, 1000.), A * pow(1000., b) * exp(-E/1000.)); +} + +TEST(ctmath, invalid) +{ + // exceptions return -1 + int fcn0 = func13_newBasic("sin", 1.); + int fcn1 = func13_newBasic("cos", 1.); + ASSERT_EQ(func13_newCompound("foo", fcn0, fcn1), -2); + ASSERT_EQ(func13_newModified("bar", fcn0, 0.), -2); +} + +TEST(ctmath, sum) +{ + double omega = 2.; + int fcn0 = func13_newBasic("sin", omega); + int fcn1 = func13_newBasic("cos", omega); + int fcn = func13_newCompound("sum", fcn0, fcn1); + ASSERT_GE(fcn, 0); + EXPECT_DOUBLE_EQ(func13_eval(fcn, 0.5), sin(omega * 0.5) + cos(omega * 0.5)); + int fcn2 = func13_newSum(fcn0, fcn1); + EXPECT_DOUBLE_EQ(func13_eval(fcn, 0.5), func13_eval(fcn2, 0.5)); +} + +TEST(ctmath, diff) +{ + double omega = 2.; + int fcn0 = func13_newBasic("sin", omega); + int fcn1 = func13_newBasic("cos", omega); + int fcn = func13_newCompound("diff", fcn0, fcn1); + ASSERT_GE(fcn, 0); + EXPECT_DOUBLE_EQ(func13_eval(fcn, 0.5), sin(omega * 0.5) - cos(omega * 0.5)); + int fcn2 = func13_newDiff(fcn0, fcn1); + EXPECT_DOUBLE_EQ(func13_eval(fcn, 0.5), func13_eval(fcn2, 0.5)); +} + +TEST(ctmath, prod) +{ + double omega = 2.; + int fcn0 = func13_newBasic("sin", omega); + int fcn1 = func13_newBasic("cos", omega); + int fcn = func13_newCompound("product", fcn0, fcn1); + ASSERT_GE(fcn, 0); + EXPECT_DOUBLE_EQ(func13_eval(fcn, 0.5), sin(omega * 0.5) * cos(omega * 0.5)); + int fcn2 = func13_newProd(fcn0, fcn1); + EXPECT_DOUBLE_EQ(func13_eval(fcn, 0.5), func13_eval(fcn2, 0.5)); +} + +TEST(ctmath, ratio) +{ + double omega = 2.; + int fcn0 = func13_newBasic("sin", omega); + int fcn1 = func13_newBasic("cos", omega); + int fcn = func13_newCompound("ratio", fcn0, fcn1); + ASSERT_GE(fcn, 0); + EXPECT_DOUBLE_EQ(func13_eval(fcn, 0.5), sin(omega * 0.5) / cos(omega * 0.5)); + int fcn2 = func13_newRatio(fcn0, fcn1); + EXPECT_DOUBLE_EQ(func13_eval(fcn, 0.5), func13_eval(fcn2, 0.5)); +} + +TEST(ctmath, composite) +{ + double omega = 2.; + int fcn0 = func13_newBasic("sin", omega); + int fcn1 = func13_newBasic("cos", omega); + int fcn = func13_newCompound("composite", fcn0, fcn1); + ASSERT_GE(fcn, 0); + EXPECT_DOUBLE_EQ(func13_eval(fcn, 0.5), sin(omega * cos(omega * 0.5))); +} + +TEST(ctmath, times_constant) +{ + double omega = 2.; + int fcn0 = func13_newBasic("sin", omega); + double A = 1.234; + int fcn = func13_newModified("times-constant", fcn0, A); + ASSERT_GE(fcn, 0); + EXPECT_DOUBLE_EQ(func13_eval(fcn, 0.5), sin(omega * 0.5) * A); +} + +TEST(ctmath, plus_constant) +{ + double omega = 2.; + int fcn0 = func13_newBasic("sin", omega); + double A = 1.234; + int fcn = func13_newModified("plus-constant", fcn0, A); + ASSERT_GE(fcn, 0); + EXPECT_DOUBLE_EQ(func13_eval(fcn, 0.5), sin(omega * 0.5) + A); +} + +TEST(ctmath, periodic) +{ + double omega = 2.; + int fcn0 = func13_newBasic("sin", omega); + double A = 1.234; + int fcn = func13_newModified("periodic", fcn0, A); + ASSERT_GE(fcn, 0); + EXPECT_DOUBLE_EQ(func13_eval(fcn, 0.5), func13_eval(fcn, 0.5 + A)); +} + +TEST(ctmath, generic) +{ + // Composite function reported in issue #752 + + vector params = {0., 0., 1., 1.}; + int fs = func13_newAdvanced("Fourier", params.size(), params.data()); // sin(x) + auto func13_s = newFunc1("sin", 1.); + EXPECT_DOUBLE_EQ(func13_eval(fs, 0.5), func13_s->eval(0.5)); + + int fs2 = func13_newCompound("product", fs, fs); // (sin(x)^2) + auto func13_s2 = newFunc1("product", func13_s, func13_s); + EXPECT_DOUBLE_EQ(func13_eval(fs2, 0.5), func13_s2->eval(0.5)); + + double one = 1.; + int fs1 = func13_newAdvanced("polynomial3", 1, &one); // 1. + auto func13_s1 = newFunc1("constant", 1.); + EXPECT_DOUBLE_EQ(func13_eval(fs1, 0.5), func13_s1->eval(0.5)); + EXPECT_DOUBLE_EQ(func13_eval(fs1, 0.5), 1.); + + int fs2_1 = func13_newCompound("diff", fs1, fs2); // 1-(sin(x))^2 + auto func13_s2_1 = newFunc1("diff", func13_s1, func13_s2); + EXPECT_DOUBLE_EQ(func13_eval(fs2_1, 0.5), func13_s2_1->eval(0.5)); + + vector p_arr = {1., .5, 0.}; + int fs3 = func13_newAdvanced("Arrhenius", 3, p_arr.data()); // sqrt function + auto func13_s3 = newFunc1("Arrhenius", p_arr); + EXPECT_DOUBLE_EQ(func13_eval(fs3, 0.5), func13_s3->eval(0.5)); + + // overall composite function + int fs4 = func13_newCompound("composite", fs3, fs2_1); // sqrt(1-(sin(x))^2) + auto func13_s4 = newFunc1("composite", func13_s3, func13_s2_1); + EXPECT_DOUBLE_EQ(func13_eval(fs4, 0.5), func13_s4->eval(0.5)); + + // an easier equivalent expression (using trigonometry) + auto func13_s5 = newFunc1("cos", 1.); // missing the absolute value + EXPECT_DOUBLE_EQ(func13_eval(fs4, 0.5), func13_s5->eval(0.5)); + EXPECT_DOUBLE_EQ(func13_eval(fs4, 3.5), -func13_s5->eval(3.5)); +}