LALSimulation  5.4.0.1-fe68b98
LALSimIMRSEOBNRv2ChirpTime.c
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2014 Michael Puerrer
3  *
4  * This program is free software; you can redistribute it and/or modify
5  * it under the terms of the GNU General Public License as published by
6  * the Free Software Foundation; either version 2 of the License, or
7  * (at your option) any later version.
8  *
9  * This program is distributed in the hope that it will be useful,
10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12  * GNU General Public License for more details.
13  *
14  * You should have received a copy of the GNU General Public License
15  * along with with program; see the file COPYING. If not, write to the
16  * Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
17  * MA 02110-1301 USA
18  */
19 
20 #ifdef __GNUC__
21 #define UNUSED __attribute__ ((unused))
22 #else
23 #define UNUSED
24 #endif
25 
26 #include <stdio.h>
27 #include <stdlib.h>
28 #include <math.h>
29 #include <complex.h>
30 #include <time.h>
31 #include <unistd.h>
32 #include <getopt.h>
33 #include <stdbool.h>
34 #include <alloca.h>
35 #include <string.h>
36 #include <libgen.h>
37 
38 #include <gsl/gsl_errno.h>
39 #include <gsl/gsl_bspline.h>
40 #include <gsl/gsl_blas.h>
41 #include <gsl/gsl_min.h>
42 #include <gsl/gsl_spline.h>
43 #include <lal/Units.h>
44 #include <lal/SeqFactories.h>
45 #include <lal/LALConstants.h>
46 #include <lal/XLALError.h>
47 #include <lal/FrequencySeries.h>
48 #include <lal/Date.h>
49 #include <lal/StringInput.h>
50 #include <lal/LALStdio.h>
51 #include <lal/FileIO.h>
52 
53 #include <lal/LALSimInspiral.h>
54 #include <lal/LALSimIMR.h>
55 #include "LALSimIMRSEOBNRROMUtilities.c"
56 
57 #include <lal/LALConfig.h>
58 #ifdef LAL_PTHREAD_LOCK
59 #include <pthread.h>
60 #endif
61 
62 
63 /******* B-spline knots over the parameter space *******/
64 static const double eta_pts[] = {0.01, 0.02, 0.03, 0.04, 0.05, 0.08264, 0.1, 0.16, 0.2, 0.25};
65 static const double chi_pts[] = {-1, -0.5, 0, 0.5, 0.99};
66 static const double fmin_pts[] = {10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, \
67 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, \
68 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, \
69 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, \
70 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, \
71 95, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, \
72 110, 111, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, \
73 124, 125, 126, 127, 128, 129, 130, 131, 132, 133, 134, 135, 136, 137, \
74 138, 139, 140, 141, 142, 143, 144, 145, 146, 147, 148, 149, 150, 151, \
75 152, 153, 154, 155, 156, 157, 158, 159, 160, 161, 162, 163, 164, 165, \
76 166, 167, 168, 169, 170, 171, 172, 173, 174, 175, 176, 177, 178, 179, \
77 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, \
78 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207, \
79 208, 209, 210, 211, 212, 213, 214, 215, 216, 217, 218, 219, 220, 221, \
80 222, 223, 224, 225, 226, 227, 228, 229, 230, 231, 232, 233, 234, 235, \
81 236, 237, 238, 239, 240, 241, 242, 243, 244, 245, 246, 247, 248, 249, \
82 250, 251, 252, 253, 254, 255, 256, 257, 258, 259, 260, 261, 262, 263, \
83 264, 265, 266, 267, 268, 269, 270, 271, 272, 273, 274, 275, 276, 277, \
84 278, 279, 280, 281, 282, 283, 284, 285, 286, 287, 288, 289, 290, 291, \
85 292, 293, 294, 295, 296, 297, 298, 299, 300, 301, 302, 303, 304, 305, \
86 306, 307, 308, 309, 310, 311, 312, 313, 314, 315, 316, 317, 318, 319, \
87 320, 321, 322, 323, 324, 325, 326, 327, 328, 329, 330, 331, 332, 333, \
88 334, 335, 336, 337, 338, 339, 340, 341, 342, 343, 344, 345, 346, 347, \
89 348, 349, 350, 351, 352, 353, 354, 355, 356, 357, 358, 359, 360, 361, \
90 362, 363, 364, 365, 366, 367, 368, 369, 370, 371, 372, 373, 374, 375, \
91 376, 377, 378, 379, 380, 381, 382, 383, 384, 385, 386, 387, 388, 389, \
92 390, 391, 392, 393, 394, 395, 396, 397, 398, 399, 400, 401, 402, 403, \
93 404, 405, 406, 407, 408, 409, 410, 411, 412, 413, 414, 415, 416, 417, \
94 418, 419, 420, 421, 422, 423, 424, 425, 426, 427, 428, 429, 430, 431, \
95 432, 433, 434, 435, 436, 437, 438, 439, 440, 441, 442, 443, 444, 445, \
96 446, 447, 448, 449, 450, 451, 452, 453, 454, 455, 456, 457, 458, 459, \
97 460, 461, 462, 463, 464, 465, 466, 467, 468, 469, 470, 471, 472, 473, \
98 474, 475, 476, 477, 478, 479, 480, 481, 482, 483, 484, 485, 486, 487, \
99 488, 489, 490, 491, 492, 493, 494, 495, 496, 497, 498, 499, 500, 501, \
100 502, 503, 504, 505, 506, 507, 508, 509, 510, 511, 512, 513, 514, 515, \
101 516, 517, 518, 519, 520, 521, 522, 523, 524, 525, 526, 527, 528, 529, \
102 530, 531, 532, 533, 534, 535, 536, 537, 538, 539, 540, 541, 542, 543, \
103 544, 545, 546, 547, 548, 549, 550, 551, 552, 553, 554, 555, 556, 557, \
104 558, 559, 560, 561, 562, 563, 564, 565, 566, 567, 568, 569, 570, 571, \
105 572, 573, 574, 575, 576, 577, 578, 579, 580, 581, 582, 583, 584, 585, \
106 586, 587, 588, 589, 590, 591, 592, 593, 594, 595, 596, 597, 598, 599, \
107 600, 601, 602, 603, 604, 605, 606, 607, 608, 609, 610, 611, 612, 613, \
108 614, 615, 616, 617, 618, 619, 620, 621, 622, 623, 624, 625, 626, 627, \
109 628, 629, 630, 631, 632, 633, 634, 635, 636, 637, 638, 639, 640, 641, \
110 642, 643, 644, 645, 646, 647, 648, 649, 650, 651, 652, 653, 654, 655, \
111 656, 657, 658, 659, 660, 661, 662, 663, 664, 665, 666, 667, 668, 669, \
112 670, 671, 672, 673, 674, 675, 676, 677, 678, 679, 680, 681, 682, 683, \
113 684, 685, 686, 687, 688, 689, 690, 691, 692, 693, 694, 695, 696, 697, \
114 698, 699, 700, 701, 702, 703, 704, 705, 706, 707, 708, 709, 710, 711, \
115 712, 713, 714, 715, 716, 717, 718, 719, 720, 721, 722, 723, 724, 725, \
116 726, 727, 728, 729, 730, 731, 732, 733, 734, 735, 736, 737, 738, 739, \
117 740, 741, 742, 743, 744, 745, 746, 747, 748, 749, 750, 751, 752, 753, \
118 754, 755, 756, 757, 758, 759, 760, 761, 762, 763, 764, 765, 766, 767, \
119 768, 769, 770, 771, 772, 773, 774, 775, 776, 777, 778, 779, 780, 781, \
120 782, 783, 784, 785, 786, 787, 788, 789, 790, 791, 792, 793, 794, 795, \
121 796, 797, 798, 799, 800, 801, 802, 803, 804, 805, 806, 807, 808, 809, \
122 810, 811, 812, 813, 814, 815, 816, 817, 818, 819, 820, 821, 822, 823, \
123 824, 825, 826, 827, 828, 829, 830, 831, 832, 833, 834, 835, 836, 837, \
124 838, 839, 840, 841, 842, 843, 844, 845, 846, 847, 848, 849, 850, 851, \
125 852, 853, 854, 855, 856, 857, 858, 859, 860, 861, 862, 863, 864, 865, \
126 866, 867, 868, 869, 870, 871, 872, 873, 874, 875, 876, 877, 878, 879, \
127 880, 881, 882, 883, 884, 885, 886, 887, 888, 889, 890, 891, 892, 893, \
128 894, 895, 896, 897, 898, 899, 900, 901, 902, 903, 904, 905, 906, 907, \
129 908, 909, 910, 911, 912, 913, 914, 915, 916, 917, 918, 919, 920, 921, \
130 922, 923, 924, 925, 926, 927, 928, 929, 930, 931, 932, 933, 934, 935, \
131 936, 937, 938, 939, 940, 941, 942, 943, 944, 945, 946, 947, 948, 949, \
132 950, 951, 952, 953, 954, 955, 956, 957, 958, 959, 960, 961, 962, 963, \
133 964, 965, 966, 967, 968, 969, 970, 971, 972, 973, 974, 975, 976, 977, \
134 978, 979, 980, 981, 982, 983, 984, 985, 986, 987, 988, 989, 990, 991, \
135 992, 993, 994, 995, 996, 997, 998, 999, 1000, 1001, 1002, 1003, 1004, \
136 1005, 1006, 1007, 1008, 1009, 1010, 1011, 1012, 1013, 1014, 1015, \
137 1016, 1017, 1018, 1019, 1020, 1021, 1022, 1023, 1024, 1025, 1026, \
138 1027, 1028, 1029, 1030, 1031, 1032, 1033, 1034, 1035, 1036, 1037, \
139 1038, 1039, 1040, 1041, 1042, 1043, 1044, 1045, 1046, 1047, 1048, \
140 1049, 1050, 1051, 1052, 1053, 1054, 1055, 1056, 1057, 1058, 1059, \
141 1060, 1061, 1062, 1063, 1064, 1065, 1066, 1067, 1068, 1069, 1070, \
142 1071, 1072, 1073, 1074, 1075, 1076, 1077, 1078, 1079, 1080, 1081, \
143 1082, 1083, 1084, 1085, 1086, 1087, 1088, 1089, 1090, 1091, 1092, \
144 1093, 1094, 1095, 1096, 1097, 1098, 1099, 1100, 1101, 1102, 1103, \
145 1104, 1105, 1106, 1107, 1108, 1109, 1110, 1111, 1112, 1113, 1114, \
146 1115, 1116, 1117, 1118, 1119, 1120, 1121, 1122, 1123, 1124, 1125, \
147 1126, 1127, 1128, 1129, 1130, 1131, 1132, 1133, 1134, 1135, 1136, \
148 1137, 1138, 1139, 1140, 1141, 1142, 1143, 1144, 1145, 1146, 1147, \
149 1148, 1149, 1150, 1151, 1152, 1153, 1154, 1155, 1156, 1157, 1158, \
150 1159, 1160, 1161, 1162, 1163, 1164, 1165, 1166, 1167, 1168, 1169, \
151 1170, 1171, 1172, 1173, 1174, 1175, 1176, 1177, 1178, 1179, 1180, \
152 1181, 1182, 1183, 1184, 1185, 1186, 1187, 1188, 1189, 1190, 1191, \
153 1192, 1193, 1194, 1195, 1196, 1197, 1198, 1199, 1200, 1201, 1202, \
154 1203, 1204, 1205, 1206, 1207, 1208, 1209, 1210, 1211, 1212, 1213, \
155 1214, 1215, 1216, 1217, 1218, 1219, 1220, 1221, 1222, 1223, 1224, \
156 1225, 1226, 1227, 1228, 1229, 1230, 1231, 1232, 1233, 1234, 1235, \
157 1236, 1237, 1238, 1239, 1240, 1241, 1242, 1243, 1244, 1245, 1246, \
158 1247, 1248, 1249, 1250, 1251, 1252, 1253, 1254, 1255, 1256, 1257, \
159 1258, 1259, 1260, 1261, 1262, 1263, 1264, 1265, 1266, 1267, 1268, \
160 1269, 1270, 1271, 1272, 1273, 1274, 1275, 1276, 1277, 1278, 1279, \
161 1280, 1281, 1282, 1283, 1284, 1285, 1286, 1287, 1288, 1289, 1290, \
162 1291, 1292, 1293, 1294, 1295, 1296, 1297, 1298, 1299, 1300, 1301, \
163 1302, 1303, 1304, 1305, 1306, 1307, 1308, 1309, 1310, 1311, 1312, \
164 1313, 1314, 1315, 1316, 1317, 1318, 1319, 1320, 1321, 1322, 1323, \
165 1324, 1325, 1326, 1327, 1328, 1329, 1330, 1331, 1332, 1333, 1334, \
166 1335, 1336, 1337, 1338, 1339, 1340, 1341, 1342, 1343, 1344, 1345, \
167 1346, 1347, 1348, 1349, 1350, 1351, 1352, 1353, 1354, 1355, 1356, \
168 1357, 1358, 1359, 1360, 1361, 1362, 1363, 1364, 1365, 1366, 1367, \
169 1368, 1369, 1370, 1371, 1372, 1373, 1374, 1375, 1376, 1377, 1378, \
170 1379, 1380, 1381, 1382, 1383, 1384, 1385, 1386, 1387, 1388, 1389, \
171 1390, 1391, 1392, 1393, 1394, 1395, 1396, 1397, 1398, 1399, 1400, \
172 1401, 1402, 1403, 1404, 1405, 1406, 1407, 1408, 1409, 1410, 1411, \
173 1412, 1413, 1414, 1415, 1416, 1417, 1418, 1419, 1420, 1421, 1422, \
174 1423, 1424, 1425, 1426, 1427, 1428, 1429, 1430, 1431, 1432, 1433, \
175 1434, 1435, 1436, 1437, 1438, 1439, 1440, 1441, 1442, 1443, 1444, \
176 1445, 1446, 1447, 1448, 1449, 1450, 1451, 1452, 1453, 1454, 1455, \
177 1456, 1457, 1458, 1459, 1460, 1461, 1462, 1463, 1464, 1465, 1466, \
178 1467, 1468, 1469, 1470, 1471, 1472, 1473, 1474, 1475, 1476, 1477, \
179 1478, 1479, 1480, 1481, 1482, 1483, 1484, 1485, 1486, 1487, 1488, \
180 1489, 1490, 1491, 1492, 1493, 1494, 1495, 1496, 1497, 1498, 1499, \
181 1500, 1501, 1502, 1503, 1504, 1505, 1506, 1507, 1508, 1509, 1510, \
182 1511, 1512, 1513, 1514, 1515, 1516, 1517, 1518, 1519, 1520, 1521, \
183 1522, 1523, 1524, 1525, 1526, 1527, 1528, 1529, 1530, 1531, 1532, \
184 1533, 1534, 1535, 1536, 1537, 1538, 1539, 1540, 1541, 1542, 1543, \
185 1544, 1545, 1546, 1547, 1548, 1549, 1550, 1551, 1552, 1553, 1554, \
186 1555, 1556, 1557, 1558, 1559, 1560, 1561, 1562, 1563, 1564, 1565, \
187 1566, 1567, 1568, 1569, 1570, 1571, 1572, 1573, 1574, 1575, 1576, \
188 1577, 1578, 1579, 1580, 1581, 1582, 1583, 1584, 1585, 1586, 1587, \
189 1588, 1589, 1590, 1591, 1592, 1593, 1594, 1595, 1596, 1597, 1598, \
190 1599, 1600, 1601, 1602, 1603, 1604, 1605, 1606, 1607, 1608, 1609, \
191 1610, 1611, 1612, 1613, 1614, 1615, 1616, 1617, 1618, 1619, 1620, \
192 1621, 1622, 1623, 1624, 1625, 1626, 1627, 1628, 1629, 1630, 1631, \
193 1632, 1633, 1634, 1635, 1636, 1637, 1638, 1639, 1640, 1641, 1642, \
194 1643, 1644, 1645, 1646, 1647, 1648, 1649, 1650, 1651, 1652, 1653, \
195 1654, 1655, 1656, 1657, 1658, 1659, 1660, 1661, 1662, 1663, 1664, \
196 1665, 1666, 1667, 1668, 1669, 1670, 1671, 1672, 1673, 1674, 1675, \
197 1676, 1677, 1678, 1679, 1680, 1681, 1682, 1683, 1684, 1685, 1686, \
198 1687, 1688, 1689, 1690, 1691, 1692, 1693, 1694, 1695, 1696, 1697, \
199 1698, 1699, 1700, 1701, 1702, 1703, 1704, 1705, 1706, 1707, 1708, \
200 1709, 1710, 1711, 1712, 1713, 1714, 1715, 1716, 1717, 1718, 1719, \
201 1720, 1721, 1722, 1723, 1724, 1725, 1726, 1727, 1728, 1729, 1730, \
202 1731, 1732, 1733, 1734, 1735, 1736, 1737, 1738, 1739, 1740, 1741, \
203 1742, 1743, 1744, 1745, 1746, 1747, 1748, 1749, 1750, 1751, 1752, \
204 1753, 1754, 1755, 1756, 1757, 1758, 1759, 1760, 1761, 1762, 1763, \
205 1764, 1765, 1766, 1767, 1768, 1769, 1770, 1771, 1772, 1773, 1774, \
206 1775, 1776, 1777, 1778, 1779, 1780, 1781, 1782, 1783, 1784, 1785, \
207 1786, 1787, 1788, 1789, 1790, 1791, 1792, 1793, 1794, 1795, 1796, \
208 1797, 1798, 1799, 1800, 1801, 1802, 1803, 1804, 1805, 1806, 1807, \
209 1808, 1809, 1810, 1811, 1812, 1813, 1814, 1815, 1816, 1817, 1818, \
210 1819, 1820, 1821, 1822, 1823};
211 
212 static const int ncx_eta = 10+2; // points in eta + 2
213 static const int ncy_chi = 5+2; // points in chi + 2
214 static const int ncz_fmin = 1814+2; // points in fmin + 2
215 
216 static const int N = 152544; // Length of flattened chirp time data vector
217 static const double Mtot0 = 12; // Total mass used for interpolant
218 static gsl_vector *ChirpTimeData = NULL;
219 
220 #ifdef LAL_PTHREAD_LOCK
221 static pthread_once_t ChirpTime_is_initialized = PTHREAD_ONCE_INIT;
222 #endif
223 
224 /*************** type definitions ******************/
225 
226 typedef struct tagSplineData
227 {
228  gsl_bspline_workspace *bwx;
229  gsl_bspline_workspace *bwy;
230  gsl_bspline_workspace *bwz;
231 } SplineData;
232 
233 static void Init_LALDATA(void);
234 
235 static int SplineData_Init(
236  SplineData **splinedata,
237  int ncx, // Number of points in x + 2
238  int ncy, // Number of points in y + 2
239  int ncz, // Number of points in z + 2
240  const double *xec, // B-spline knots in x
241  const double *yvec, // B-spline knots in y
242  const double *zvec // B-spline knots in z
243 );
244 
245 static void SplineData_Destroy(SplineData *splinedata);
246 
247 static REAL8 TP_Spline_interpolation_3d(
248  REAL8 x, // Input: x-value for which coefficient should be evaluated
249  REAL8 y, // Input: y-value for which coefficient should be evaluated
250  REAL8 z, // Input: z-value for which coefficient should be evaluated
251  gsl_vector *cvec, // Input: data for spline coefficients
252  int ncx, // Number of points in x + 2
253  int ncy, // Number of points in y + 2
254  int ncz, // Number of points in z + 2
255  const double *xvec, // B-spline knots in x
256  const double *yvec, // B-spline knots in y
257  const double *zvec // B-spline knots in z
258 );
259 
260 
261 /**************** Internal functions **********************/
262 
263 // Setup B-spline basis functions for given points
264 static int SplineData_Init(
265  SplineData **splinedata,
266  int ncx, // Number of points in x + 2
267  int ncy, // Number of points in y + 2
268  int ncz, // Number of points in z + 2
269  const double *xvec, // B-spline knots in x
270  const double *yvec, // B-spline knots in y
271  const double *zvec // B-spline knots in z
272 )
273 {
274  if(!splinedata) XLAL_ERROR(XLAL_EFAULT);
275  if(*splinedata) SplineData_Destroy(*splinedata);
276 
277  (*splinedata)=XLALCalloc(1,sizeof(SplineData));
278 
279  // Set up B-spline basis for desired knots
280  const size_t nbreak_x = ncx-2; // must have nbreak = n-2 for cubic splines
281  const size_t nbreak_y = ncy-2; // must have nbreak = n-2 for cubic splines
282  const size_t nbreak_z = ncz-2; // must have nbreak = n-2 for cubic splines
283 
284  // Allocate a cubic bspline workspace (k = 4)
285  gsl_bspline_workspace *bwx = gsl_bspline_alloc(4, nbreak_x);
286  gsl_bspline_workspace *bwy = gsl_bspline_alloc(4, nbreak_y);
287  gsl_bspline_workspace *bwz = gsl_bspline_alloc(4, nbreak_z);
288 
289  // Set breakpoints (and thus knots by hand)
290  gsl_vector *breakpts_x = gsl_vector_alloc(nbreak_x);
291  gsl_vector *breakpts_y = gsl_vector_alloc(nbreak_y);
292  gsl_vector *breakpts_z = gsl_vector_alloc(nbreak_z);
293  for (UINT4 i=0; i<nbreak_x; i++)
294  gsl_vector_set(breakpts_x, i, xvec[i]);
295  for (UINT4 j=0; j<nbreak_y; j++)
296  gsl_vector_set(breakpts_y, j, yvec[j]);
297  for (UINT4 k=0; k<nbreak_z; k++)
298  gsl_vector_set(breakpts_z, k, zvec[k]);
299 
300  gsl_bspline_knots(breakpts_x, bwx);
301  gsl_bspline_knots(breakpts_y, bwy);
302  gsl_bspline_knots(breakpts_z, bwz);
303 
304  gsl_vector_free(breakpts_x);
305  gsl_vector_free(breakpts_y);
306  gsl_vector_free(breakpts_z);
307 
308  (*splinedata)->bwx=bwx;
309  (*splinedata)->bwy=bwy;
310  (*splinedata)->bwz=bwz;
311 
312  return XLAL_SUCCESS;
313 }
314 
315 static void SplineData_Destroy(SplineData *splinedata) {
316  if(!splinedata) return;
317  if(splinedata->bwx) gsl_bspline_free(splinedata->bwx);
318  if(splinedata->bwy) gsl_bspline_free(splinedata->bwy);
319  if(splinedata->bwz) gsl_bspline_free(splinedata->bwz);
320  XLALFree(splinedata);
321 }
322 
323 // Interpolate coefficients for amplitude and phase over the parameter space.
324 // The multi-dimensional interpolation is carried out via a tensor product decomposition.
325 static REAL8 TP_Spline_interpolation_3d(
326  REAL8 x, // Input: x-value for which coefficient should be evaluated
327  REAL8 y, // Input: y-value for which coefficient should be evaluated
328  REAL8 z, // Input: z-value for which coefficient should be evaluated
329  gsl_vector *cvec, // Input: data for spline coefficients
330  int ncx, // Number of points in x + 2
331  int ncy, // Number of points in y + 2
332  int ncz, // Number of points in z + 2
333  const double *xvec, // B-spline knots in x
334  const double *yvec, // B-spline knots in y
335  const double *zvec // B-spline knots in z
336 ) {
337  SplineData *splinedata=NULL;
338  int ret = SplineData_Init(&splinedata, ncx, ncy, ncz, xvec, yvec, zvec);
339  if (ret != XLAL_SUCCESS) XLAL_ERROR(ret, "SplineData_Init() failed.");
340 
341  gsl_bspline_workspace *bwx=splinedata->bwx;
342  gsl_bspline_workspace *bwy=splinedata->bwy;
343  gsl_bspline_workspace *bwz=splinedata->bwz;
344 
345  // Evaluate the TP spline
346  REAL8 c = Interpolate_Coefficent_Tensor(cvec, x, y, z, ncy, ncz, bwx, bwy, bwz);
347 
348  SplineData_Destroy(splinedata);
349 
350  return(c);
351 }
352 
353 /** Load data file installed in $LAL_DATA_PATH
354  */
355 static void Init_LALDATA(void)
356 {
357  if (ChirpTimeData != NULL) return;
358 
359  char datafile[] = "SEOBNRv2ChirpTimeSS.dat";
360  char *path = XLAL_FILE_RESOLVE_PATH(datafile);
361  if (path==NULL)
362  XLAL_ERROR_VOID(XLAL_EIO, "Unable to resolve data file %s in $LAL_DATA_PATH\n", datafile);
363  char *dir = dirname(path);
364  ChirpTimeData = gsl_vector_alloc(N);
365  int ret = read_vector(dir, datafile, ChirpTimeData);
366  XLALFree(path);
367 
368  if(ret != XLAL_SUCCESS) {
369  gsl_vector_free(ChirpTimeData);
370  ChirpTimeData = NULL;
371  XLAL_ERROR_VOID(ret, "Unable to find data file %s in $LAL_DATA_PATH\n", datafile);
372  }
373 }
374 
375 /**
376  * @addtogroup LALSimIMRSEOBNRv2ChirpTime_c
377  *
378  * @author Michael Puerrer
379  *
380  * @brief C code for SEOBNRv2 chirp time interpolant
381  *
382  * Parameter ranges:
383  * 0.01 <= eta <= 0.25
384  * -1 <= chi <= 0.99
385  * Mtot >= 12Msun
386  * 10 <= fmin <= 1823
387  *
388  * Download the data file SEOBNRv2ChirpTimeSS.dat from
389  * https://dcc.ligo.org/LIGO-G1500097 and point LAL_DATA_PATH to it.
390  *
391  * @{
392  */
393 
394 /**
395  * Compute SEOBNRv2 chirp time from an interpolant assuming a single-spin.
396  *
397  * The chirp time is currently measured from the starting frequency f_min
398  * until the location of the amplitude peak in time.
399  */
400 REAL8 XLALSimIMRSEOBNRv2ChirpTimeSingleSpin(
401  const REAL8 m1_SI, /**< Mass of companion 1 [kg] */
402  const REAL8 m2_SI, /**< Mass of companion 2 [kg] */
403  const REAL8 chi_in, /**< Effective aligned spin */
404  const REAL8 f_min /**< Starting frequency [Hz] */
405 ) {
406  const REAL8 m1 = m1_SI / LAL_MSUN_SI;
407  const REAL8 m2 = m2_SI / LAL_MSUN_SI;
408  const REAL8 Mtot = m1 + m2;
409  REAL8 eta = m1*m2 / (Mtot*Mtot);
410  REAL8 chi = chi_in;
411 
412  // 'Nudge' parameter values to allowed boundary values if close by
413  nudge(&eta, 0.25, 1e-6);
414  nudge(&eta, 0.01, 1e-6);
415  nudge(&chi, -1.0, 1e-6);
416  nudge(&chi, 0.99, 1e-6);
417 
418  XLAL_PRINT_INFO("XLALSimIMRSEOBNRv2ChirpTimeSingleSpin(): (Mtot / Mtot0) * f_min = %g\n", (Mtot / Mtot0) * f_min);
419 
420  /* Check input domains */
421  if (eta < eta_pts[0] || eta > eta_pts[ncx_eta-3])
422  XLAL_ERROR(XLAL_EDOM);
423  if (chi < chi_pts[0] || chi > chi_pts[ncy_chi-3])
424  XLAL_ERROR(XLAL_EDOM);
425  /* Check whether total mass and f_min are inside the domain: Mtot0/Mtot fmin_min <= fmin <= Mtot0/Mtot fmin_max */
426  // For the moment hard cutoff at max(fmin_pts)
427  if (f_min * Mtot/Mtot0 > fmin_pts[ncz_fmin-3])
428  XLAL_ERROR(XLAL_EDOM, "XLALSimIMRSEOBNRv2ChirpTimeSingleSpin(): f_min * Mtot/Mtot0 = %g is outside of allowed range [%g, %g]\n", f_min * Mtot/Mtot0, fmin_pts[0], fmin_pts[ncz_fmin-3]);
429 
430  /* This is too low for us, let's fall back to PN */
431  if (f_min * Mtot/Mtot0 < fmin_pts[0]) {
432  XLAL_PRINT_WARNING("XLALSimIMRSEOBNRv2ChirpTimeSingleSpin(): f_min * Mtot/Mtot0 = %g < %g\n", f_min * Mtot/Mtot0, fmin_pts[0]);
433  XLAL_PRINT_WARNING("XLALSimIMRSEOBNRv2ChirpTimeSingleSpin(): Total mass or frequency is too low: Falling back to XLALSimInspiralTaylorF2ReducedSpinChirpTime().\n");
434  return XLALSimInspiralTaylorF2ReducedSpinChirpTime(f_min, m1_SI, m2_SI, chi, -1);
435  }
436 
437 #ifdef LAL_PTHREAD_LOCK
438  (void) pthread_once(&ChirpTime_is_initialized, Init_LALDATA);
439 #else
440  Init_LALDATA();
441 #endif
442 
443  return (Mtot / Mtot0) * TP_Spline_interpolation_3d(
444  eta,
445  chi,
446  (Mtot / Mtot0) * f_min,
447  ChirpTimeData,
448  ncx_eta,
449  ncy_chi,
450  ncz_fmin,
451  eta_pts,
452  chi_pts,
453  fmin_pts
454  );
455 }
456 
457 /** @} */
nudge
static void nudge(REAL8 *x, REAL8 X, REAL8 epsilon)
Definition: LALSimIMRPhenomP.c:1912
LALSimIMRSEOBNRROMUtilities.c
Auxiliary functions for SEOBNRv1/v2 reduced order modeling codes in LALSimIMRSEOBNRv2ROMDoubleSpin....
Interpolate_Coefficent_Tensor
static UNUSED REAL8 Interpolate_Coefficent_Tensor(gsl_vector *v, REAL8 eta, REAL8 chi1, REAL8 chi2, int ncy, int ncz, gsl_bspline_workspace *bwx, gsl_bspline_workspace *bwy, gsl_bspline_workspace *bwz)
Definition: LALSimIMRSEOBNRROMUtilities.c:432
read_vector
static UNUSED int read_vector(const char dir[], const char fname[], gsl_vector *v)
Definition: LALSimIMRSEOBNRROMUtilities.c:120
TP_Spline_interpolation_3d
static REAL8 TP_Spline_interpolation_3d(REAL8 x, REAL8 y, REAL8 z, gsl_vector *cvec, int ncx, int ncy, int ncz, const double *xvec, const double *yvec, const double *zvec)
Definition: LALSimIMRSEOBNRv2ChirpTime.c:325
Init_LALDATA
static void Init_LALDATA(void)
Load data file installed in $LAL_DATA_PATH.
Definition: LALSimIMRSEOBNRv2ChirpTime.c:355
SplineData_Destroy
static void SplineData_Destroy(SplineData *splinedata)
Definition: LALSimIMRSEOBNRv2ChirpTime.c:315
fmin_pts
static const double fmin_pts[]
Definition: LALSimIMRSEOBNRv2ChirpTime.c:66
ncy_chi
static const int ncy_chi
Definition: LALSimIMRSEOBNRv2ChirpTime.c:213
N
static const int N
Definition: LALSimIMRSEOBNRv2ChirpTime.c:216
ncx_eta
static const int ncx_eta
Definition: LALSimIMRSEOBNRv2ChirpTime.c:212
chi_pts
static const double chi_pts[]
Definition: LALSimIMRSEOBNRv2ChirpTime.c:65
ncz_fmin
static const int ncz_fmin
Definition: LALSimIMRSEOBNRv2ChirpTime.c:214
Mtot0
static const double Mtot0
Definition: LALSimIMRSEOBNRv2ChirpTime.c:217
eta_pts
static const double eta_pts[]
Definition: LALSimIMRSEOBNRv2ChirpTime.c:64
ChirpTimeData
static gsl_vector * ChirpTimeData
Definition: LALSimIMRSEOBNRv2ChirpTime.c:218
SplineData_Init
static int SplineData_Init(SplineData **splinedata, int ncx, int ncy, int ncz, const double *xec, const double *yvec, const double *zvec)
Definition: LALSimIMRSEOBNRv2ChirpTime.c:264
c
#define c
Definition: LALSimInspiralHGimri.c:35
i
double i
Definition: bh_ringdown.c:118
e
double e
Definition: bh_ringdown.c:117
XLAL_FILE_RESOLVE_PATH
#define XLAL_FILE_RESOLVE_PATH(fname)
LAL_MSUN_SI
#define LAL_MSUN_SI
REAL8
double REAL8
UINT4
uint32_t UINT4
XLALCalloc
void * XLALCalloc(size_t m, size_t n)
XLALFree
void XLALFree(void *p)
XLALSimIMRSEOBNRv2ChirpTimeSingleSpin
REAL8 XLALSimIMRSEOBNRv2ChirpTimeSingleSpin(const REAL8 m1_SI, const REAL8 m2_SI, const REAL8 chi_in, const REAL8 f_min)
Compute SEOBNRv2 chirp time from an interpolant assuming a single-spin.
Definition: LALSimIMRSEOBNRv2ChirpTime.c:400
XLALSimInspiralTaylorF2ReducedSpinChirpTime
REAL8 XLALSimInspiralTaylorF2ReducedSpinChirpTime(const REAL8 fStart, const REAL8 m1_SI, const REAL8 m2_SI, const REAL8 chi, const INT4 O)
Compute the chirp time of the "reduced-spin" templates.
Definition: LALSimInspiralTaylorF2ReducedSpin.c:307
XLAL_ERROR_VOID
#define XLAL_ERROR_VOID(...)
XLAL_PRINT_INFO
#define XLAL_PRINT_INFO(...)
XLAL_ERROR
#define XLAL_ERROR(...)
XLAL_PRINT_WARNING
#define XLAL_PRINT_WARNING(...)
XLAL_SUCCESS
XLAL_SUCCESS
XLAL_EFAULT
XLAL_EFAULT
XLAL_EDOM
XLAL_EDOM
XLAL_EIO
XLAL_EIO
x
list x
y
list y
path
path
pthread_once
int pthread_once(pthread_once_t *once_control, void(*init_routine)(void))
SplineData::bwx
gsl_bspline_workspace * bwx
Definition: LALSimIMRSEOBNRv1ROMDoubleSpin.c:139
SplineData::bwz
gsl_bspline_workspace * bwz
Definition: LALSimIMRSEOBNRv1ROMDoubleSpin.c:141
SplineData::bwy
gsl_bspline_workspace * bwy
Definition: LALSimIMRSEOBNRv1ROMDoubleSpin.c:140
f_min
double f_min
Definition: unicorn.c:22