C
cupiemayo
Hello. I'm trying to make a sort of generic integral class, holding
the boundary values and the integrand. Eventually I'd like to use it
inside a class, like the example below. I think the problem is that
in order to pass a member function, I have to qualify the name a bit
differently, for example:
typedef double (astro::universe::*pfn)(double);
(by the way, I don't think this actually worked...) But doing this
sort of thing of course doesn't make the integral class generic enough.
I'd like to get some advices as to how I can make such an integral
class that even a member function can be passed easily as the
integrand. How can you do this using C++?
Thanks in advance...
----------------------------------------------------------------
#include <cmath>
#include <iostream>
#include <string>
#include <vector>
namespace math {
typedef double (*pfn)(double);
class integral
{
public:
integral(double a, double b, pfn f)
: lower_(a), upper_(b), integrand_(f)
{}
double lower_bound() const { return lower_; }
double upper_bound() const { return upper_; }
void change_bounds(double a, double b)
{
lower_ = a;
upper_ = b;
}
double simpson(unsigned int const n) const
{
double h = (upper_ - lower_) / n;
double sum = integrand_(lower_) * 0.5;
for (int i = 1; i < n; ++i)
sum += integrand_(lower_ + i * h);
sum += integrand_(upper_) * 0.5;
double summid = 0.0;
for (int i = 1; i <= n; ++i)
summid += integrand_(lower_ + (i - 0.5)*h);
return (sum + 2 * summid) * h / 3.0;
}
private:
double lower_;
double upper_;
pfn integrand_;
};
} // namespace math
namespace astro {
class universe
{
public:
universe(double const omega_matter,
double const omega_lambda,
double const hubble_const)
: omega_matter_(omega_matter),
omega_lambda_(omega_lambda),
hubble_const_(hubble_const)
{}
double comoving_distance(double z)
{
//
// Want to use a generic integral class above with
// the integrand astro::universe::inverse_e(double),
// but does not compile....
//
math::integral func(0.0, z, inverse_e);
double speed_of_light = 3.0e8;
double c_cgs = speed_of_light * (100.0);
double H0_cgs = hubble_const_ * (1.0 / 3.0856776e19);
return (c_cgs / H0_cgs) * func.simpson(1000);
}
private:
double hubble_const_;
double omega_matter_;
double omega_lambda_;
double inverse_e(double z)
{
return 1.0 / std::sqrt(omega_matter_ * std:ow(1 + z, 3)
+ omega_lambda_);
}
};
} // namespace astro
int main()
{
//
// Create a universe with a particular cosmology.
//
astro::universe uni(0.3, 0.7, 70.0);
//
// Compute a comoving distance to an object at redshift 0.4.
//
std::cout << uni.comoving_distance(0.4) << '\n';
return 0;
}
the boundary values and the integrand. Eventually I'd like to use it
inside a class, like the example below. I think the problem is that
in order to pass a member function, I have to qualify the name a bit
differently, for example:
typedef double (astro::universe::*pfn)(double);
(by the way, I don't think this actually worked...) But doing this
sort of thing of course doesn't make the integral class generic enough.
I'd like to get some advices as to how I can make such an integral
class that even a member function can be passed easily as the
integrand. How can you do this using C++?
Thanks in advance...
----------------------------------------------------------------
#include <cmath>
#include <iostream>
#include <string>
#include <vector>
namespace math {
typedef double (*pfn)(double);
class integral
{
public:
integral(double a, double b, pfn f)
: lower_(a), upper_(b), integrand_(f)
{}
double lower_bound() const { return lower_; }
double upper_bound() const { return upper_; }
void change_bounds(double a, double b)
{
lower_ = a;
upper_ = b;
}
double simpson(unsigned int const n) const
{
double h = (upper_ - lower_) / n;
double sum = integrand_(lower_) * 0.5;
for (int i = 1; i < n; ++i)
sum += integrand_(lower_ + i * h);
sum += integrand_(upper_) * 0.5;
double summid = 0.0;
for (int i = 1; i <= n; ++i)
summid += integrand_(lower_ + (i - 0.5)*h);
return (sum + 2 * summid) * h / 3.0;
}
private:
double lower_;
double upper_;
pfn integrand_;
};
} // namespace math
namespace astro {
class universe
{
public:
universe(double const omega_matter,
double const omega_lambda,
double const hubble_const)
: omega_matter_(omega_matter),
omega_lambda_(omega_lambda),
hubble_const_(hubble_const)
{}
double comoving_distance(double z)
{
//
// Want to use a generic integral class above with
// the integrand astro::universe::inverse_e(double),
// but does not compile....
//
math::integral func(0.0, z, inverse_e);
double speed_of_light = 3.0e8;
double c_cgs = speed_of_light * (100.0);
double H0_cgs = hubble_const_ * (1.0 / 3.0856776e19);
return (c_cgs / H0_cgs) * func.simpson(1000);
}
private:
double hubble_const_;
double omega_matter_;
double omega_lambda_;
double inverse_e(double z)
{
return 1.0 / std::sqrt(omega_matter_ * std:ow(1 + z, 3)
+ omega_lambda_);
}
};
} // namespace astro
int main()
{
//
// Create a universe with a particular cosmology.
//
astro::universe uni(0.3, 0.7, 70.0);
//
// Compute a comoving distance to an object at redshift 0.4.
//
std::cout << uni.comoving_distance(0.4) << '\n';
return 0;
}