Skip to content
Snippets Groups Projects

Additional sanity checks and renormalization in NR waveform generation.

Merged Ian Harry requested to merge ian-harry/lalsuite:nrwaveform_accuracy_fix into master
1 file
+ 43
2
Compare changes
  • Side-by-side
  • Inline
@@ -361,8 +361,8 @@ UNUSED static UINT4 XLALSimInspiralNRWaveformGetRotationAnglesFromH5File(
REAL8 n_hat_x, n_hat_y, n_hat_z, n_hat_norm;
REAL8 pos1x, pos1y, pos1z, pos2x, pos2y, pos2z;
REAL8 corb_phase, sorb_phase, sinclination, cinclination;
REAL8 ln_cross_n_x, ln_cross_n_y, ln_cross_n_z;
REAL8 z_wave_x, z_wave_y, z_wave_z;
REAL8 ln_cross_n_x, ln_cross_n_y, ln_cross_n_z, ln_cross_n_norm;
REAL8 z_wave_x, z_wave_y, z_wave_z, z_wave_norm;
REAL8 stheta, ctheta, spsi, cpsi, theta_hat_x, theta_hat_y, theta_hat_z;
REAL8 psi_hat_x, psi_hat_y, psi_hat_z;
REAL8 n_dot_theta, ln_cross_n_dot_theta, n_dot_psi, ln_cross_n_dot_psi;
@@ -420,12 +420,24 @@ UNUSED static UINT4 XLALSimInspiralNRWaveformGetRotationAnglesFromH5File(
ln_hat_norm = sqrt(ln_hat_x * ln_hat_x + ln_hat_y * ln_hat_y + ln_hat_z * ln_hat_z);
if (fabs(ln_hat_norm - 1) > 0.001)
{
XLAL_PRINT_ERROR("ln_hat_mag = %.16f\n", ln_hat_norm);
XLAL_ERROR(XLAL_EDOM, "The size of the LN hat vector in the supplied HDF file is not equal to unity");
}
ln_hat_x = ln_hat_x / ln_hat_norm;
ln_hat_y = ln_hat_y / ln_hat_norm;
ln_hat_z = ln_hat_z / ln_hat_norm;
n_hat_norm = sqrt(n_hat_x * n_hat_x + n_hat_y * n_hat_y + n_hat_z * n_hat_z);
if (fabs(n_hat_norm - 1) > 0.001)
{
XLAL_PRINT_ERROR("n_hat_mag = %.16f\n", n_hat_norm);
XLAL_ERROR(XLAL_EDOM, "The size of the n hat vector in the supplied HDF file is not equal to unity");
}
n_hat_x = n_hat_x / n_hat_norm;
n_hat_y = n_hat_y / n_hat_norm;
n_hat_z = n_hat_z / n_hat_norm;
@@ -441,6 +453,17 @@ UNUSED static UINT4 XLALSimInspiralNRWaveformGetRotationAnglesFromH5File(
ln_cross_n_y = ln_hat_z * n_hat_x - ln_hat_x * n_hat_z;
ln_cross_n_z = ln_hat_x * n_hat_y - ln_hat_y * n_hat_x;
ln_cross_n_norm = sqrt(ln_cross_n_x*ln_cross_n_x + ln_cross_n_y*ln_cross_n_y + ln_cross_n_z*ln_cross_n_z);
ln_cross_n_x = ln_cross_n_x / ln_cross_n_norm;
ln_cross_n_y = ln_cross_n_y / ln_cross_n_norm;
ln_cross_n_z = ln_cross_n_z / ln_cross_n_norm;
if (fabs(ln_cross_n_norm - 1) > 0.001)
{
XLAL_PRINT_ERROR("ln_cross_n mag = %.16f\n", ln_cross_n_norm);
XLAL_ERROR(XLAL_EDOM, "The size of the LN x n vector computed here is not equal to unity. This shouldn't be possible, please email lalsimulation developers.");
}
z_wave_x = sinclination * (sorb_phase * n_hat_x + corb_phase * ln_cross_n_x);
z_wave_y = sinclination * (sorb_phase * n_hat_y + corb_phase * ln_cross_n_y);
z_wave_z = sinclination * (sorb_phase * n_hat_z + corb_phase * ln_cross_n_z);
@@ -449,6 +472,18 @@ UNUSED static UINT4 XLALSimInspiralNRWaveformGetRotationAnglesFromH5File(
z_wave_y += cinclination * ln_hat_y;
z_wave_z += cinclination * ln_hat_z;
z_wave_norm = sqrt(z_wave_x * z_wave_x + z_wave_y * z_wave_y + z_wave_z * z_wave_z);
z_wave_x = z_wave_x / z_wave_norm;
z_wave_y = z_wave_y / z_wave_norm;
z_wave_z = z_wave_z / z_wave_norm;
if (fabs(z_wave_norm - 1) > 0.001)
{
XLAL_PRINT_ERROR("Z mag = %.16f\n", z_wave_norm);
XLAL_ERROR(XLAL_EDOM, "The size of the Z vector computed here is not equal to unity. This shouldn't be possible, please email lalsimulation developers.");
}
/* Step 3.1: Extract theta and psi from Z in the lal wave frame
* NOTE: Theta can only run between 0 and pi, so no problem with arccos here
*/
@@ -480,6 +515,12 @@ UNUSED static UINT4 XLALSimInspiralNRWaveformGetRotationAnglesFromH5File(
*psi = 0.;
}
}
else
{
XLAL_PRINT_ERROR("z_wave_x = %.16f\n", z_wave_x);
XLAL_PRINT_ERROR("sin(theta) = %.16f\n", sin(*theta));
XLAL_ERROR(XLAL_EDOM, "Z_x cannot be bigger than sin(theta). Please email lalsimulation developers.");
}
}
else
{
Loading