Loading [MathJax]/extensions/TeX/AMSsymbols.js
LALSimulation 6.2.0.1-8a6b96f
All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages
PrecessingHlmsTest.c
Go to the documentation of this file.
1/*
2 * Copyright (C) 2018 Riccardo Sturani
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//#include <stdlib.h>
21#include <math.h>
22#include <complex.h>
23#include <lal/Date.h>
24#include <lal/Units.h>
25#include <lal/TimeSeries.h>
26#include <lal/LALConstants.h>
27#include <lal/LALStdlib.h>
28#include <lal/LALSimInspiral.h>
29#include <lal/LALSimInspiralPrecess.h>
30#include <lal/LALSimSphHarmMode.h>
31
32/*
33#include <lal/LALSimInspiral.h>
34#include <lal/LALConstants.h>
35#include <lal/LALSimInspiralWaveformFlags.h>
36#include <lal/Units.h>
37#include <lal/XLALError.h>*/
38
39#define EPSILON 1.e-8
40
41#define ROTATEZ(angle, vx, vy, vz)\
42 tmp1 = vx*cos(angle) - vy*sin(angle);\
43 tmp2 = vx*sin(angle) + vy*cos(angle);\
44 vx = tmp1;\
45 vy = tmp2
46
47#define ROTATEY(angle, vx, vy, vz)\
48 tmp1 = vx*cos(angle) + vz*sin(angle);\
49 tmp2 = - vx*sin(angle) + vz*cos(angle);\
50 vx = tmp1;\
51 vz = tmp2
52
53static int c_compare(COMPLEX16 val1,
54 COMPLEX16 val2)
55{
56 if ( (fabs(creal(val1) - creal(val2)) > EPSILON) || (fabs(cimag(val1) - cimag(val2)) > EPSILON) )
57 return 1;
58 else
59 return 0;
60}
61
62static int compare(REAL8 val1,
63 REAL8 val2)
64{
65 if ( ( (fabs(val1 - val2)) > EPSILON) > EPSILON )
66 return 1;
67 else
68 return 0;
69}
70
71int main (int argc, char **argv)
72{
73 /* Ignore unused parameters. */
74 (void)argc;
75 (void)argv;
76
77 const UINT4 Ntest=10;
78 const REAL8 NtestR=(REAL8) Ntest;
79 int errCode=0;
80 int ret=0;
81 REAL8 idxr;
82 REAL8 alpha,iota,psi,iota2,psi2,v;
83 REAL8 tmp1,tmp2,tmpx,tmpy,tmpz;
84
85 const INT4 ampO = -1;
86 const REAL8 m1 = 10.*LAL_MSUN_SI;
87 const REAL8 m2 = 5.*LAL_MSUN_SI;
88 const REAL8 eta = m1/(m1+m2)*m2/(m1+m2);
89 const REAL8 dist = LAL_PC_SI*1.e6;
90 const REAL8 norm = eta*(m1+m2)*LAL_G_SI/pow(LAL_C_SI,2.) / dist;
91 const REAL8 chi1 = 0.8;
92 const REAL8 chi2 = 0.3;
93 const REAL8 dt = 1.;
94 const LIGOTimeGPS time0=LIGOTIMEGPSZERO;
95
133 REAL8TimeSeries *hpF = XLALCreateREAL8TimeSeries("",&time0,0.,dt,&lalStrainUnit,Ntest);
134 REAL8TimeSeries *hcF = XLALCreateREAL8TimeSeries("",&time0,0.,dt,&lalStrainUnit,Ntest);
135 REAL8TimeSeries *hpFs = XLALCreateREAL8TimeSeries("",&time0,0.,dt,&lalStrainUnit,Ntest);
136 REAL8TimeSeries *hcFs = XLALCreateREAL8TimeSeries("",&time0,0.,dt,&lalStrainUnit,Ntest);
137 REAL8TimeSeries *hpFc = XLALCreateREAL8TimeSeries("",&time0,0.,dt,&lalStrainUnit,Ntest);
138 REAL8TimeSeries *hcFc = XLALCreateREAL8TimeSeries("",&time0,0.,dt,&lalStrainUnit,Ntest);
139 REAL8TimeSeries *hpFc2 = XLALCreateREAL8TimeSeries("",&time0,0.,dt,&lalStrainUnit,Ntest);
140 REAL8TimeSeries *hcFc2 = XLALCreateREAL8TimeSeries("",&time0,0.,dt,&lalStrainUnit,Ntest);
141 REAL8TimeSeries *hpFc3 = XLALCreateREAL8TimeSeries("",&time0,0.,dt,&lalStrainUnit,Ntest);
142 REAL8TimeSeries *hcFc3 = XLALCreateREAL8TimeSeries("",&time0,0.,dt,&lalStrainUnit,Ntest);
143
144 SphHarmTimeSeries *hStart, *hFin, *hFin2, *hFinRotm1, *hFinCheck, *hFinCheck2;
145 hStart=NULL;
146 hFin=NULL;
147 hFin2=NULL;
148 hFinRotm1=NULL;
149 hFinCheck=NULL;
150 hFinCheck2=NULL;
151 COMPLEX16TimeSeries *hStartm,*hFinm,*hFinRotm1m,*hFin2m,*hFinChkm,*hFinChk2m;
152
153 for (INT4 l=2;l<=4;l++) {
154
155 for (UINT4 idx=0;idx<Ntest;idx++) {
156
157 idxr=((REAL8)idx)/((REAL8) Ntest);
158 alpha = idxr*2.*LAL_PI/NtestR;
159 iota = idxr*LAL_PI/NtestR;
160 iota2 = idxr*LAL_PI/NtestR*0.4;
161 psi = (NtestR-idxr-1.)*2.*LAL_PI/NtestR;
162 psi2 = (NtestR-idxr-1.)*2.*LAL_PI/NtestR*0.4;
163 v = (idxr+1.)/NtestR;
164
165 zts->data->data[idx]=0.;
166 vts->data->data[idx]=v;
167 psits->data->data[idx]=psi;
168 psi2ts->data->data[idx]=psi2;
169 mpsits->data->data[idx]=-psi;
170 mpsi2ts->data->data[idx]=-psi2;
171 psi0ts->data->data[idx]=0.;
172 iotats->data->data[idx]=iota;
173 iota2ts->data->data[idx]=iota2;
174 alphats->data->data[idx]=alpha;
175 miotats->data->data[idx]=-iota;
176 miota2ts->data->data[idx]=-iota2;
177 malphats->data->data[idx]=-alpha;
178 LNhx0->data->data[idx]=0.;
179 LNhy0->data->data[idx]=0.;
180 LNhz0->data->data[idx]=1.;
181 e1x0->data->data[idx]=0.;
182 e1y0->data->data[idx]=1.;
183 e1z0->data->data[idx]=0.;
184 LNhx->data->data[idx]=sin(iota)*cos(alpha);
185 LNhy->data->data[idx]=sin(iota)*sin(alpha);
186 LNhz->data->data[idx]=cos(iota);
187 e1x->data->data[idx]=-sin(alpha);
188 e1y->data->data[idx]=cos(alpha);
189 e1z->data->data[idx]=0.;
190 hpF->data->data[idx]=0.;
191 hcF->data->data[idx]=0.;
192 hpFs->data->data[idx]=0.;
193 hcFs->data->data[idx]=0.;
194 hpFc->data->data[idx]=0.;
195 hcFc->data->data[idx]=0.;
196 hpFc2->data->data[idx]=0.;
197 hcFc2->data->data[idx]=0.;
198 hpFc3->data->data[idx]=0.;
199 hcFc3->data->data[idx]=0.;
200
201 S1x0->data->data[idx]=0.;
202 S1y0->data->data[idx]=0.;
203 S1z0->data->data[idx]=chi1;
204 S2x0->data->data[idx]=0.;
205 S2y0->data->data[idx]=0.;
206 S2z0->data->data[idx]=chi2;
207
208 tmpx=S1x0->data->data[idx];
209 tmpy=S1y0->data->data[idx];
210 tmpz=S1z0->data->data[idx];
211 ROTATEZ(psi, tmpx,tmpy,tmpz);
212 ROTATEY(iota, tmpx,tmpy,tmpz);
213 ROTATEZ(alpha,tmpx,tmpy,tmpz);
214 S1x->data->data[idx]=tmpx;
215 S1y->data->data[idx]=tmpy;
216 S1z->data->data[idx]=tmpz;
217
218 tmpx=S2x0->data->data[idx];
219 tmpy=S2y0->data->data[idx];
220 tmpz=S2z0->data->data[idx];
221 ROTATEZ(psi, tmpx,tmpy,tmpz);
222 ROTATEY(iota, tmpx,tmpy,tmpz);
223 ROTATEZ(alpha,tmpx,tmpy,tmpz);
224 S2x->data->data[idx]=tmpx;
225 S2y->data->data[idx]=tmpy;
226 S2z->data->data[idx]=tmpz;
227
228 }
229
230 //We create the mode with trivial and non-trivial argument
231 if (l==2) {
232 errCode += XLALSimInspiralSpinPNMode2m(&hStart,vts,psi0ts,LNhx0,LNhy0,LNhz0,e1x0,e1y0,e1z0,S1x0,S1y0,S1z0,S2x0,S2y0,S2z0,m1,m2,dist,ampO);
233 errCode += XLALSimInspiralSpinPNMode2m(&hFin,vts,psits,LNhx,LNhy,LNhz,e1x,e1y,e1z,S1x,S1y,S1z,S2x,S2y,S2z,m1,m2,dist,ampO);
234 }
235 else if (l==3) {
236 errCode += XLALSimInspiralSpinPNMode3m(&hStart,vts,psi0ts,LNhx0,LNhy0,LNhz0,e1x0,e1y0,e1z0,S1x0,S1y0,S1z0,S2x0,S2y0,S2z0,m1,m2,dist,ampO);
237 errCode += XLALSimInspiralSpinPNMode3m(&hFin,vts,psits,LNhx,LNhy,LNhz,e1x,e1y,e1z,S1x,S1y,S1z,S2x,S2y,S2z,m1,m2,dist,ampO);
238 }
239 else if (l==4) {
240 errCode += XLALSimInspiralSpinPNMode4m(&hStart,vts,psi0ts,LNhx0,LNhy0,LNhz0,e1x0,e1y0,e1z0,S1x0,S1y0,S1z0,S2x0,S2y0,S2z0,m1,m2,dist,ampO);
241 errCode += XLALSimInspiralSpinPNMode4m(&hFin,vts,psits,LNhx,LNhy,LNhz,e1x,e1y,e1z,S1x,S1y,S1z,S2x,S2y,S2z,m1,m2,dist,ampO);
242 }
243
244 //and we check that a rotation takes one into the other and back
245 errCode = XLALSimInspiralPrecessionRotateModesOut(&hFin2, hStart, psits, iotats, alphats);
246 errCode += XLALSimInspiralPrecessionRotateModesOut(&hFinRotm1, hFin, malphats, miotats, mpsits);
247 //Two more checks with additional rotations
248 errCode += XLALSimInspiralPrecessionRotateModesOut(&hFinCheck, hFin, zts, iota2ts, psi2ts);
249 errCode += XLALSimInspiralPrecessionRotateModesOut(&hFinCheck2, hFin, psi2ts, iota2ts, zts);
250
251 for(INT4 m=-l;m<=l;m++) {
252 hStartm = XLALSphHarmTimeSeriesGetMode(hStart, l, m );
253 hFinRotm1m = XLALSphHarmTimeSeriesGetMode(hFinRotm1, l, m );
254 hFinm = XLALSphHarmTimeSeriesGetMode(hFin, l, m );
255 hFin2m = XLALSphHarmTimeSeriesGetMode(hFin2, l, m );
256 hFinChkm = XLALSphHarmTimeSeriesGetMode(hFinCheck, l, m );
257 hFinChk2m = XLALSphHarmTimeSeriesGetMode(hFinCheck2, l, m );
258
259 errCode += XLALSimAddMode(hpF, hcF, hFinm, 0., 0., l, m, 0);
260 errCode += XLALSimAddModeAngleTimeSeries(hpFs, hcFs, hStartm, miotats, mpsits, l, m, 0);
261 errCode += XLALSimAddModeAngleTimeSeries(hpFc, hcFc, hFinChkm, iota2ts, psi2ts, l, m, 0);
262 errCode += XLALSimAddMode(hpFc2, hcFc2, hFinChk2m, 0., 0., l, m, 0);
263 errCode += XLALSimAddModeAngleTimeSeries(hpFc3, hcFc3, hFinm, miota2ts , mpsi2ts, l, m, 0);
264
265 for(UINT4 idx=0;idx<Ntest;idx++) {
266 ret+=c_compare(hStartm->data->data[idx]/norm,hFinRotm1m->data->data[idx]/norm);
267 ret+=c_compare(hFinm->data->data[idx]/norm,hFin2m->data->data[idx]/norm);
268 }
269
271 hStartm=NULL;
273 hFinRotm1m=NULL;
275 hFinm=NULL;
277 hFinChkm=NULL;
279 hFinChk2m=NULL;
280 }
281
282 hStart=NULL;
283 hFin=NULL;
284 hFin2=NULL;
285 hFinRotm1=NULL;
286 hFinCheck=NULL;
287 hFinCheck2=NULL;
288
289 for(UINT4 idx=0;idx<Ntest;idx++) {
290 ret+=compare(hpF->data->data[idx]/norm,(cos(2.*alphats->data->data[idx])*hpFs->data->data[idx]-sin(2.*alphats->data->data[idx])*hcFs->data->data[idx])/norm);
291 ret+=compare(hcF->data->data[idx]/norm,(cos(2.*alphats->data->data[idx])*hcFs->data->data[idx]+sin(2.*alphats->data->data[idx])*hpFs->data->data[idx])/norm);
292 ret+=compare(hpF->data->data[idx]/norm,hpFc->data->data[idx]/norm);
293 ret+=compare(hcF->data->data[idx]/norm,hcFc->data->data[idx]/norm);
294 ret+=compare(hpFc2->data->data[idx]/norm,hpFc3->data->data[idx]/norm);
295 ret+=compare(hcFc2->data->data[idx]/norm,hcFc3->data->data[idx]/norm);
296 }
297
298 }
299
300 if ( (ret == 0) && (errCode == 0) ) {
301 fprintf(stdout, "\n Precessing modes test passed.\n");
302 }
303 else {
304 fprintf(stderr, "\nFAILURE: %u Precessing modes test failed.\n", ret+errCode);
305 }
306
307 return ret + errCode ;
308}
INT4 XLALSimInspiralSpinPNMode4m(SphHarmTimeSeries **hlm, REAL8TimeSeries *V, REAL8TimeSeries *Phi, REAL8TimeSeries *LNhx, REAL8TimeSeries *LNhy, REAL8TimeSeries *LNhz, REAL8TimeSeries *e1x, REAL8TimeSeries *e1y, REAL8TimeSeries *e1z, REAL8TimeSeries *S1x, REAL8TimeSeries *S1y, REAL8TimeSeries *S1z, REAL8TimeSeries *S2x, REAL8TimeSeries *S2y, REAL8TimeSeries *S2z, REAL8 m1, REAL8 m2, REAL8 distance, int ampO)
#define EPSILON
#define ROTATEZ(angle, vx, vy, vz)
int main(int argc, char **argv)
static int c_compare(COMPLEX16 val1, COMPLEX16 val2)
#define ROTATEY(angle, vx, vy, vz)
static int compare(REAL8 val1, REAL8 val2)
REAL8 tmp1
#define fprintf
int l
Definition: bh_qnmode.c:135
double dt
Definition: bh_ringdown.c:113
#define LAL_C_SI
#define LAL_MSUN_SI
#define LAL_PI
#define LAL_PC_SI
#define LAL_G_SI
#define LIGOTIMEGPSZERO
double complex COMPLEX16
double REAL8
uint32_t UINT4
int32_t INT4
INT4 XLALSimInspiralSpinPNMode2m(SphHarmTimeSeries **hlm, REAL8TimeSeries *V, REAL8TimeSeries *Phi, REAL8TimeSeries *LNhx, REAL8TimeSeries *LNhy, REAL8TimeSeries *LNhz, REAL8TimeSeries *e1x, REAL8TimeSeries *e1y, REAL8TimeSeries *e1z, REAL8TimeSeries *S1x, REAL8TimeSeries *S1y, REAL8TimeSeries *S1z, REAL8TimeSeries *S2x, REAL8TimeSeries *S2y, REAL8TimeSeries *S2z, REAL8 m1, REAL8 m2, REAL8 distance, int ampO)
Computes the 5 l=2 modes of spherical harmonic decomposition of the post-Newtonian inspiral waveform ...
INT4 XLALSimInspiralSpinPNMode3m(SphHarmTimeSeries **hlm, REAL8TimeSeries *V, REAL8TimeSeries *Phi, REAL8TimeSeries *LNhx, REAL8TimeSeries *LNhy, REAL8TimeSeries *LNhz, REAL8TimeSeries *e1x, REAL8TimeSeries *e1y, REAL8TimeSeries *e1z, REAL8TimeSeries *S1x, REAL8TimeSeries *S1y, REAL8TimeSeries *S1z, REAL8TimeSeries *S2x, REAL8TimeSeries *S2y, REAL8TimeSeries *S2z, REAL8 m1, REAL8 m2, REAL8 distance, int ampO)
Computes all l=3 modes of spherical harmonic decomposition of the post-Newtonian inspiral waveform fo...
int XLALSimInspiralPrecessionRotateModesOut(SphHarmTimeSeries **hlm_out, SphHarmTimeSeries *hlm_in, const REAL8TimeSeries *alpha, const REAL8TimeSeries *beta, const REAL8TimeSeries *gam)
Takes in the h_lm spherical harmonic decomposed modes and rotates the modes by Euler angles alpha,...
int XLALSimAddModeAngleTimeSeries(REAL8TimeSeries *hplus, REAL8TimeSeries *hcross, COMPLEX16TimeSeries *hmode, REAL8TimeSeries *theta, REAL8TimeSeries *phi, int l, int m, int sym)
Multiplies a mode h(l,m) by a spin-2 weighted spherical harmonic to obtain hplus - i hcross,...
int XLALSimAddMode(REAL8TimeSeries *hplus, REAL8TimeSeries *hcross, COMPLEX16TimeSeries *hmode, REAL8 theta, REAL8 phi, int l, int m, int sym)
Multiplies a mode h(l,m) by a spin-2 weighted spherical harmonic to obtain hplus - i hcross,...
COMPLEX16TimeSeries * XLALSphHarmTimeSeriesGetMode(SphHarmTimeSeries *ts, UINT4 l, INT4 m)
Get the time series of a waveform's (l,m) spherical harmonic mode from a SphHarmTimeSeries linked lis...
static const INT4 m
void XLALDestroyCOMPLEX16TimeSeries(COMPLEX16TimeSeries *series)
REAL8TimeSeries * XLALCreateREAL8TimeSeries(const CHAR *name, const LIGOTimeGPS *epoch, REAL8 f0, REAL8 deltaT, const LALUnit *sampleUnits, size_t length)
const LALUnit lalStrainUnit
const LALUnit lalDimensionlessUnit
double alpha
Definition: sgwb.c:106
COMPLEX16Sequence * data
COMPLEX16 * data
REAL8Sequence * data
REAL8 * data
Structure to carry a collection of spherical harmonic modes in COMPLEX16 time series.