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::univers e::*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(d ouble a, double b)

{

lower_ = a;

upper_ = b;

}

double simpson(unsigne d int const n) const

{

double h = (upper_ - lower_) / n;

double sum = integrand_(lowe r_) * 0.5;

for (int i = 1; i < n; ++i)

sum += integrand_(lowe r_ + i * h);

sum += integrand_(uppe r_) * 0.5;

double summid = 0.0;

for (int i = 1; i <= n; ++i)

summid += integrand_(lowe r_ + (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_(o mega_matter),

omega_lambda_(o mega_lambda),

hubble_const_(h ubble_const)

{}

double comoving_distan ce(double z)

{

//

// Want to use a generic integral class above with

// the integrand astro::universe ::inverse_e(dou ble),

// 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(10 00);

}

private:

double hubble_const_;

double omega_matter_;

double omega_lambda_;

double inverse_e(doubl e z)

{

return 1.0 / std::sqrt(omega _matter_ * std::pow(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_di stance(0.4) << '\n';

return 0;

}