34#pragma GCC diagnostic push
35#pragma GCC diagnostic ignored "-Wcast-align"
37#pragma GCC diagnostic pop
39#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION
40#include <numpy/arrayobject.h>
45#include <gsl/gsl_math.h>
46#include <gsl/gsl_interp.h>
47#include <gsl/gsl_errno.h>
48#include <lal/cs_cosmo.h>
49#include <lal/cs_lambda_cosmo.h>
51#define CUSPS_PER_LOOP 1.0
52#define LOOP_RAD_POWER 50.0
74 PyArrayObject *Numpy_zofA;
76 double Gmu,
alpha, f, Gamma, *zofA, *dRdz;
82 if (!PyArg_ParseTuple(
args,
"ddddO!", &Gmu, &
alpha, &f, &Gamma, &PyArray_Type, &Numpy_zofA))
85 Numpy_zofA = PyArray_GETCONTIGUOUS(Numpy_zofA);
88 Namp = PyArray_DIM(Numpy_zofA, 0);
89 zofA = PyArray_DATA(Numpy_zofA);
92 npy_intp dims[1] = {Namp};
93 Numpy_dRdz = PyArray_SimpleNew(1, dims, NPY_DOUBLE);
95 dRdz = PyArray_DATA((PyArrayObject *) Numpy_dRdz);
99 for ( j = 0; j < Namp; j++ )
106 dRdz[j] = 0.5 *
H0 * pow(f/
H0,-2.0/3.0) * pow(
alpha, -5.0/3.0) / (Gamma*Gmu) * pow(cosmofns.
phit[j],-14.0/3.0) * cosmofns.
phiV[j] * pow(1+cosmofns.
z[j],-5.0/3.0);
107 if(gsl_isnan(dRdz[j])) {
108 Py_DECREF(Numpy_dRdz);
115 Py_DECREF(Numpy_zofA);
135 PyArrayObject *Numpy_amp;
136 PyObject *Numpy_zofA;
137 double Gmu,
alpha, *zofA, *amp;
138 unsigned long int Namp;
141 double z_min = 1
e-20, z_max = 1e10;
143 unsigned numz = floor( (log(z_max) - log(z_min)) / dlnz );
148 gsl_interp *zofa_interp;
149 gsl_interp_accel *acc_zofa = gsl_interp_accel_alloc();
151 if (!PyArg_ParseTuple(
args,
"ddO!", &Gmu, &
alpha, &PyArray_Type, &Numpy_amp))
154 Numpy_amp = PyArray_GETCONTIGUOUS(Numpy_amp);
157 Namp = PyArray_DIM(Numpy_amp, 0);
158 amp = PyArray_DATA(Numpy_amp);
161 npy_intp dims[1] = {Namp};
162 Numpy_zofA = PyArray_SimpleNew(1, dims, NPY_DOUBLE);
164 zofA = PyArray_DATA((PyArrayObject *) Numpy_zofA);
168 zofa_interp = gsl_interp_alloc (gsl_interp_linear, cosmofns.
n);
170 fz = calloc( cosmofns.
n,
sizeof( *fz ) );
171 z = calloc( cosmofns.
n,
sizeof( *z ) );
175 for (
i = cosmofns.
n ;
i > 0;
i-- )
177 unsigned long int j = cosmofns.
n -
i;
178 z[j] = cosmofns.
z[
i-1];
179 fz[j] = pow(cosmofns.
phit[
i-1], 2.0/3.0) * pow(1+z[j], -1.0/3.0) / cosmofns.
phiA[
i-1];
182 gsl_interp_init (zofa_interp, fz, z, cosmofns.
n);
185 for (
i = 0;
i < Namp;
i++ )
187 a = amp[
i] * pow(
H0,-1.0/3.0) * pow(
alpha,-2.0/3.0) / Gmu;
189 zofA[
i] = gsl_interp_eval (zofa_interp, fz, z,
a, acc_zofa );
190 if(gsl_isnan(zofA[
i])) {
191 Py_DECREF(Numpy_zofA);
198 Py_DECREF(Numpy_amp);
201 gsl_interp_free (zofa_interp);
202 gsl_interp_accel_free(acc_zofa);
215static double nu(
double Gmu,
double Gamma,
double A,
double z,
double phit,
double phiA,
const char *model)
217 double l = pow( A / Gmu /
H0 * pow(1.0+z,1.0/3.0)* phiA, 3.0/2.0);
218 double alpha, nuR, nuM;
219 double P_R=1.60, P_M=1.41, a_index_R=1.0/2.0, a_index_M=2.0/3.0;
220 double Upsilon = 10.0;
221 double chi_R=1.0-P_R/2.0;
222 double chi_M=1.0-P_M/2.0;
223 double gamma_c_R, gamma_c_M;
226 double teq=8.122570474611143e+11;
227 double crateR, crateRadStragglers, crateM;
230 if (strcmp(model,
"Siemens06") == 0) {
232 nuR=0.4*15*sqrt(
alpha);
236 if( (
l <
alpha*t) && (t < teq) )
237 crateR = nuR * pow(t,-1.5) * pow(
l + Gamma*Gmu/
H0 * phit, -2.5 );
240 crateRadStragglers = 0.0;
241 if ( (
l <
alpha*teq-Gamma*Gmu*(t-teq) ) && ( t > teq ) )
242 crateRadStragglers = nuR*pow(teq, 0.5)/t/t* pow(
l + Gamma*Gmu/
H0 * phit, -2.5);
246 if( (
l <
alpha*t) && (t > teq) && (
l >
alpha*teq-Gamma*Gmu*(t-teq)) )
247 crateM = nuM / t / t / (
l + Gamma*Gmu/
H0 * phit) / (
l + Gamma*Gmu/
H0 * phit);
249 crate = crateR + crateRadStragglers + crateM;
250 }
else if (strcmp(model,
"Blanco-Pillado14") == 0) {
255 if( (
l <
alpha*t) && (t < teq) )
256 crateR = nuR * pow(t,-1.5) * pow(
l + Gamma*Gmu/
H0 * phit, -2.5 );
259 crateRadStragglers = 0.0;
260 if ( (
l <
alpha*teq-Gamma*Gmu*(t-teq) ) && ( t > teq ) )
261 crateRadStragglers = nuR*pow(teq, 0.5)/t/t* pow(
l + Gamma*Gmu/
H0 * phit, -2.5);
265 if( (
l < 0.18*t) && (t > teq) && (
l >
alpha*teq-Gamma*Gmu*(t-teq)) )
266 crateM = (0.27-0.45*pow(
l/t,0.31)) / t / t / (
l + Gamma*Gmu/
H0 * phit) / (
l + Gamma*Gmu/
H0 * phit);
268 crate = crateR + crateRadStragglers + crateM;
269 }
else if (strcmp(model,
"Ringeval07") == 0) {
270 nuR = 0.21 * pow(1.0-a_index_R, 3.0-P_R);
271 nuM = 0.09 * pow(1.0-a_index_M, 3.0-P_M);
274 gamma_c_R = Upsilon * pow(Gmu,1.0+2.0*chi_R);
275 gamma_c_M = Upsilon * pow(Gmu,1.0+2.0*chi_M);
278 if( (
l > Gamma * Gmu * t) && (
l <= t/(1.0-a_index_R)) )
279 crateR = pow(t,-4.0) * nuR / pow(
l/t + Gamma*Gmu, P_R+1);
281 if( (
l > gamma_c_R * t) && (
l <= Gamma * Gmu * t) )
282 crateR = pow(t,-4.0) * nuR * (3.0*a_index_R-2.0*chi_R-1.0) / (2.0-2.0*chi_R) / (Gamma*Gmu) / pow(
l/t,P_R);
284 if(
l <= gamma_c_R * t)
285 crateR = pow(t,-4.0) * nuR * (3.0*a_index_R-2.0*chi_R-1.0) / (2.0-2.0*chi_R) / (Gamma*Gmu) / pow(gamma_c_R,P_R);
289 if( (
l > Gamma * Gmu * t) && (
l <= t/(1.0-a_index_M)) )
290 crateM = pow(t,-4.0) * nuM / pow(
l/t + Gamma*Gmu, P_M+1);
292 if( (
l > gamma_c_M * t) && (
l <= Gamma * Gmu * t) )
293 crateM = pow(t,-4.0) * nuM * (3.0*a_index_M-2.0*chi_M-1.0) / (2.0-2.0*chi_M) / (Gamma*Gmu) / pow(
l/t,P_M);
295 if(
l <= gamma_c_M * t)
296 crateM = pow(t,-4.0) * nuM * (3.0*a_index_M-2.0*chi_M-1.0) / (2.0-2.0*chi_M) / (Gamma*Gmu) / pow(gamma_c_M,P_M);
299 crate = crateR + crateM;
302 cuspdist = 3.0/A * crate;
316 PyArrayObject *Numpy_amp, *Numpy_z;
317 PyObject *Numpy_dRdzdA;
318 double Gmu, f, Gamma, *amp, *z;
325 if (!PyArg_ParseTuple(
args,
"dddO!O!s", &Gmu, &f, &Gamma, &PyArray_Type, &Numpy_amp, &PyArray_Type, &Numpy_z, &model))
328 Numpy_amp = PyArray_GETCONTIGUOUS(Numpy_amp);
329 Numpy_z = PyArray_GETCONTIGUOUS(Numpy_z);
330 if(!Numpy_amp || !Numpy_z){
331 Py_DECREF(Numpy_amp);
335 Namp = PyArray_DIM(Numpy_amp, 0);
336 amp = PyArray_DATA(Numpy_amp);
337 Nz = PyArray_DIM(Numpy_z, 0);
338 z = PyArray_DATA(Numpy_z);
341 npy_intp dims[] = {Namp,Nz};
342 Numpy_dRdzdA = PyArray_SimpleNew(2, dims, NPY_DOUBLE);
348 for (
i = 0;
i < Namp;
i++ ){
351 for ( j = 0; j < Nz; j++ ){
352 double phit = cosmofns.
phit[j];
353 double phiA = cosmofns.
phiA[j];
354 double phiV = cosmofns.
phiV[j];
355 double l = pow ( A / Gmu /
H0 * pow(1+z[j],1.0/3.0)* phiA, 3.0/2.0);
356 double theta = pow(f*(1+z[j])*
l, -1.0/3.0);
363 dRdzdA = pow(
H0,-3.0) * phiV / (1.0+z[j]) *
nu(Gmu, Gamma, A, z[j], phit, phiA, model) * Delta;
364 if(gsl_isnan(dRdzdA)) {
365 Py_DECREF(Numpy_dRdzdA);
368 Py_DECREF(Numpy_amp);
373 *(
double *) PyArray_GETPTR2((PyArrayObject *) Numpy_dRdzdA,
i, j) = dRdzdA;
378 Py_DECREF(Numpy_amp);
390 "Function to find z(A). From cs_gamma.c; modified to be called from Python."},
392 "Function to find dR/dz. From cs_gamma.c; modified to be called from Python."},
394 "Function to find d^2R/dzdA. From cs_gammaLargeLoops.c; modified to be called from Python."},
395 {NULL, NULL, 0, NULL}
402 static PyModuleDef moduledef = {
403 PyModuleDef_HEAD_INIT,
405 NULL, NULL, NULL, NULL
407 PyObject *module = PyModule_Create(&moduledef);
cs_cosmo_functions_t XLALCSCosmoFunctionsAlloc(double zmin, double dlnz, size_t n)
void XLALCSCosmoFunctionsFree(cs_cosmo_functions_t)
cs_cosmo_functions_t XLALCSCosmoFunctions(double *z, size_t n)
static double nu(double Gmu, double Gamma, double A, double z, double phit, double phiA, const char *model)
static PyObject * cs_gamma_finddRdzdA(PyObject *self, PyObject *args)
static PyObject * cs_gamma_finddRdz(PyObject *self, PyObject *args)
static PyMethodDef cs_gammaMethods[]
static PyObject * cs_gamma_findzofA(PyObject *self, PyObject *args)
PyMODINIT_FUNC PyInit_cs_gamma(void)