Commit 4efb09cc authored by Duncan Macleod's avatar Duncan Macleod
Browse files

Merge branch 'o3-liv-final' into 'master'

Adding code to constrain Lorentz Invariance Violation (LIV)

See merge request lscsoft/lalsuite!1292
parents a8948d21 86dcc34b
......@@ -47,7 +47,7 @@ __version__= "git id %s"%git_version.id
__date__= git_version.date
#List of parameters to plot/bin . Need to match (converted) column names.
oneDMenu=['mtotal','m1','m2','mchirp','mc','chirpmass','distance','distMPC','dist','iota','psi','eta','q','asym_massratio','spin1','spin2','a1','a2','phi1','theta1','phi2','theta2','costilt1','costilt2','costhetas','cosbeta','phi_orb', 'lambdat', 'dlambdat', 'lambda1', 'lambda2', 'lam_tilde', 'dlam_tilde','theta_jn','a1z','a2z'] + bppu.snrParams + bppu.spinParams + bppu.cosmoParam + bppu.calParams
oneDMenu=['mtotal','m1','m2','mchirp','mc','chirpmass','distance','distMPC','dist','iota','psi','eta','q','asym_massratio','spin1','spin2','a1','a2','phi1','theta1','phi2','theta2','costilt1','costilt2','costhetas','cosbeta','phi_orb', 'lambdat', 'dlambdat', 'lambda1', 'lambda2', 'lam_tilde', 'dlam_tilde','theta_jn','a1z','a2z'] + bppu.snrParams + bppu.spinParams + bppu.cosmoParam + bppu.calParams + bppu.lorentzInvarianceViolationParams
#List of parameter pairs to bin . Need to match (converted) column names.
twoDGreedyMenu=[['mc','eta'],['mchirp','eta'],['chirpmass','eta'],['mc','q'],['mchirp','q'],['chirpmass','q'],['mc','asym_massratio'],['mchirp','asym_massratio'],['chirpmass','asym_massratio'],['m1','m2'],['mtotal','eta'],['distance','iota'],['dist','iota'],['dist','m1'],['ra','dec'],['dist','cos(iota)'],['phi_orb','iota'],['theta_jn','dist'],['spin1','spin2'],['spin1','mchirp'],['spin1','m1'],['a1','a2'],['a1','mchirp'],['a1','m1'],['tilt1','tilt2'],['tilt1','mchirp'],['tilt1','m1'],['a1z','a2z']]
#Bin size/resolution for binning. Need to match (converted) column names.
......@@ -62,10 +62,10 @@ paramNameLatexMap = {'m1': 'm_1', 'm2' : 'm_2', 'mtotal' : r'M_{\rm tot}', 'mchi
'cos(iota)': r'\cos \iota', 'tilt1': r't_1', 'tilt2': r't_2', 'ra': r'\alpha', 'dec': r'\delta',
'lambdat' : r'\tilde{\Lambda}', 'dlambdat': r'\delta \tilde{\Lambda}',
'lambda1' : r'\lambda_1', 'lambda2': r'\lambda_2',
'lam_tilde' : r'\tilde{\Lambda}', 'dlam_tilde': r'\delta \tilde{\Lambda}','dchi0':r'\delta\chi_0','dchi1':r'\delta\chi_1','dchi2':r'\delta\chi_2','dchi3':r'\delta\chi_3','dchi4':r'\delta\chi_4','dchi5':r'\delta\chi_5','dchi5l':r'\delta\chi_{5l}','dchi6':r'\delta\chi_6','dchi6l':r'\delta\chi_{6l}','dchi7':r'\delta\chi_7','dbeta2':r'\delta\beta_2','dbeta3':r'\delta\beta_3','dsigma2':r'\delta\sigma_2','dsigma3':r'\delta\sigma_3','dsigma4':r'\delta\sigma_4','dbeta2':r'\delta\beta_2','dbeta3':r'\delta\beta_3' }
'lam_tilde' : r'\tilde{\Lambda}', 'dlam_tilde': r'\delta \tilde{\Lambda}','dchi0':r'\delta\chi_0','dchi1':r'\delta\chi_1','dchi2':r'\delta\chi_2','dchi3':r'\delta\chi_3','dchi4':r'\delta\chi_4','dchi5':r'\delta\chi_5','dchi5l':r'\delta\chi_{5l}','dchi6':r'\delta\chi_6','dchi6l':r'\delta\chi_{6l}','dchi7':r'\delta\chi_7','dbeta2':r'\delta\beta_2','dbeta3':r'\delta\beta_3','dsigma2':r'\delta\sigma_2','dsigma3':r'\delta\sigma_3','dsigma4':r'\delta\sigma_4','dbeta2':r'\delta\beta_2','dbeta3':r'\delta\beta_3' ,'log10lambda_a':r'$\log\lambda_{\mathbb{A}}$','log10lambda_eff':r'$\log\lambda_{eff}$','log10livamp':r'$\log \mathbb{A}$','lambda_a':r'$\lambda_{\mathbb{A}}$','lambda_eff':r'$\lambda_{eff}$','liv_amp':r'$\mathbb{A}$'}
# Only these parameters, in this order appear in confidence level table.
clTableParams = ['mchirp', 'mc', 'chirpmass', 'eta', 'q', 'm1', 'm2', 'distance', 'distMPC', 'dist', 'cos(iota)', 'iota', 'theta_jn', 'psi', 'ra', 'dec', 'time', 'phase', 'a1', 'a2', 'costilt1', 'costilt2','dchi0','dchi1','dchi2','dchi3','dchi4','dchi5','dchi5l','dchi6','dchi6l','dchi7','dbeta2','dbeta3','dsigma2','dsigma3','dsigma4','dbeta2','dbeta3']
clTableParams = ['mchirp', 'mc', 'chirpmass', 'eta', 'q', 'm1', 'm2', 'distance', 'distMPC', 'dist', 'cos(iota)', 'iota', 'theta_jn', 'psi', 'ra', 'dec', 'time', 'phase', 'a1', 'a2', 'costilt1', 'costilt2','dchi0','dchi1','dchi2','dchi3','dchi4','dchi5','dchi5l','dchi6','dchi6l','dchi7','dbeta2','dbeta3','dsigma2','dsigma3','dsigma4','dbeta2','dbeta3', 'log10lambda_eff','log10lambda_a','log10livamp','lambda_eff','lambda_a','liv_amp']
greedyBinSizes={'mc':0.001,'m1':0.1,'m2':0.1,'mass1':0.1,'mass2':0.1,'mtotal':0.1,'eta':0.001,'q':0.001,'asym_massratio':0.001,'iota':0.05,'time':1e-4,'distance':5.0,'dist':1.0,'mchirp':0.01,'chirpmass':0.01,'a1':0.02,'a2':0.02,'phi1':0.05,'phi2':0.05,'theta1':0.05,'theta2':0.05,'ra':0.05,'dec':0.005,'psi':0.1,'cos(iota)':0.01, 'cos(tilt1)':0.01, 'cos(tilt2)':0.01, 'tilt1':0.05, 'tilt2':0.05, 'cos(thetas)':0.01, 'cos(beta)':0.01,'phi_orb':0.2,'inclination':0.05,'theta_jn':0.05,'spin1':0.02,'spin2':0.02}
......@@ -73,7 +73,7 @@ for s in bppu.snrParams:
greedyBinSizes[s]=0.02
for s in bppu.calParams:
greedyBinSizes[s]=0.02
for s in bppu.tigerParams:
for s in bppu.tigerParams + bppu.lorentzInvarianceViolationParams:
greedyBinSizes[s]=0.01
for s in bppu.spinParams:
greedyBinSizes[s]=0.02
......
......@@ -1192,6 +1192,21 @@ if __name__=='__main__':
for mp in massParams:
for dchi in bppu.tigerParams:
twoDGreedyMenu.append([mp,dchi])
for mp in massParams:
for lvp in bppu.lorentzInvarianceViolationParams:
twoDGreedyMenu.append([mp,lvp])
for sp in spinParams:
for lvp in bppu.lorentzInvarianceViolationParams:
twoDGreedyMenu.append([sp,lvp])
for dchi in bppu.tigerParams:
for lvp in bppu.lorentzInvarianceViolationParams:
twoDGreedyMenu.append([dchi, lvp])
for eg in bppu.energyParams:
for lvp in bppu.lorentzInvarianceViolationParams:
twoDGreedyMenu.append([eg, lvp])
for dp in distParams:
for lvp in bppu.lorentzInvarianceViolationParams:
twoDGreedyMenu.append([dp, lvp])
for dp in distParams:
for sp in snrParams:
twoDGreedyMenu.append([dp,sp])
......
......@@ -754,6 +754,9 @@ LALInferenceModel *LALInferenceInitCBCModel(LALInferenceRunState *state) {
(--singleSpin) template will assume only the spin of the most massive binary component exists.\n\
(--noSpin, --disable-spin) template will assume no spins (giving this will void spinOrder!=0) \n\
(--no-detector-frame) model will NOT use detector-centred coordinates and instead RA,dec\n\
(--nonGR_alpha value) this is a LIV parameter which should only be passed when log10lambda_eff/lambda_eff is passed as a grtest-parameter for LIV test\n\
(--LIV_A_sign) this is a LIV parameter determining if +A or -A is being tested; A occurs in the modified dispersion relation. LIV_A_sign has to be either +1 or -1 \n\
(--liv) this flag is now needed for launching a LIV run\n\
(--grtest-parameters dchi0,..,dxi1,..,dalpha1,..) template will assume deformations in the corresponding phase coefficients.\n\
(--ppe-parameters aPPE1,.... template will assume the presence of an arbitrary number of PPE parameters. They must be paired correctly.\n\
(--modeList lm,l-m...,lm,l-m) List of modes to be used by the model. The chosen modes ('lm') should be passed as a ',' seperated list.\n\
......@@ -1373,6 +1376,16 @@ LALInferenceModel *LALInferenceInitCBCModel(LALInferenceRunState *state) {
{
LALInferenceInitNonGRParams(state, model);
}
if (LALInferenceGetProcParamVal(commandLine,"--grtest-parameters"))
{
ppt=LALInferenceGetProcParamVal(commandLine,"--grtest-parameters");
if ((checkParamInList(ppt->value,"log10lambda_eff")) || (checkParamInList(ppt->value,"lambda_eff")) ){
if ( (approx != IMRPhenomPv2) && (approx != IMRPhenomD) && (approx != SEOBNRv4_ROM) && (approx != IMRPhenomPv2_NRTidal) && (approx != IMRPhenomPv2_NRTidalv2) && (approx != IMRPhenomPv3HM) && (approx != IMRPhenomHM) )
{
XLALPrintWarning("LIV parameters not compatible with approximant %s. Can be used only with IMRPhenomPv2, IMRPhenomD, SEOBNRv4_ROM, IMRPhenomPv2_NRTidalv2, IMRPhenomPv3HM, IMRPhenomHM, and IMRPhenomPv2_NRTidal.\n",XLALSimInspiralGetStringFromApproximant(approx));
}
}
}
/* PPE parameters */
ppt=LALInferenceGetProcParamVal(commandLine, "--TaylorF2ppE");
......@@ -1458,6 +1471,8 @@ LALInferenceModel *LALInferenceInitCBCModel(LALInferenceRunState *state) {
XLALSimInspiralWaveformParamsInsertNumRelData(model->LALpars, ppt->value);
fprintf(stdout,"Template will use %s.\n",ppt->value);
}
if (LALInferenceGetProcParamVal(commandLine, "--liv"))
XLALSimInspiralWaveformParamsInsertEnableLIV(model->LALpars, 1);
/* Pass custom mode array list to waveform generator */
if((ppt=LALInferenceGetProcParamVal(commandLine,"--modeList"))) {
......@@ -2282,6 +2297,8 @@ static void LALInferenceInitNonGRParams(LALInferenceRunState *state, LALInferenc
{
ProcessParamsTable *commandLine = state->commandLine;
ProcessParamsTable *ppt=NULL;
ProcessParamsTable *ppta=NULL;
ProcessParamsTable *pptb=NULL;
/* check that the user does not request both a TaylorF2Test and a PPE waveform model */
if (LALInferenceGetProcParamVal(commandLine,"--grtest-parameters") && LALInferenceGetProcParamVal(commandLine,"--ppe-parameters"))
{
......@@ -2302,6 +2319,24 @@ static void LALInferenceInitNonGRParams(LALInferenceRunState *state, LALInferenc
REAL8 dsigma_max=1.;
REAL8 dsigma_min=-1.;
REAL8 tmpVal=0.0;
if ((pptb=LALInferenceGetProcParamVal(commandLine,"--LIV_A_sign"))) {
REAL8 LIV_A_sign;
LIV_A_sign = atof(pptb->value);
if ((LIV_A_sign != -1) && (LIV_A_sign != 1)) {
XLALPrintError("LIV_A_sign can only take either +1 or -1.\n");
XLAL_ERROR_VOID(XLAL_EFAULT);
}
else LALInferenceAddVariable(model->params,"LIV_A_sign", &LIV_A_sign, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_FIXED);
}
if ((ppta=LALInferenceGetProcParamVal(commandLine,"--nonGR_alpha"))) {
REAL8 nonGR_alpha;
nonGR_alpha = atof(ppta->value);
if (nonGR_alpha==2.0) {
XLALPrintError("nonGR_alpha=2 is degenerate with p^2c^2 term in dispersion relation and cannot be tested. Please choose a different value.\n");
XLAL_ERROR_VOID(XLAL_EFAULT);
}
else LALInferenceAddVariable(model->params,"nonGR_alpha", &nonGR_alpha, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_FIXED);
}
/* Relative shifts for inspiral phase PN coefficients (absolute value for dchi1) */
if (checkParamInList(ppt->value,"dchi0")) LALInferenceRegisterUniformVariableREAL8(state, model->params, "dchi0", tmpVal, dchi_min, dchi_max, LALINFERENCE_PARAM_LINEAR);
if (checkParamInList(ppt->value,"dchi1")) LALInferenceRegisterUniformVariableREAL8(state, model->params, "dchi1", tmpVal, dchi_min, dchi_max, LALINFERENCE_PARAM_LINEAR);
......@@ -2335,6 +2370,36 @@ static void LALInferenceInitNonGRParams(LALInferenceRunState *state, LALInferenc
if (checkParamInList(ppt->value,"dbeta1")) LALInferenceRegisterUniformVariableREAL8(state, model->params, "dbeta1", tmpVal, dbeta_min, dbeta_max, LALINFERENCE_PARAM_LINEAR);
if (checkParamInList(ppt->value,"dbeta2")) LALInferenceRegisterUniformVariableREAL8(state, model->params, "dbeta2", tmpVal, dbeta_min, dbeta_max, LALINFERENCE_PARAM_LINEAR);
if (checkParamInList(ppt->value,"dbeta3")) LALInferenceRegisterUniformVariableREAL8(state, model->params, "dbeta3", tmpVal, dbeta_min, dbeta_max, LALINFERENCE_PARAM_LINEAR);
if (checkParamInList(ppt->value,"lambda_eff")) {
if (ppta==NULL) {
XLALPrintError("A value for nonGR_alpha has to be passed with lambda_eff.\n");
XLAL_ERROR_VOID(XLAL_EFAULT);
}
else if (pptb==NULL) {
XLALPrintError("A value of +1 or -1 for LIV_A_sign has to be passed with lambda_eff, respectively determing if +A or -A in the MDR is being tested.\n");
XLAL_ERROR_VOID(XLAL_EFAULT);
}
else {
REAL8 lambda_eff_min = 1E-6;
REAL8 lambda_eff_max = 1E13;
LALInferenceRegisterUniformVariableREAL8(state, model->params, "lambda_eff", 1E11, lambda_eff_min, lambda_eff_max, LALINFERENCE_PARAM_LINEAR);
}
}
if (checkParamInList(ppt->value,"log10lambda_eff")) {
if (ppta==NULL) {
XLALPrintError("A value for nonGR_alpha has to be passed with log10lambda_eff. A value for LIV_A_sign also has to be passed!\n");
XLAL_ERROR_VOID(XLAL_EFAULT);
}
else if (pptb==NULL) {
XLALPrintError("A value of +1 or -1 for LIV_A_sign has to be passed with log10lambda_eff, respectively determing if +A or -A in the MDR is being tested.\n");
XLAL_ERROR_VOID(XLAL_EFAULT);
}
else {
REAL8 log10lambda_eff_min = -6.0;
REAL8 log10lambda_eff_max = 13.0;
LALInferenceRegisterUniformVariableREAL8(state, model->params, "log10lambda_eff", 3.0, log10lambda_eff_min, log10lambda_eff_max, LALINFERENCE_PARAM_LINEAR);
}
}
}
ppt=LALInferenceGetProcParamVal(commandLine,"--ppe-parameters");
if (ppt)
......
......@@ -96,7 +96,7 @@ const char *const splineCalibrationProposalName = "SplineCalibration";
const char *const distanceLikelihoodProposalName = "DistanceLikelihood";
static const char *intrinsicNames[] = {"chirpmass", "q", "eta", "mass1", "mass2", "a_spin1", "a_spin2",
"tilt_spin1", "tilt_spin2", "phi12", "phi_jl", "frequency", "quality", "duration","polar_angle", "phase", "polar_eccentricity","dchi0","dchi1","dchi2","dchi3","dchi4","dchi5","dchi5l","dchi6","dchi6l","dchi7","aPPE","alphaPPE","bPPE","betaPPE","betaStep","fStep","dxi1","dxi2","dxi3","dxi4","dxi5","dxi6","dalpha1","dalpha2","dalpha3","dalpha4","dalpha5","dbeta1","dbeta2","dbeta3","dsigma1","dsigma2","dsigma3","dsigma4","lambda1","lambda2","lambdaT","dlambdaT","logp1", "gamma1", "gamma2", "gamma3", "SDgamma0","SDgamma1","SDgamma2","SDgamma3",NULL};
"tilt_spin1", "tilt_spin2", "phi12", "phi_jl", "frequency", "quality", "duration","polar_angle", "phase", "polar_eccentricity","dchi0","dchi1","dchi2","dchi3","dchi4","dchi5","dchi5l","dchi6","dchi6l","dchi7","aPPE","alphaPPE","bPPE","betaPPE","betaStep","fStep","dxi1","dxi2","dxi3","dxi4","dxi5","dxi6","dalpha1","dalpha2","dalpha3","dalpha4","dalpha5","dbeta1","dbeta2","dbeta3","dsigma1","dsigma2","dsigma3","dsigma4","lambda1","lambda2","lambdaT","dlambdaT","logp1", "gamma1", "gamma2", "gamma3", "SDgamma0","SDgamma1","SDgamma2","SDgamma3","log10lambda_eff","lambda_eff","nonGR_alpha","LIV_A_sign",NULL};
static const char *extrinsicNames[] = {"rightascension", "declination", "cosalpha", "azimuth", "polarisation", "distance",
......
......@@ -70,10 +70,11 @@
static void q2masses(double mc, double q, double *m1, double *m2);
/* list of testing GR parameters to be passed to the waveform */
/* the first batch of parameters dchis through dsigmas refer to the parameterised tests for generation (TIGER) while the parameters log10lambda_eff through LIV_A_sign are testing coefficients for the parameterised tests for propagation using a deformed dispersion relation (LIV); new parameters may be added at the end for readability although the order of parameters in this list does not matter */
const char list_extra_parameters[34][16] = {"dchi0","dchi1","dchi2","dchi3","dchi4","dchi5","dchi5l","dchi6","dchi6l","dchi7","aPPE","alphaPPE","bPPE","betaPPE","betaStep","fStep","dxi1","dxi2","dxi3","dxi4","dxi5","dxi6","dalpha1","dalpha2","dalpha3","dalpha4","dalpha5","dbeta1","dbeta2","dbeta3","dsigma1","dsigma2","dsigma3","dsigma4"};
const char list_extra_parameters[38][16] = {"dchi0","dchi1","dchi2","dchi3","dchi4","dchi5","dchi5l","dchi6","dchi6l","dchi7","aPPE","alphaPPE","bPPE","betaPPE","betaStep","fStep","dxi1","dxi2","dxi3","dxi4","dxi5","dxi6","dalpha1","dalpha2","dalpha3","dalpha4","dalpha5","dbeta1","dbeta2","dbeta3","dsigma1","dsigma2","dsigma3","dsigma4","log10lambda_eff","lambda_eff","nonGR_alpha","LIV_A_sign"};
const UINT4 N_extra_params = 34;
const UINT4 N_extra_params = 38;
/* Return the quadrupole moment of a neutron star given its lambda
* We use the relations defined here. https://arxiv.org/pdf/1302.4499.pdf.
......
......@@ -56,6 +56,7 @@ from scipy import special
from scipy import signal
from scipy.optimize import newton
from scipy import interpolate
from scipy import integrate
from numpy import linspace
import random
import socket
......@@ -141,11 +142,12 @@ ppEParams=['ppEalpha','ppElowera','ppEupperA','ppEbeta','ppElowerb','ppEupperB',
tigerParams=['dchi%i'%(i) for i in range(8)] + ['dchi%il'%(i) for i in [5,6] ] + ['dxi%d'%(i+1) for i in range(6)] + ['dalpha%i'%(i+1) for i in range(5)] + ['dbeta%i'%(i+1) for i in range(3)] + ['dsigma%i'%(i+1) for i in range(4)]
bransDickeParams=['omegaBD','ScalarCharge1','ScalarCharge2']
massiveGravitonParams=['lambdaG']
lorentzInvarianceViolationParams=['log10lambda_a','lambda_a','log10lambda_eff','lambda_eff','log10livamp','liv_amp']
tidalParams=['lambda1','lambda2','lam_tilde','dlam_tilde','lambdat','dlambdat','lambdas','bluni']
fourPiecePolyParams=['logp1','gamma1','gamma2','gamma3']
spectralParams=['sdgamma0','sdgamma1','sdgamma2','sdgamma3']
energyParams=['e_rad', 'e_rad_evol', 'e_rad_nonevol', 'l_peak', 'l_peak_evol', 'l_peak_nonevol', 'e_rad_maxldist', 'e_rad_maxldist_evol', 'e_rad_maxldist_nonevol']
strongFieldParams=ppEParams+tigerParams+bransDickeParams+massiveGravitonParams+tidalParams+fourPiecePolyParams+spectralParams+energyParams
strongFieldParams=ppEParams+tigerParams+bransDickeParams+massiveGravitonParams+tidalParams+fourPiecePolyParams+spectralParams+energyParams+lorentzInvarianceViolationParams
#Extrinsic
distParams=['distance','distMPC','dist','distance_maxl']
......@@ -170,7 +172,7 @@ for derived_time in ['h1_end_time','l1_end_time','v1_end_time','h1l1_delay','l1v
greedyBinSizes[derived_time]=greedyBinSizes['time']
for derived_phase in relativePhaseParams:
greedyBinSizes[derived_phase]=0.05
for param in tigerParams + bransDickeParams + massiveGravitonParams:
for param in tigerParams + bransDickeParams + massiveGravitonParams + lorentzInvarianceViolationParams:
greedyBinSizes[param]=0.01
for param in tidalParams:
greedyBinSizes[param]=2.5
......@@ -555,7 +557,13 @@ def plot_label(param):
'v1_optimal_snr':r'$\rho^{opt}_{V1}$',
'matched_filter_snr':r'$\rho^{MF}$',
'lambdas':r'$\Lambda_S$',
'bluni' : r'$BL_{uniform}$'
'bluni' : r'$BL_{uniform}$',
'log10lambda_a':r'$\log\lambda_{\mathbb{A}} [\mathrm{m}]$',
'log10lambda_eff':r'$\log\lambda_{eff} [\mathrm{m}]$',
'lambda_eff':r'$\lambda_{eff} [\mathrm{m}]$',
'lambda_a':r'$\lambda_{\mathbb{A}} [\mathrm{m}]$',
'liv_amp':r'$\mathbb{A} [\mathrm{{eV}^{2-\alpha}}]$' ,
'log10livamp':r'$\log \mathbb{A}[\mathrm{{eV}^{2-\alpha}}]$'
}
# Handle cases where multiple names have been used
......@@ -1137,6 +1145,16 @@ class Posterior(object):
if ('mc' in pos.names) and ('redshift_maxldist' in pos.names):
pos.append_mapping('mc_source_maxldist', source_mass, ['mc', 'redshift_maxldist'])
# Calling functions testing Lorentz invariance violation
if ('log10lambda_eff' in pos.names) and ('redshift' in pos.names):
pos.append_mapping('log10lambda_a', lambda z,nonGR_alpha,wl,dist:np.log10(lambda_a(z, nonGR_alpha, 10**wl, dist)), ['redshift', 'nonGR_alpha', 'log10lambda_eff', 'dist'])
if ('log10lambda_eff' in pos.names) and ('redshift' in pos.names):
pos.append_mapping('log10livamp', lambda z,nonGR_alpha,wl,dist:np.log10(amplitudeMeasure(z, nonGR_alpha, 10**wl, dist)), ['redshift','nonGR_alpha','log10lambda_eff', 'dist'])
if ('lambda_eff' in pos.names) and ('redshift' in pos.names):
pos.append_mapping('lambda_a', lambda_a, ['redshift', 'nonGR_alpha', 'log10lambda_eff', 'dist'])
if ('lambda_eff' in pos.names) and ('redshift' in pos.names):
pos.append_mapping('liv_amp', amplitudeMeasure, ['redshift','nonGR_alpha','lambda_eff', 'dist'])
#Calculate new tidal parameters
new_tidal_params = ['lam_tilde','dlam_tilde']
old_tidal_params = ['lambda1','lambda2','q']
......@@ -3942,6 +3960,52 @@ def source_mass(mass, redshift):
"""
return mass / (1.0 + redshift)
## Following functions added for testing Lorentz violations
def integrand_distance(redshift,nonGR_alpha):
"""
Calculate D_alpha integral; multiplicative factor put later
D_alpha = integral{ ((1+z')^(alpha-2))/sqrt(Omega_m*(1+z')^3 +Omega_lambda) dz'} # eq.15 of arxiv 1110.2720
"""
omega = lal.CreateCosmologicalParameters(0.6790,0.3065,0.6935,-1.0,0.0,0.0) ## Planck 2015 values
omega_m = omega.om # matter density
omega_l = omega.ol # dark energy density
#lal.DestroyCosmologicalParameters(omega)
return (1.0+redshift)**(nonGR_alpha-2.0)/(np.sqrt(omega_m*(1.0+redshift)**3.0 + omega_l))
def DistanceMeasure(redshift,nonGR_alpha):
"""
D_alpha = ((1+z)^(1-alpha))/H_0 * D_alpha # from eq.15 of arxiv 1110.2720
D_alpha calculated from integrand in above function
"""
omega = lal.CreateCosmologicalParameters(0.6790,0.3065,0.6935,-1.0,0.0,0.0) ## Planck 2015 values
H0 = omega.h*lal.H0FAC_SI ## Hubble constant in SI units
dist = integrate.quad(integrand_distance, 0, redshift ,args=(nonGR_alpha))[0]
dist *= (1.0 + redshift)**(1.0 - nonGR_alpha)
dist /= H0
#lal.DestroyCosmologicalParameters(omega)
return dist*lal.C_SI ## returns D_alpha in metres
def lambda_a(redshift, nonGR_alpha, lambda_eff, distance):
"""
Converting from the effective wavelength-like parameter to lambda_A:
lambda_A = lambda_{eff}*(D_alpha/D_L)^(1/(2-alpha))*(1/(1+z)^((1-alpha)/(2-alpha)))
"""
Dfunc = np.vectorize(DistanceMeasure)
D_alpha = Dfunc(redshift, nonGR_alpha)
dl = distance*lal.PC_SI*1e6 ## luminosity distane in metres
return lambda_eff*(D_alpha/(dl*(1.0+redshift)**(1.0-nonGR_alpha)))**(1./(2.0-nonGR_alpha))
def amplitudeMeasure(redshift, nonGR_alpha, lambda_eff, distance):
"""
Converting to Lorentz violating parameter "A" in dispersion relation from lambda_A:
A = (lambda_A/h)^(alpha-2) # eqn. 13 of arxiv 1110.2720
"""
hPlanck = 4.13567e-15 # Planck's constant in eV.s
ampFunc = np.vectorize(lambda_a)
lambdaA = ampFunc(redshift, nonGR_alpha, lambda_eff, distance)/lal.C_SI # convert to seconds
return (lambdaA/hPlanck)**(nonGR_alpha-2.0)
############################ changes for testing Lorentz violations made till here
def physical2radiationFrame(theta_jn, phi_jl, tilt1, tilt2, phi12, a1, a2, m1, m2, fref,phiref):
"""
Wrapper function for SimInspiralTransformPrecessingNewInitialConditions().
......@@ -6080,7 +6144,7 @@ class PEOutputParser(object):
llines.append(np.array(list(map(lambda a:float(a.text),row))))
flines=np.array(llines)
for i in range(0,len(header)):
if header[i].lower().find('log')!=-1 and header[i].lower() not in logParams and re.sub('log', '', header[i].lower()) not in [h.lower() for h in header]:
if header[i].lower().find('log')!=-1 and header[i].lower() not in logParams and re.sub('log', '', header[i].lower()) not in [h.lower() for h in header] and header[i].lower() not in lorentzInvarianceViolationParams:
print('exponentiating %s'%(header[i]))
flines[:,i]=np.exp(flines[:,i])
......@@ -6152,7 +6216,7 @@ class PEOutputParser(object):
for param in params:
param_low = param.lower()
if param_low.find('log') != -1 and param_low not in logParams and re.sub('log', '', param_low) not in [p.lower() for p in params]:
if param_low.find('log') != -1 and param_low not in logParams and re.sub('log', '', param_low) not in [p.lower() for p in params] and param_low not in lorentzInvarianceViolationParams:
print('exponentiating %s' % param)
new_param = re.sub('log', '', param, flags=re.IGNORECASE)
samples[new_param] = np.exp(samples[param])
......@@ -6280,7 +6344,7 @@ class PEOutputParser(object):
for i in range(0,len(header)):
if header[i].lower().find('log')!=-1 and header[i].lower() not in logParams and re.sub('log', '', header[i].lower()) not in [h.lower() for h in header]:
if header[i].lower().find('log')!=-1 and header[i].lower() not in logParams and re.sub('log', '', header[i].lower()) not in [h.lower() for h in header] and header[i].lower() not in lorentzInvarianceViolationParams:
print('exponentiating %s'%(header[i]))
flines[:,i]=np.exp(flines[:,i])
......
......@@ -2271,6 +2271,9 @@ int XLALSimInspiralChooseFDWaveform(
}
if (ret == XLAL_FAILURE) XLAL_ERROR(XLAL_EFUNC);
if (XLALSimInspiralWaveformParamsLookupEnableLIV(LALparams))
ret = XLALSimLorentzInvarianceViolationTerm(hptilde, hctilde, m1/LAL_MSUN_SI, m2/LAL_MSUN_SI, distance, LALparams);
if (ret == XLAL_FAILURE) XLAL_ERROR(XLAL_EFUNC);
return ret;
}
......@@ -6167,7 +6170,6 @@ int XLALSimInspiralApproximantAcceptTestGRParams(Approximant approx){
case SEOBNRv2_ROM_DoubleSpin:
case SEOBNRv2_ROM_DoubleSpin_HI:
case Lackey_Tidal_2013_SEOBNRv2_ROM:
case SEOBNRv4_ROM:
case SEOBNRv4HM_ROM:
case SEOBNRv4_ROM_NRTidal:
case SEOBNRv4_ROM_NRTidalv2:
......@@ -6206,6 +6208,7 @@ int XLALSimInspiralApproximantAcceptTestGRParams(Approximant approx){
case PhenSpinTaylor:
case PhenSpinTaylorRD:
case EccentricTD:
case SEOBNRv4_ROM:
case IMRPhenomC:
case IMRPhenomD:
case IMRPhenomP:
......@@ -6226,6 +6229,79 @@ int XLALSimInspiralApproximantAcceptTestGRParams(Approximant approx){
return testGR_accept;
};
/* Function for introducing Lorentz violating changes in FD phase; calculates eqns. 30 & 32 of arxiv 1110.2720 for the LV phase term in FD and multiplies to h+ and hx */
int XLALSimLorentzInvarianceViolationTerm(
COMPLEX16FrequencySeries **hptilde, /**< Frequency-domain waveform h+ */
COMPLEX16FrequencySeries **hctilde, /**< Frequency-domain waveform hx */
REAL8 m1, /**< Mass 1 in solar masses */
REAL8 m2, /**< Mass 2 in solar masses */
REAL8 r, /**< distance in metres*/
LALDict *LALparams /**< LAL dictionary containing accessory parameters */
)
{
REAL8 f0, f, df;
COMPLEX16 hplus, hcross, tmpExp;
REAL8 M, eta, zeta, dPhiPref, Mc, tmpVal;
UINT4 len, i;
M = m1+m2;
eta = m1*m2/(M*M);
Mc = M*pow(eta, 0.6);
len = (*hptilde)->data->length;
REAL8 lambda_eff = pow(10,XLALSimInspiralWaveformParamsLookupNonGRLIVLogLambdaEff(LALparams)); /* Effective wavelength-like parameter in phase in metres */
REAL8 nonGR_alpha = XLALSimInspiralWaveformParamsLookupNonGRLIVAlpha(LALparams); /* Exponent defined in terms of PN order characterising LIV*/
REAL8 LIV_A_sign = XLALSimInspiralWaveformParamsLookupNonGRLIVASign(LALparams); /* Sign of A determining the sign of LV phase */
if ((*hctilde)->data->length != len) {
XLALPrintError("Lengths of plus and cross polarization series do not agree \n");
XLAL_ERROR(XLAL_EBADLEN);
}
f0 = (*hptilde)->f0;
if ((*hctilde)->f0 != f0) {
XLALPrintError("Starting frequencies of plus and cross polarization series do not agree \n");
XLAL_ERROR(XLAL_EINVAL);
}
df = (*hptilde)->deltaF;
if ((*hctilde)->deltaF != df) {
XLALPrintError("Frequency steps of plus and cross polarization series do not agree \n");
XLAL_ERROR(XLAL_EINVAL);
}
UINT4 k = 0;
if (f0 == 0.0)
k=1;
if (nonGR_alpha == 1) {
zeta = LIV_A_sign*LAL_PI*r/lambda_eff; /*Eqn. (32) of arxiv:1110.2720*/
dPhiPref = zeta*log(LAL_PI*Mc*LAL_MTSUN_SI); /*Eqn. (31) of arxiv:1110.2720;the frequency dependence is treated below*/
for (i=k; i<len; i++) {
f = f0 + i*df;
tmpExp = cexp(I*(dPhiPref + zeta*log(f)));
hplus = (*hptilde)->data->data[i] * tmpExp;
(*hptilde)->data->data[i] = hplus;
hcross = (*hctilde)->data->data[i] * tmpExp;
(*hctilde)->data->data[i] = hcross;
}
}
else {
zeta = LIV_A_sign*pow(LAL_PI, (2. - nonGR_alpha))*r*pow(Mc*LAL_MRSUN_SI, (1. - nonGR_alpha))/((1. - nonGR_alpha)*pow(lambda_eff, (2. - nonGR_alpha))); /*Eqn. (30) of arxiv:1110.2720*/
dPhiPref = zeta*pow(LAL_PI*Mc*LAL_MTSUN_SI, (nonGR_alpha - 1.)); /*Eqn. (28) of arxiv:1110.2720;the frequency dependence is treated below*/
for (i=k; i<len; i++) {
f = f0 + i*df;
tmpVal = pow(f, (nonGR_alpha - 1.));
tmpExp=cexp(-I*dPhiPref*tmpVal);
hplus = (*hptilde)->data->data[i] * tmpExp;
(*hptilde)->data->data[i] = hplus;
hcross = (*hctilde)->data->data[i] * tmpExp;
(*hctilde)->data->data[i] = hcross;
}
}
return XLAL_SUCCESS;
}
/** @} */
/**
......
......@@ -918,6 +918,8 @@ SphHarmTimeSeries *XLALSimIMRNRHybSur3dq8Modes(
LALDict* LALparams /**< Dict with extra parameters */
);
/* routine for checking Lorentz violation */
int XLALSimLorentzInvarianceViolationTerm(COMPLEX16FrequencySeries **hptilde, COMPLEX16FrequencySeries **hctilde, REAL8 m1, REAL8 m2, REAL8 r, LALDict *LALparams);
#if 0
{ /* so that editors will match succeeding brace */
......
......@@ -135,6 +135,10 @@ DEFINE_INSERT_FUNC(NonGRAlphaPPE6, REAL8, "alphaPPE6", 0)
DEFINE_INSERT_FUNC(NonGRBetaPPE6, REAL8, "betaPPE6", 0)
DEFINE_INSERT_FUNC(NonGRAlphaPPE7, REAL8, "alphaPPE7", 0)
DEFINE_INSERT_FUNC(NonGRBetaPPE7, REAL8, "betaPPE7", 0)
DEFINE_INSERT_FUNC(EnableLIV, INT4, "liv", 0)
DEFINE_INSERT_FUNC(NonGRLIVLogLambdaEff, REAL8, "log10lambda_eff", 100)
DEFINE_INSERT_FUNC(NonGRLIVASign, REAL8, "LIV_A_sign", 1)
DEFINE_INSERT_FUNC(NonGRLIVAlpha, REAL8, "nonGR_alpha", 0)
/* NLTides parameters */
/* used within LALSimInspiralTaylorF2NLTides.c */
......@@ -273,6 +277,10 @@ DEFINE_LOOKUP_FUNC(NonGRAlphaPPE6, REAL8, "alphaPPE6", 0)
DEFINE_LOOKUP_FUNC(NonGRBetaPPE6, REAL8, "betaPPE6", 0)
DEFINE_LOOKUP_FUNC(NonGRAlphaPPE7, REAL8, "alphaPPE7", 0)
DEFINE_LOOKUP_FUNC(NonGRBetaPPE7, REAL8, "betaPPE7", 0)
DEFINE_LOOKUP_FUNC(EnableLIV, INT4, "liv", 0)
DEFINE_LOOKUP_FUNC(NonGRLIVLogLambdaEff, REAL8, "log10lambda_eff", 100)
DEFINE_LOOKUP_FUNC(NonGRLIVASign, REAL8, "LIV_A_sign", 1)
DEFINE_LOOKUP_FUNC(NonGRLIVAlpha, REAL8, "nonGR_alpha", 0)
/* NLTides parameters */
/* used within LALSimInspiralTaylorF2NLTides.c */
......@@ -402,6 +410,10 @@ DEFINE_ISDEFAULT_FUNC(NonGRAlphaPPE6, REAL8, "alphaPPE6", 0)
DEFINE_ISDEFAULT_FUNC(NonGRBetaPPE6, REAL8, "betaPPE6", 0)
DEFINE_ISDEFAULT_FUNC(NonGRAlphaPPE7, REAL8, "alphaPPE7", 0)
DEFINE_ISDEFAULT_FUNC(NonGRBetaPPE7, REAL8, "betaPPE7", 0)
DEFINE_ISDEFAULT_FUNC(EnableLIV, INT4, "liv", 0)
DEFINE_ISDEFAULT_FUNC(NonGRLIVLogLambdaEff, REAL8, "log10lambda_eff", 100)
DEFINE_ISDEFAULT_FUNC(NonGRLIVASign, REAL8, "LIV_A_sign", 1)
DEFINE_ISDEFAULT_FUNC(NonGRLIVAlpha, REAL8, "nonGR_alpha", 0)
/* SEOBNRv4P */
DEFINE_ISDEFAULT_FUNC(EOBChooseNumOrAnalHamDer, INT4, "EOBChooseNumOrAnalHamDer", 1)
......
......@@ -122,6 +122,10 @@ int XLALSimInspiralWaveformParamsInsertNonGRAlphaPPE6(LALDict *params, REAL8 val
int XLALSimInspiralWaveformParamsInsertNonGRBetaPPE6(LALDict *params, REAL8 value);
int XLALSimInspiralWaveformParamsInsertNonGRAlphaPPE7(LALDict *params, REAL8 value);
int XLALSimInspiralWaveformParamsInsertNonGRBetaPPE7(LALDict *params, REAL8 value);
int XLALSimInspiralWaveformParamsInsertEnableLIV(LALDict *params, INT4 value);
int XLALSimInspiralWaveformParamsInsertNonGRLIVLogLambdaEff(LALDict *params, REAL8 value);
int XLALSimInspiralWaveformParamsInsertNonGRLIVASign(LALDict *params, REAL8 value);
int XLALSimInspiralWaveformParamsInsertNonGRLIVAlpha(LALDict *params, REAL8 value);
/* NLTides parameters */
/* used within LALSimInspiralTaylorF2NLTides.c */
......@@ -246,6 +250,10 @@ REAL8 XLALSimInspiralWaveformParamsLookupNonGRAlphaPPE6(LALDict *params);
REAL8 XLALSimInspiralWaveformParamsLookupNonGRBetaPPE6(LALDict *params);
REAL8 XLALSimInspiralWaveformParamsLookupNonGRAlphaPPE7(LALDict *params);
REAL8 XLALSimInspiralWaveformParamsLookupNonGRBetaPPE7(LALDict *params);
INT4 XLALSimInspiralWaveformParamsLookupEnableLIV(LALDict *params);
REAL8 XLALSimInspiralWaveformParamsLookupNonGRLIVLogLambdaEff(LALDict *params);
REAL8 XLALSimInspiralWaveformParamsLookupNonGRLIVASign(LALDict *params);
REAL8 XLALSimInspiralWaveformParamsLookupNonGRLIVAlpha(LALDict *params);
/* NLTides parameters */
/* used within LALSimInspiralTaylorF2NLTides.c */
......@@ -369,6 +377,10 @@ int XLALSimInspiralWaveformParamsNonGRAlphaPPE6IsDefault(LALDict *params);
int XLALSimInspiralWaveformParamsNonGRBetaPPE6IsDefault(LALDict *params);
int XLALSimInspiralWaveformParamsNonGRAlphaPPE7IsDefault(LALDict *params);
int XLALSimInspiralWaveformParamsNonGRBetaPPE7IsDefault(LALDict *params);
int XLALSimInspiralWaveformParamsEnableLIVIsDefault(LALDict *params);
int XLALSimInspiralWaveformParamsNonGRLIVLogLambdaEffIsDefault(LALDict *params);
int XLALSimInspiralWaveformParamsNonGRLIVASignIsDefault(LALDict *params);
int XLALSimInspiralWaveformParamsNonGRLIVAlphaIsDefault(LALDict *params);
/* SEOBNRv4P */
INT4 XLALSimInspiralWaveformParamsEOBChooseNumOrAnalHamDerIsDefault(LALDict *params);
......
......@@ -14,6 +14,7 @@ test_scripts += \
test_phenomX.py \
test_SEOBNRv4HM_ROM.py \
test_enum_compatibility.py \
test_liv.py \
$(END_OF_LIST)
......
#!/usr/bin/env python
# -*- coding: utf-8 -*-
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http: //www.gnu.org/licenses/>.
# define required parameters for time domain and frequency domain
# as well as their default values
"""Simple test to see if the LIV parameters can be read and inserted correctly
"""
import sys
import pytest
import lal
import lalsimulation
import numpy as np
# -- utility functions ---------------------
def read_liv_params(LALparams):
"""
Reads LIV parameters
"""
logLeff = lalsimulation.SimInspiralWaveformParamsLookupNonGRLIVLogLambdaEff(LALparams)
signOfA = lalsimulation.SimInspiralWaveformParamsLookupNonGRLIVASign(LALparams)
alpha = lalsimulation.SimInspiralWaveformParamsLookupNonGRLIVAlpha(LALparams)
return logLeff, signOfA, alpha