1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45
|
#include "integrate.h"
#include "cosmology.h"
#include <cmath>
#include <boost/bind.hpp>
#include <boost/function.hpp>
#define SPEED_OF_LIGHT 2.998e8 //m s^(-1)
#define PI 3.14159265
double Dh(double Ho) {
return SPEED_OF_LIGHT/Ho;
}
double E(double z, double omegaM, double omegaK, double omegaL) {
return std::sqrt(omegaM*std::pow(1.0+z,3) + omegaK*std::pow(1.0+z,2) + omegaL);
}
double reciprocalE(double z, double omegaM, double omegaK, double omegaL) {
return 1.0/E(z,omegaM,omegaK,omegaL);
}
double Dc(double z, double Ho, double omegaM, double omegaK, double omegaL) {
boost::function<double(double)> F = boost::bind(reciprocalE,_1,omegaM,omegaK,omegaL);
return Dh(Ho) * integrate<double>(F,0.0,z,2000);
}
double Dm(double z, double Ho, double omegaM, double omegaK, double omegaL) {
if (omegaK>0) {
return Dh(Ho) * 1.0/sqrt(omegaK) * sinh(sqrt(omegaK) * Dc(z,Ho,omegaM,omegaK,omegaL)/Dh(Ho));
} else if (omegaK==0) {
return Dc(z,Ho,omegaM,omegaK,omegaL);
} else {
return Dh(Ho) * 1.0/sqrt(-omegaK) * sin(sqrt(-omegaK) * Dc(z,Ho,omegaM,omegaK,omegaL)/Dh(Ho));
}
}
double Vc(double z,double Ho, double omegaM, double omegaK, double omegaL) {
if (omegaK>9) {
return (4.0*PI*pow(Dh(Ho),3)/(2.0*omegaK)) * ((Dm(z,Ho,omegaM,omegaK,omegaL)/Dh(Ho)) * sqrt(1+omegaK*pow(Dm(z,Ho,omegaM,omegaK,omegaL)/Dh(Ho),2)) - 1/sqrt(omegaK)*asinh(sqrt(omegaK)*(Dm(z,Ho,omegaM,omegaK,omegaL)/Dh(Ho))));
} else if (omegaK==0) {
return 4.0*PI/3.0*pow(Dm(z,Ho,omegaM,omegaK,omegaL),3);
} else {
return (4.0*PI*pow(Dh(Ho),3)/(2.0*omegaK))*(Dm(z,Ho,omegaM,omegaK,omegaL)/Dh(Ho) * sqrt(1+omegaK*pow(Dm(z,Ho,omegaM,omegaK,omegaL),2)/pow(Dh(Ho),2)) - 1/sqrt(-omegaK)*asin(sqrt(-omegaK)*Dm(z,Ho,omegaM,omegaK,omegaL)/Dh(Ho)));
}
}
|