Loading [MathJax]/extensions/TeX/AMSsymbols.js
LALBurst 2.0.7.1-b246709
All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Macros Modules Pages
cs_gamma.c
Go to the documentation of this file.
1/*
2* Copyright (C) 2007 Xavier Siemens
3* Copyright (C) 2010 Andrew Mergl
4* Copyright (C) 2018 Daichi Tsuna
5*
6* This program is free software; you can redistribute it and/or modify
7* it under the terms of the GNU General Public License as published by
8* the Free Software Foundation; either version 2 of the License, or
9* (at your option) any later version.
10*
11* This program is distributed in the hope that it will be useful,
12* but WITHOUT ANY WARRANTY; without even the implied warranty of
13* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14* GNU General Public License for more details.
15*
16* You should have received a copy of the GNU General Public License
17* along with with program; see the file COPYING. If not, write to the
18* Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
19* MA 02110-1301 USA
20*/
21
22/*********************************************************************************/
23/* Cosmic string burst rate computation code for small/large loops */
24/* */
25/* Xavier Siemens, Jolien Creighton, Irit Maor */
26/* */
27/* UWM/Caltech - September 2006 */
28/*********************************************************************************/
29/*Modified June 2010 by Andrew Mergl for use with Python*/
30/* Modified November 2018 by Daichi Tsuna to include large loop scenarios */
31/* Updated to be able to test 3 models of large loop distributions */
32
33/* Ignore warnings in Python API itself */
34#pragma GCC diagnostic push
35#pragma GCC diagnostic ignored "-Wcast-align"
36#include <Python.h>
37#pragma GCC diagnostic pop
38
39#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION
40#include <numpy/arrayobject.h>
41
42#include <math.h>
43#include <stdlib.h>
44#include <string.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>
50
51#define CUSPS_PER_LOOP 1.0 /* c */
52#define LOOP_RAD_POWER 50.0 /* Gamma */
53
54#define H0 LAMBDA_H_0
55
56
57/*****************************************************************************/
58/*int finddRdz(double Gmu, double alpha, double f, double Gamma, int Namp, double *zofA, double *dRdz)
59 * Find dR/dz given Gmu, alpha, f, and Gamma. This is a C function that was taken from cs_gamma.c
60 * and modified so that it could be called from Python. The global variables
61 * that the function used to use are now all passed through arguments.
62 *
63 * Arguments:
64 * Gmu, alpha, gamma: parameters calculated by the main program. See technical
65 * document for what they represent cosmologically.
66 * f: the frequency that is passed to the main program with --frequency opt
67 * Namp: The size of the data arrays. This is set when opening the data file
68 * *zofA, *dRdz: 1D arrays of length Namp. See technical document for what
69 * they represent cosmologically.
70 */
71/*****************************************************************************/
72static PyObject *cs_gamma_finddRdz(PyObject *self, PyObject *args)
73{
74 PyArrayObject *Numpy_zofA;
75 PyObject *Numpy_dRdz;
76 double Gmu, alpha, f, Gamma, *zofA, *dRdz;
77 long int Namp;
78 cs_cosmo_functions_t cosmofns;
79 int j;
80 (void)self; /* silence unused parameter warning */
81
82 if (!PyArg_ParseTuple(args, "ddddO!", &Gmu, &alpha, &f, &Gamma, &PyArray_Type, &Numpy_zofA))
83 return NULL;
84
85 Numpy_zofA = PyArray_GETCONTIGUOUS(Numpy_zofA);
86 if(!Numpy_zofA)
87 return NULL;
88 Namp = PyArray_DIM(Numpy_zofA, 0);
89 zofA = PyArray_DATA(Numpy_zofA);
90
91 {
92 npy_intp dims[1] = {Namp};
93 Numpy_dRdz = PyArray_SimpleNew(1, dims, NPY_DOUBLE);
94 }
95 dRdz = PyArray_DATA((PyArrayObject *) Numpy_dRdz);
96
97 cosmofns = XLALCSCosmoFunctions( zofA, Namp);
98
99 for ( j = 0; j < Namp; j++ )
100 {
101 /*double theta = pow((1+cosmofns.z[j]) * f * alpha * cosmofns.phit[j] / H0, -1.0/3.0);
102
103 if (theta > 1.0)
104 dRdz[j] = 0.0;
105 else*/
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);
109 Numpy_dRdz = NULL;
110 break;
111 }
112 }
113
114 XLALCSCosmoFunctionsFree( cosmofns );
115 Py_DECREF(Numpy_zofA);
116
117 return Numpy_dRdz;
118}
119/*****************************************************************************/
120/*int findzofA(double Gmu, double alpha, int Namp, double *zofA, double *amp)
121 *Find z(A) given Gmu and alpha. This function was taken from cs_gamma.c and
122 * modified so that it could be called from Python. The global variables it
123 * used to use are now passed to it.
124 *
125 * Arugments:
126 * Gmu, alpha: values calculated by the main program. See S4 technical
127 * documentation for their cosmological meanings.
128 * Namp: The length of the data arrays. This is set when reading in the data
129 * zofA, amp: 1D data arrays of length Namp. See S4 technical documentation
130 * for their cosmological meanings.
131 */
132/*****************************************************************************/
133static PyObject *cs_gamma_findzofA(PyObject *self, PyObject *args)
134{
135 PyArrayObject *Numpy_amp;
136 PyObject *Numpy_zofA;
137 double Gmu, alpha, *zofA, *amp;
138 unsigned long int Namp;
139 (void)self; /* silence unused parameter warning */
140
141 double z_min = 1e-20, z_max = 1e10;
142 double dlnz = 0.05;
143 unsigned numz = floor( (log(z_max) - log(z_min)) / dlnz );
144 unsigned long int i;
145 cs_cosmo_functions_t cosmofns;
146 double *fz,*z;
147 double a;
148 gsl_interp *zofa_interp;
149 gsl_interp_accel *acc_zofa = gsl_interp_accel_alloc();
150
151 if (!PyArg_ParseTuple(args, "ddO!", &Gmu, &alpha, &PyArray_Type, &Numpy_amp))
152 return NULL;
153
154 Numpy_amp = PyArray_GETCONTIGUOUS(Numpy_amp);
155 if(!Numpy_amp)
156 return NULL;
157 Namp = PyArray_DIM(Numpy_amp, 0);
158 amp = PyArray_DATA(Numpy_amp);
159
160 {
161 npy_intp dims[1] = {Namp};
162 Numpy_zofA = PyArray_SimpleNew(1, dims, NPY_DOUBLE);
163 }
164 zofA = PyArray_DATA((PyArrayObject *) Numpy_zofA);
165
166 cosmofns = XLALCSCosmoFunctionsAlloc( z_min, dlnz, numz );
167
168 zofa_interp = gsl_interp_alloc (gsl_interp_linear, cosmofns.n);
169
170 fz = calloc( cosmofns.n, sizeof( *fz ) );
171 z = calloc( cosmofns.n, sizeof( *z ) );
172
173 /* first compute the function that relates A and z */
174 /* invert order; b/c fz is a monotonically decreasing func of z */
175 for ( i = cosmofns.n ; i > 0; i-- )
176 {
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];
180 }
181
182 gsl_interp_init (zofa_interp, fz, z, cosmofns.n);
183
184 /* now compute the amplitudes (suitably multiplied) that are equal to fz for some z*/
185 for ( i = 0; i < Namp; i++ )
186 {
187 a = amp[i] * pow(H0,-1.0/3.0) * pow(alpha,-2.0/3.0) / Gmu;
188 /* evaluate z(fz) at fz=a */
189 zofA[i] = gsl_interp_eval (zofa_interp, fz, z, a, acc_zofa );
190 if(gsl_isnan(zofA[i])) {
191 Py_DECREF(Numpy_zofA);
192 Numpy_zofA = NULL;
193 break;
194 }
195 }
196
197 XLALCSCosmoFunctionsFree( cosmofns );
198 Py_DECREF(Numpy_amp);
199 free(fz);
200 free(z);
201 gsl_interp_free (zofa_interp);
202 gsl_interp_accel_free(acc_zofa);
203
204 return Numpy_zofA;
205}
206
207
208/*******************************************************************************/
209/* nu(double Gmu, double Gamma, double *amp, double *z, char model)
210 * Find loop distribution givem Gmu, Gamma, amplitude, z, and a particular large loop model.
211 * The loop distribution nu is calculated as a function of amp and z, and is
212 * returned as an array having a size of Namp * Nz.
213 */
214/*****************************************************************************/
215static double nu(double Gmu, double Gamma, double A, double z, double phit, double phiA, const char *model)
216{
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;
224 double t=phit/H0;
225 double cuspdist;
226 double teq=8.122570474611143e+11;
227 double crateR, crateRadStragglers, crateM;
228 double crate;
229
230 if (strcmp(model,"Siemens06") == 0) {
231 alpha=1e-1;
232 nuR=0.4*15*sqrt(alpha);
233 nuM=0.12*4;
234 /* Radiation era loops */
235 crateR = 0.0;
236 if( (l < alpha*t) && (t < teq) )
237 crateR = nuR * pow(t,-1.5) * pow( l + Gamma*Gmu/H0 * phit, -2.5 );
238
239 /* Radiation stragglers */
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);
243
244 /* matter era loops */
245 crateM = 0.0;
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);
248
249 crate = crateR + crateRadStragglers + crateM;
250 } else if (strcmp(model,"Blanco-Pillado14") == 0) {
251 alpha = 0.1;
252 nuR = 0.18;
253 /* Radiation era loops */
254 crateR = 0.0;
255 if( (l < alpha*t) && (t < teq) )
256 crateR = nuR * pow(t,-1.5) * pow( l + Gamma*Gmu/H0 * phit, -2.5 );
257
258 /* Radiation stragglers */
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);
262
263 /* matter era loops */
264 crateM = 0.0;
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);
267
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);
272 crateR = 0.0;
273 crateM = 0.0;
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);
276 /* Radiation era loops */
277 if (t <= teq){
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);
280
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);
283
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);
286 }
287 /* Matter era loops */
288 if (t > teq){
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);
291
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);
294
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);
297 }
298
299 crate = crateR + crateM;
300 }
301
302 cuspdist = 3.0/A * crate;
303
304 return cuspdist;
305}
306
307/*******************************************************************************/
308/* finddRdzdA(double Gmu, double f, double Gamma, double *amp, double *z, double *nu)
309 * Find d^2R/dzdA givem Gmu, f, Gamma, amplitude, z, and the loop distribution.
310 * Using the equations formulated in e.g. arXiv:1712.01168 Appendix B, d^2R/dzdA is calculated
311 * as a function of amp and z, and is returned as an array having a size of Namp * Nz.
312 */
313/*****************************************************************************/
314static PyObject *cs_gamma_finddRdzdA(PyObject *self, PyObject *args)
315{
316 PyArrayObject *Numpy_amp, *Numpy_z;
317 PyObject *Numpy_dRdzdA;
318 double Gmu, f, Gamma, *amp, *z;
319 long int Namp, Nz;
320 int i,j;
321 char *model;
322 cs_cosmo_functions_t cosmofns;
323 (void)self; /* silence unused parameter warning */
324
325 if (!PyArg_ParseTuple(args, "dddO!O!s", &Gmu, &f, &Gamma, &PyArray_Type, &Numpy_amp, &PyArray_Type, &Numpy_z, &model))
326 return NULL;
327
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);
332 Py_DECREF(Numpy_z);
333 return NULL;
334 }
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);
339
340 {
341 npy_intp dims[] = {Namp,Nz};
342 Numpy_dRdzdA = PyArray_SimpleNew(2, dims, NPY_DOUBLE);
343 }
344
345 cosmofns = XLALCSCosmoFunctions(z, Nz);
346
347 /* loop over amplitudes */
348 for ( i = 0; i < Namp; i++ ){
349 double A = amp[i];
350 /* loop over redshifts */
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);
357 double dRdzdA;
358
359 if(theta > 1.0){
360 dRdzdA = 0.0;
361 } else {
362 double Delta = 0.25*theta*theta;
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);
366 Numpy_dRdzdA = NULL;
367 XLALCSCosmoFunctionsFree( cosmofns );
368 Py_DECREF(Numpy_amp);
369 Py_DECREF(Numpy_z);
370 return NULL;
371 }
372 }
373 *(double *) PyArray_GETPTR2((PyArrayObject *) Numpy_dRdzdA, i, j) = dRdzdA;
374 }
375 }
376
377 XLALCSCosmoFunctionsFree( cosmofns );
378 Py_DECREF(Numpy_amp);
379 Py_DECREF(Numpy_z);
380
381 return Numpy_dRdzdA;
382}
383
384
385/*******************************************************************************/
386
387//List of functions available to the Python module.
388static PyMethodDef cs_gammaMethods[] = {
389 {"findzofA", cs_gamma_findzofA, METH_VARARGS,
390 "Function to find z(A). From cs_gamma.c; modified to be called from Python."},
391 {"finddRdz", cs_gamma_finddRdz, METH_VARARGS,
392 "Function to find dR/dz. From cs_gamma.c; modified to be called from Python."},
393 {"finddRdzdA", cs_gamma_finddRdzdA, METH_VARARGS,
394 "Function to find d^2R/dzdA. From cs_gammaLargeLoops.c; modified to be called from Python."},
395 {NULL, NULL, 0, NULL}
396};
397
398//Then Python module initialization function.
399PyMODINIT_FUNC PyInit_cs_gamma(void); /* silence -Wmissing-prototypes */
400PyMODINIT_FUNC PyInit_cs_gamma(void)
401{
402 static PyModuleDef moduledef = {
403 PyModuleDef_HEAD_INIT,
404 "cs_gamma", NULL, -1, cs_gammaMethods,
405 NULL, NULL, NULL, NULL
406 };
407 PyObject *module = PyModule_Create(&moduledef);
408
409 import_array();
410
411 return module;
412}
int l
double i
double e
double theta
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)
#define H0
Definition: cs_gamma.c:54
static double nu(double Gmu, double Gamma, double A, double z, double phit, double phiA, const char *model)
Definition: cs_gamma.c:215
static PyObject * cs_gamma_finddRdzdA(PyObject *self, PyObject *args)
Definition: cs_gamma.c:314
static PyObject * cs_gamma_finddRdz(PyObject *self, PyObject *args)
Definition: cs_gamma.c:72
static PyMethodDef cs_gammaMethods[]
Definition: cs_gamma.c:388
static PyObject * cs_gamma_findzofA(PyObject *self, PyObject *args)
Definition: cs_gamma.c:133
PyMODINIT_FUNC PyInit_cs_gamma(void)
Definition: cs_gamma.c:400
static const INT4 a
double alpha