1 from __future__ import division
2 import os,sys,math,copy
3 import numpy
4 import matplotlib
5 matplotlib.use('Agg')
6 import pylab
7 from lal import PI as LAL_PI
8 from lal import MTSUN_SI as LAL_MTSUN_SI
9
10
11 mtsun = LAL_MTSUN_SI
12 LAL_GAMMA = 0.5772156649015328606065120900824024
13
14
15
17 def is_exe(fpath):
18 return os.path.isfile(fpath) and os.access(fpath, os.X_OK)
19
20 fpath, fname = os.path.split(program)
21 if fpath:
22 if is_exe(program):
23 return program
24 else:
25 for path in os.environ["PATH"].split(os.pathsep):
26 exe_file = os.path.join(path, program)
27 if is_exe(exe_file):
28 return exe_file
29
30 return None
31
32 -def determine_eigen_directions(psd,order,f0,f_low,f_upper,delta_f,\
33 return_gs=False,verbose=False,elapsed_time=None,\
34 vary_fmax=False,vary_density=25):
35
36 evals = {}
37 evecs = {}
38 metric = {}
39 if verbose:
40 print >>sys.stdout,"Beginning to calculate moments at %d." %(elapsed_time())
41
42 moments = get_moments(psd,f0,f_low,f_upper,delta_f,vary_fmax=vary_fmax,\
43 vary_density=vary_density)
44
45 if verbose:
46 print >>sys.stdout,"Moments calculated, transforming to metric at %d." \
47 %(elapsed_time())
48
49 list = []
50 if vary_fmax:
51 for t_fmax in numpy.arange(f_low+vary_density,f_upper,vary_density):
52 list.append(t_fmax)
53 else:
54 list.append('fixed')
55
56 for item in list:
57 Js = []
58 for i in range(18):
59 Js.append(moments['J%d'%(i)][item])
60 Js.append(moments['J%d'%(-1)][item])
61
62 logJs = []
63 for i in range(18):
64 logJs.append(moments['log%d'%(i)][item])
65 logJs.append(moments['log%d'%(-1)][item])
66
67 loglogJs = []
68 for i in range(18):
69 loglogJs.append(moments['loglog%d'%(i)][item])
70 loglogJs.append(moments['loglog%d'%(-1)][item])
71
72 logloglogJs = []
73 for i in range(18):
74 logloglogJs.append(moments['logloglog%d'%(i)][item])
75 logloglogJs.append(moments['logloglog%d'%(-1)][item])
76
77 loglogloglogJs = []
78 for i in range(18):
79 loglogloglogJs.append(moments['loglogloglog%d'%(i)][item])
80 loglogloglogJs.append(moments['loglogloglog%d'%(-1)][item])
81
82 mapping = {}
83
84 if order == 'twoPN':
85 maxLen = 4
86 gs = numpy.matrix(numpy.zeros(shape=(maxLen,maxLen),dtype=float))
87 mapping['Lambda0'] = 0
88 mapping['Lambda2'] = 1
89 mapping['Lambda3'] = 2
90 mapping['Lambda4'] = 3
91 elif order == 'threePointFivePN':
92 maxLen = 8
93 gs = numpy.matrix(numpy.zeros(shape=(maxLen,maxLen),dtype=float))
94 mapping['Lambda0'] = 0
95 mapping['Lambda2'] = 1
96 mapping['Lambda3'] = 2
97 mapping['Lambda4'] = 3
98 mapping['LogLambda5'] = 4
99 mapping['Lambda6'] = 5
100 mapping['Lambda7'] = 6
101 mapping['LogLambda6'] = 7
102 elif order == 'taylorF4_45PN':
103 maxLen = 12
104 gs = numpy.matrix(numpy.zeros(shape=(maxLen,maxLen),dtype=float))
105 mapping['Lambda0'] = 0
106 mapping['Lambda2'] = 1
107 mapping['Lambda3'] = 2
108 mapping['Lambda4'] = 3
109 mapping['LogLambda5'] = 4
110 mapping['Lambda6'] = 5
111 mapping['Lambda7'] = 6
112 mapping['LogLambda6'] = 7
113 mapping['LogLambda8'] = 8
114 mapping['LogLogLambda8'] = 9
115 mapping['Lambda9'] = 10
116 mapping['LogLambda9'] = 11
117 else:
118 raise BrokenError
119
120 for i in range(16):
121 for j in range(16):
122
123 if mapping.has_key('Lambda%d'%i) and mapping.has_key('Lambda%d'%j):
124 gs[mapping['Lambda%d'%i],mapping['Lambda%d'%j]] = 0.5 * (Js[17-i-j] - Js[12-i]*Js[12-j] - (Js[9-i] - Js[4]*Js[12-i]) * (Js[9-j] - Js[4] * Js[12-j])/(Js[1] - Js[4]*Js[4]))
125
126 if mapping.has_key('Lambda%d'%i) and mapping.has_key('LogLambda%d'%j):
127 gammaij = logJs[17-i-j] - logJs[12-j] * Js[12-i]
128 gamma0i = (Js[9-i] - Js[4] * Js[12-i])
129 gamma0j = logJs[9-j] - logJs[12-j] * Js[4]
130 gs[mapping['Lambda%d'%i],mapping['LogLambda%d'%j]] = \
131 gs[mapping['LogLambda%d'%j],mapping['Lambda%d'%i]] = \
132 0.5 * (gammaij - gamma0i*gamma0j/(Js[1] - Js[4]*Js[4]))
133
134 if mapping.has_key('LogLambda%d'%i) and mapping.has_key('LogLambda%d'%j):
135 gammaij = loglogJs[17-i-j] - logJs[12-j] * logJs[12-i]
136 gamma0i = (logJs[9-i] - Js[4] * logJs[12-i])
137 gamma0j = logJs[9-j] - logJs[12-j] * Js[4]
138 gs[mapping['LogLambda%d'%i],mapping['LogLambda%d'%j]] = \
139 0.5 * (gammaij - gamma0i*gamma0j/(Js[1] - Js[4]*Js[4]))
140
141 if mapping.has_key('Lambda%d'%i) and mapping.has_key('LogLogLambda%d'%j):
142 gammaij = loglogJs[17-i-j] - loglogJs[12-j] * Js[12-i]
143 gamma0i = (Js[9-i] - Js[4] * Js[12-i])
144 gamma0j = loglogJs[9-j] - loglogJs[12-j] * Js[4]
145 gs[mapping['Lambda%d'%i],mapping['LogLogLambda%d'%j]] = \
146 gs[mapping['LogLogLambda%d'%j],mapping['Lambda%d'%i]] = \
147 0.5 * (gammaij - gamma0i*gamma0j/(Js[1] - Js[4]*Js[4]))
148
149 if mapping.has_key('LogLambda%d'%i) and mapping.has_key('LogLogLambda%d'%j):
150 gammaij = logloglogJs[17-i-j] - loglogJs[12-j] * logJs[12-i]
151 gamma0i = (logJs[9-i] - Js[4] * logJs[12-i])
152 gamma0j = loglogJs[9-j] - loglogJs[12-j] * Js[4]
153 gs[mapping['LogLambda%d'%i],mapping['LogLogLambda%d'%j]] = \
154 gs[mapping['LogLogLambda%d'%j],mapping['LogLambda%d'%i]] = \
155 0.5 * (gammaij - gamma0i*gamma0j/(Js[1] - Js[4]*Js[4]))
156
157 if mapping.has_key('LogLogLambda%d'%i) and mapping.has_key('LogLogLambda%d'%j):
158 gammaij = loglogloglogJs[17-i-j] - loglogJs[12-j] * loglogJs[12-i]
159 gamma0i = (loglogJs[9-i] - Js[4] * loglogJs[12-i])
160 gamma0j = loglogJs[9-j] - loglogJs[12-j] * Js[4]
161 gs[mapping['LogLogLambda%d'%i],mapping['LogLogLambda%d'%j]] = \
162 0.5 * (gammaij - gamma0i*gamma0j/(Js[1] - Js[4]*Js[4]))
163
164 evals[item],evecs[item] = numpy.linalg.eig(gs)
165 metric[item] = numpy.matrix(gs)
166
167 for i in range(len(evals[item])):
168 if evals[item][i] < 0:
169 print "WARNING: Negative eigenvalue %e. Setting as positive." %(evals[item][i])
170 evals[item][i] = -evals[item][i]
171 if evecs[item][i,i] < 0:
172
173
174 evecs[item][:,i] = - evecs[item][:,i]
175
176 if verbose:
177 print >>sys.stdout,"Metric and eigenvalues calculated at %d." \
178 %(elapsed_time())
179
180 if return_gs:
181 return evals,evecs,gs
182
183 return evals,evecs
184
185 -def get_moments(psd_file,f0,f_low,f_high,deltaF,vary_fmax=False,vary_density=25):
186 psd = numpy.loadtxt(psd_file)
187 psd_f = psd[:,0]
188 psd_amp = psd[:,1]
189 psd_amp = psd_amp * psd_amp
190 new_f,new_amp = interpolate_psd(psd_f,psd_amp,deltaF)
191
192
193 funct = lambda x: 1
194 I7 = calculate_moment(new_f,new_amp,f_low,f_high,f0,funct,vary_fmax=vary_fmax,vary_density=vary_density)
195
196
197 moments = {}
198 for i in range(-1,18):
199 funct = lambda x: x**((-i+7)/3.)
200 moments['J%d' %(i)] = calculate_moment(new_f,new_amp,f_low,f_high,f0,funct,norm=I7,vary_fmax=vary_fmax,vary_density=vary_density)
201
202
203 for i in range(-1,18):
204 funct = lambda x: (numpy.log(x**(1./3.))) * x**((-i+7)/3.)
205 moments['log%d' %(i)] = calculate_moment(new_f,new_amp,f_low,f_high,f0,funct,norm=I7,vary_fmax=vary_fmax,vary_density=vary_density)
206
207
208 for i in range(-1,18):
209 funct = lambda x: (numpy.log(x**(1./3.))) * (numpy.log(x**(1./3.))) * x**((-i+7)/3.)
210 moments['loglog%d' %(i)] = calculate_moment(new_f,new_amp,f_low,f_high,f0,funct,norm=I7,vary_fmax=vary_fmax,vary_density=vary_density)
211
212
213 for i in range(-1,18):
214 funct = lambda x: (numpy.log(x**(1./3.))) * (numpy.log(x**(1./3.))) * (numpy.log(x**(1./3.))) * x**((-i+7)/3.)
215 moments['logloglog%d' %(i)] = calculate_moment(new_f,new_amp,f_low,f_high,f0,funct,norm=I7,vary_fmax=vary_fmax,vary_density=vary_density)
216
217
218 for i in range(-1,18):
219 funct = lambda x: (numpy.log(x**(1./3.))) * (numpy.log(x**(1./3.))) * (numpy.log(x**(1./3.))) * (numpy.log(x**(1./3.))) * x**((-i+7)/3.)
220 moments['loglogloglog%d' %(i)] = calculate_moment(new_f,new_amp,f_low,f_high,f0,funct,norm=I7,vary_fmax=vary_fmax,vary_density=vary_density)
221
222 return moments
223
225 new_psd_f = []
226 new_psd_amp = []
227 fcurr = psd_f[0]
228
229 for i in range(len(psd_f) - 1):
230 f_low = psd_f[i]
231 f_high = psd_f[i+1]
232 amp_low = psd_amp[i]
233 amp_high = psd_amp[i+1]
234 while(1):
235 if fcurr > f_high:
236 break
237 new_psd_f.append(fcurr)
238 gradient = (amp_high - amp_low) / (f_high - f_low)
239 fDiff = fcurr - f_low
240 new_psd_amp.append(amp_low + fDiff * gradient)
241 fcurr = fcurr + deltaF
242 return numpy.asarray(new_psd_f),numpy.asarray(new_psd_amp)
243
244
245 -def calculate_moment(psd_f,psd_amp,fmin,fmax,f0,funct,norm=None,vary_fmax=False,vary_density=25):
246
247 psd_x = psd_f / f0
248 deltax = psd_x[1] - psd_x[0]
249
250 comps = (psd_x)**(-7./3.) * funct(psd_x) * deltax/ psd_amp
251 moment = {}
252 logica = numpy.logical_and(psd_f > fmin, psd_f < fmax)
253 comps_red = comps[logica]
254 psdf_red = psd_f[logica]
255 moment['fixed'] = comps_red.sum()
256 if norm:
257 moment['fixed'] = moment['fixed']/norm['fixed']
258 if vary_fmax:
259 for t_fmax in numpy.arange(fmin+vary_density,fmax,vary_density):
260 comps_red2 = comps_red[psdf_red < t_fmax]
261 moment[t_fmax] = comps_red2.sum()
262 if norm:
263 moment[t_fmax] = moment[t_fmax]/norm[t_fmax]
264 return moment
265
266 -def estimate_mass_range_slimline(numiters,order,evals,evecs,maxmass1,minmass1,maxmass2,minmass2,maxspin,f0,covary=True,maxBHspin=None,evecsCV=None,vary_fmax=False,maxmass=None):
267 out = []
268 valsF = get_random_mass_slimline(numiters,minmass1,maxmass1,minmass2,maxmass2,maxspin,maxBHspin = maxBHspin,return_spins=True,maxmass=maxmass)
269 valsF = numpy.array(valsF)
270 mass = valsF[0]
271 eta = valsF[1]
272 beta = valsF[2]
273 sigma = valsF[3]
274 gamma = valsF[4]
275 chis = 0.5*(valsF[5] + valsF[6])
276 if covary:
277 lambdas = get_cov_params(mass,eta,beta,sigma,gamma,chis,f0,evecs,evals,evecsCV,order)
278 else:
279 lambdas = get_conv_params(mass,eta,beta,sigma,gamma,chis,f0,evecs,evals,order,vary_fmax=vary_fmax)
280
281 return numpy.array(lambdas)
282
283 -def get_random_mass_slimline(N,minmass1,maxmass1,minmass2,maxmass2,maxspin,maxBHspin = None,return_spins=False,qm_scalar_fac=1,maxmass=None):
284
285 minmass = minmass1 + minmass2
286 if not maxmass:
287 maxmass = maxmass1 + maxmass2
288 mincompmass = minmass2
289 maxcompmass = maxmass1
290
291 mass = numpy.random.random(N) * (minmass**(-5./3.)-maxmass**(-5./3.)) + maxmass**(-5./3.)
292 mass = mass**(-3./5.)
293 maxmass2 = numpy.minimum(mass/2.,maxmass2)
294 minmass1 = numpy.maximum(minmass1,mass/2.)
295 mineta = numpy.maximum(mincompmass * (mass-mincompmass)/(mass*mass), maxcompmass*(mass-maxcompmass)/(mass*mass))
296 maxeta = numpy.minimum(0.25,maxmass2 * (mass - maxmass2) / (mass*mass))
297 maxeta = numpy.minimum(maxeta,minmass1 * (mass - minmass1) / (mass*mass))
298 if (maxeta < mineta).any():
299 print "WARNING: Max eta is smaller than min eta!!"
300 eta = numpy.random.random(N) * (maxeta - mineta) + mineta
301 diff = (mass*mass * (1-4*eta))**0.5
302 mass1 = (mass + diff)/2.
303 mass2 = (mass - diff)/2.
304 if (mass1 > maxmass1).any() or (mass1 < minmass1).any():
305 print "WARNING: Mass1 outside of mass range"
306 if (mass2 > maxmass2).any() or (mass2 < minmass2).any():
307 print "WARNING: Mass1 outside of mass range"
308 if maxspin > 0:
309 mspin = numpy.zeros(len(mass1))
310 mspin += maxspin
311 if maxBHspin:
312 mspin[mass1 > 3] = maxBHspin
313 spin1z = numpy.random.random(N) * mspin*2 - mspin
314 mspin = numpy.zeros(len(mass2))
315 mspin += maxspin
316 if maxBHspin:
317 mspin[mass2 > 3] = maxBHspin
318 spin2z = numpy.random.random(N) * mspin*2 - mspin
319
320 spinspin = spin1z*spin2z
321 else:
322 spinspin = numpy.zeros(N,dtype=float)
323 spin1z = numpy.zeros(N,dtype=float)
324 spin2z = numpy.zeros(N,dtype=float)
325
326 chiS = 0.5 * (spin1z + spin2z)
327 chiA = 0.5 * (spin1z - spin2z)
328 delta = (mass1 - mass2) / (mass1 + mass2)
329
330 beta = (113. / 12. - 19./3. * eta) * chiS
331 beta += 113. / 12. * delta * chiA
332 sigma = eta / 48. * (474 * spinspin)
333 gamma = numpy.zeros(len(sigma))
334 for spinA,massA in zip([spin1z,spin2z],[mass1,mass2]):
335 sigmaFac = 1. / 96. * (massA / mass)**2
336 sigmaFac2 = (720 * qm_scalar_fac -1) * spinA * spinA
337 sigmaFac3 = (240 * qm_scalar_fac - 7) * spinA * spinA
338 sigma += sigmaFac * (sigmaFac2 - sigmaFac3)
339 gamma = (732985./2268. - 24260./81.*eta - 340./9.*eta*eta)*chiS
340 gamma += (732985. / 2268. + 140./9. * eta) * delta * chiA
341
342 if return_spins:
343 return mass,eta,beta,sigma,gamma,spin1z,spin2z
344 else:
345 return mass,eta,beta,sigma,gamma
346
348 temp = 0
349 for i in range(length):
350 temp += evecs[i,index] * old_vector[i]
351 temp *= rescale_factor
352 return temp
353
355 temp = 0
356 for i in range(length):
357 temp += evecs[index,i] * old_vector[i]
358 temp *= rescale_factor
359 return temp
360
361 -def get_conv_params(totmass,eta,beta,sigma,gamma,chis,f0,evecs,evals,order,vary_fmax=False):
362
363 lambdas = get_chirp_params(totmass,eta,beta,sigma,gamma,chis,f0,order)
364
365 lams = []
366 if not vary_fmax:
367 length = len(evals)
368 for i in range(length):
369 lams.append(rotate_vector(evecs,lambdas,math.sqrt(evals[i]),i,length))
370 return lams
371 else:
372
373 fs = numpy.array(evals.keys(),dtype=float)
374 fs.sort()
375
376 fISCO = (1/6.)**(3./2.) / (LAL_PI * totmass * LAL_MTSUN_SI)
377
378
379 length = len(evals[fs[0]])
380 output=numpy.zeros([length,len(totmass)])
381 lambdas = numpy.array(lambdas)
382
383 for i in range(len(fs)):
384 if (i == 0):
385 logicArr = fISCO < ((fs[0] + fs[1])/2.)
386 if (i == (len(fs)-1)):
387 logicArr = fISCO > ((fs[-1] + fs[-1])/2.)
388 else:
389 logicArrA = fISCO > ((fs[i-1] + fs[i])/2.)
390 logicArrB = fISCO < ((fs[i] + fs[i+1])/2.)
391 logicArr = numpy.logical_and(logicArrA,logicArrB)
392 if logicArr.any():
393 for j in range(length):
394 output[j,logicArr] = rotate_vector(evecs[fs[i]],lambdas[:,logicArr],math.sqrt(evals[fs[i]][j]),j,length)
395
396 for i in range(length):
397 lams.append(output[i])
398 return lams
399
401 lams = []
402 length = len(evals)
403 for i in range(length):
404 lams.append(rotate_vector(evecs,lambdas,math.sqrt(evals[i]),i,length))
405 return lams
406
407 -def get_cov_params(totmass,eta,beta,sigma,gamma,chis,f0,evecs,evals,evecsCV,order):
408 mus = get_conv_params(totmass,eta,beta,sigma,gamma,chis,f0,evecs,evals,order)
409 xis = get_covaried_params(mus,evecsCV)
410 return xis
411
413 length = len(evecsCV)
414 lams = []
415 for i in range(length):
416 lams.append(rotate_vector(evecsCV,lambdas,1.,i,length))
417 return lams
418
420
421 totmass = totmass * mtsun
422 pi = numpy.pi
423 lambda0 = 3 / (128 * eta * (pi * totmass * f0)**(5/3))
424 lambda2 = 5 / (96 * pi * eta * totmass * f0) * (743/336 + 11/4 * eta)
425 lambda3 = (-3 * pi**(1/3))/(8 * eta * (totmass*f0)**(2/3)) * (1 - beta/ (4 * pi))
426 lambda4 = 15 / (64 * eta * (pi * totmass * f0)**(1/3)) * (3058673/1016064 + 5429/1008 * eta + 617/144 * eta**2 - sigma)
427 if order == 'twoPN':
428 return lambda0,lambda2,lambda3,lambda4
429 elif order[0:16] == 'threePointFivePN' or order[0:8] == 'taylorF4':
430 lambda5 = 3. * (38645.*pi/756. - 65.*pi*eta/9. - gamma)
431 lambda5 = lambda5 * (3./(128.*eta))
432 lambda6 = 11583231236531./4694215680. - (640.*pi*pi)/3. - (6848.*LAL_GAMMA)/21.
433 lambda6 -= (6848./21.) * numpy.log(4 * (pi * totmass * f0)**(1./3.))
434 lambda6 += (-15737765635/3048192. + 2255.*pi*pi/12.)*eta
435 lambda6 += (76055.*eta*eta)/1728. - (127825.*eta*eta*eta)/1296.;
436 lambda6 = lambda6
437 lambda6 = lambda6 * 3./(128.*eta) * (pi * totmass * f0)**(1/3.)
438 lambda7 = (77096675.*pi)/254016. + (378515.*pi*eta)/1512. - (74045.*pi*eta*eta)/756.
439 lambda7 = lambda7
440 lambda7 = lambda7 * 3./(128.*eta) * (pi * totmass * f0)**(2/3.)
441 lambda6log = -( 6848./21)
442 lambda6log = lambda6log * 3./(128.*eta) * (pi * totmass * f0)**(1/3.)
443 if order[0:16] == 'threePointFivePN':
444 return lambda0,lambda2,lambda3,lambda4,lambda5,lambda6,lambda7,lambda6log
445 elif order[0:13] == 'taylorF4_45PN' or order[0:13] == 'taylorF4_35PN':
446
447 lambda6spin = 502.6548245743669 * beta + 88.45238095238095 * sigma + \
448 (110. * eta * sigma) - 20. * beta * beta
449 lambda6spin = lambda6spin * 3./(128.*eta) * (pi * totmass * f0)**(1/3.)
450 lambda6 += lambda6spin
451
452 lambda7spin = -510.0603994773098*beta - 368.01525846326734*beta*eta + \
453 1944.363555525455*chis*eta - 502.6548245743669*sigma + \
454 40.*beta*sigma + 241.47615535889872*beta*eta*eta + \
455 2961.654024441635*chis*eta*eta + 676.0619469026549*chis*eta*eta*eta
456 lambda7spin = lambda7spin * 3./(128.*eta) * (pi * totmass * f0)**(2/3.)
457 lambda7 += lambda7spin
458
459 lambda8 = 342.6916926002232 + 2869.024558661873*eta - \
460 3773.659169914512*eta*eta + 172.0609149438239*eta*eta*eta - \
461 24.284336419753085*eta*eta*eta*eta
462 lambda8log = -1028.0750778006693 - 8607.073675985623*eta + \
463 11320.977509743536*eta*eta - 516.1827448314717*eta*eta*eta + \
464 72.85300925925927*eta*eta*eta*eta
465 lambda8loglog = 480.7316704459562 + 597.8412698412699*eta
466
467 lambda8spin = 936.7471880419974*beta - 311.03929625364435*beta*eta - \
468 2455.4171568883194*chis*eta + 195.39588893571195*beta*chis*eta + \
469 48.491201779065534*sigma + 101.92901234567901*eta*sigma - \
470 58.81315259633844*beta*beta + \
471 8.918387413962636*eta*beta*beta - 686.5167663065837*chis*eta*eta \
472 + 54.631268436578175*beta*chis*eta*eta + \
473 71.69753086419753*sigma*eta*eta - \
474 4.444444444444445*sigma*sigma
475 lambda8logspin = -2810.241564125992*beta + 933.117888760933*beta*eta + \
476 7366.251470664957*chis*eta - 586.1876668071359*beta*chis*eta - \
477 145.4736053371966*sigma - 305.78703703703707*eta*sigma \
478 + 176.4394577890153*beta*beta - \
479 26.755162241887906*eta*beta*beta + \
480 2059.5502989197507*chis*eta*eta - \
481 163.89380530973452*beta*chis*eta*eta - \
482 215.0925925925926*sigma*eta*eta + \
483 13.333333333333334*sigma*sigma
484
485 lambda8 = (lambda8 + lambda8spin) * 3./(128.*eta) * \
486 (pi * totmass * f0)**(3./3.)
487 lambda8log = (lambda8log + lambda8logspin) * 3./(128.*eta) * \
488 (pi * totmass * f0)**(3./3.)
489 lambda8loglog = lambda8loglog * 3./(128.*eta) * \
490 (pi * totmass * f0)**(3./3.)
491
492 lambda9 = 20021.24571514093 - 42141.161261993766*eta - \
493 4047.211701119762*eta*eta - 2683.4848475303043*eta*eta*eta
494 lambda9log = -4097.833617482457
495
496 lambda9spin = 2105.9471080635244*beta + 3909.271818583914*beta*eta - \
497 2398.354686411564*chis*eta + 1278.4225104920606*sigma - \
498 198.6688790560472*beta*sigma + 589.0486225480862*eta*sigma - \
499 62.43362831858406*beta*eta*sigma + \
500 439.6407501053519*chis*eta*sigma - 376.99111843077515*beta*beta + \
501 10.*beta*beta*beta - 202.1451909795383*beta*eta*eta - \
502 5711.929102446965*chis*eta*eta + \
503 122.9203539823009*chis*sigma*eta*eta - \
504 493.00738145963066*beta*eta*eta*eta - \
505 4955.659484448894*chis*eta*eta*eta - \
506 991.4721607669617*chis*eta*eta*eta*eta
507 lambda9logspin = 326.0952380952381*beta
508 lambda9 = (lambda9 + lambda9spin) * 3./(128.*eta) * \
509 (pi * totmass * f0)**(4./3.)
510 lambda9log = (lambda9log + lambda9logspin) * 3./(128.*eta) * \
511 (pi * totmass * f0)**(4./3.)
512 if order[0:13] == 'taylorF4_45PN':
513 return lambda0,lambda2,lambda3,lambda4,lambda5,lambda6,lambda7,lambda6log,lambda8log,lambda8loglog,lambda9,lambda9log
514 elif order[0:13] == 'taylorF4_35PN':
515 return lambda0,lambda2,lambda3,lambda4,lambda5,lambda6,lambda7,lambda6log
516 else:
517 raise BrokenError
518 else:
519 raise BrokenError
520
521 -def make_plots(a,b,c,d,aname,bname,cname,dname,paper_plots=False):
522 if paper_plots:
523 paper_plot()
524 if not os.path.isdir('plots'):
525 os.makedirs('plots')
526 vals = [a,b,c,d]
527 names = [aname,bname,cname,dname]
528 for i in range(4):
529 for j in range(i+1,4):
530 pylab.plot(vals[i],vals[j],'b.')
531 pylab.xlabel(names[i])
532 pylab.ylabel(names[j])
533
534
535
536
537
538
539
540
541 pylab.savefig('plots/%s_vs_%s.png' % (names[i],names[j]))
542 pylab.clf()
543
544
545
546
547
548
549
550
551
552
553
554
555
556
558
559 v1s = [minv1]
560 v2s = [minv2]
561 initPoint = [minv1,minv2]
562
563 initLine = [initPoint]
564 tmpv1 = minv1
565 while (tmpv1 < maxv1):
566 tmpv1 = tmpv1 + (3 * mindist)**(0.5)
567 initLine.append([tmpv1,minv2])
568 v1s.append(tmpv1)
569 v2s.append(minv2)
570 initLine = numpy.array(initLine)
571 initLine2 = copy.deepcopy(initLine)
572 initLine2[:,0] += 0.5 * (3*mindist)**0.5
573 initLine2[:,1] += 1.5 * (mindist)**0.5
574 for i in xrange(len(initLine2)):
575 v1s.append(initLine2[i,0])
576 v2s.append(initLine2[i,1])
577 tmpv2_1 = initLine[0,1]
578 tmpv2_2 = initLine2[0,1]
579 while tmpv2_1 < maxv2 and tmpv2_2 < maxv2:
580 tmpv2_1 = tmpv2_1 + 3.0 * (mindist)**0.5
581 tmpv2_2 = tmpv2_2 + 3.0 * (mindist)**0.5
582 initLine[:,1] = tmpv2_1
583 initLine2[:,1] = tmpv2_2
584 for i in xrange(len(initLine)):
585 v1s.append(initLine[i,0])
586 v2s.append(initLine[i,1])
587 for i in xrange(len(initLine2)):
588 v1s.append(initLine2[i,0])
589 v2s.append(initLine2[i,1])
590 v1s = numpy.array(v1s)
591 v2s = numpy.array(v2s)
592 return v1s,v2s
593
595 import lal
596 tiling = lal.CreateFlatLatticeTiling(3)
597 lal.SetFlatLatticeConstantBound(tiling,0,minv1,maxv1)
598 lal.SetFlatLatticeConstantBound(tiling,1,minv2,maxv2)
599 lal.SetFlatLatticeConstantBound(tiling,2,minv3,maxv3)
600 lal.SetFlatLatticeGenerator(tiling,lal.AnstarLatticeGeneratorPtr)
601
602 a = lal.gsl_matrix(3,3)
603 a.data[0,0] = 1
604 a.data[1,1] = 1
605 a.data[2,2] = 1
606 lal.SetFlatLatticeMetric(tiling,a,mindist)
607
608 vs1 = []
609 vs2 = []
610 vs3 = []
611 count = 0
612 while (lal.NextFlatLatticePoint(tiling) >= 0):
613 count += 1
614 if not (count % 100000):
615 print "Now %d points" %(count)
616 p = lal.GetFlatLatticePoint(tiling)
617 vs1.append(p.data[0])
618 vs2.append(p.data[1])
619 vs3.append(p.data[2])
620 return vs1,vs2,vs3
621
622 -def get_physical_covaried_masses(xis,bestMasses,bestXis,f0,temp_number,req_match,order,evecs,evals,evecsCV,maxmass1,minmass1,maxmass2,minmass2,maxNSspin,maxBHspin,return_smaller=False,nsbh_flag=False):
623
624 origScaleFactor = 1
625
626
627 xi_size = len(xis)
628 scaleFactor = origScaleFactor
629 bestChirpmass = bestMasses[0] * (bestMasses[1])**(3./5.)
630 count = 0
631 unFixedCount = 0
632 currDist = 100000000000000000
633 while(1):
634
635
636
637 if count:
638 if currDist > 1 and scaleFactor == origScaleFactor:
639 scaleFactor = origScaleFactor*10
640 chirpmass,totmass,eta,spin1z,spin2z,diff,mass1,mass2,beta,sigma,gamma,chis,new_xis = get_mass_distribution(bestChirpmass,bestMasses[1],bestMasses[2],bestMasses[3],scaleFactor,order,evecs,evals,evecsCV,maxmass1,minmass1,minmass2,maxmass2,maxNSspin,maxBHspin,f0,nsbh_flag = nsbh_flag)
641 cDist = (new_xis[0] - xis[0])**2
642 for j in range(1,xi_size):
643 cDist += (new_xis[j] - xis[j])**2
644 if (cDist.min() < req_match):
645 idx = cDist.argmin()
646 scaleFactor = origScaleFactor
647 return mass1[idx],mass2[idx],spin1z[idx],spin2z[idx],count,cDist.min(),new_xis[0][idx],new_xis[1][idx],new_xis[2][idx],new_xis[3][idx]
648 if (cDist.min() < currDist):
649 idx = cDist.argmin()
650 bestMasses[0] = totmass[idx]
651 bestMasses[1] = eta[idx]
652 bestMasses[2] = spin1z[idx]
653 bestMasses[3] = spin2z[idx]
654 bestChirpmass = bestMasses[0] * (bestMasses[1])**(3./5.)
655 currDist = cDist.min()
656 unFixedCount = 0
657 scaleFactor = origScaleFactor
658 count += 1
659 unFixedCount += 1
660 if unFixedCount > 5000:
661 if return_smaller:
662 diff = (bestMasses[0]*bestMasses[0] * (1-4*bestMasses[1]))**0.5
663 mass1 = (bestMasses[0] + diff)/2.
664 mass2 = (bestMasses[0] - diff)/2.
665 return mass1,mass2,bestMasses[2],bestMasses[3],count,currDist,new_xis[0][0],new_xis[1][0],new_xis[2][0],new_xis[3][0]
666
667 else:
668 raise BrokenError
669 if not unFixedCount % 100:
670 scaleFactor *= 2
671 if scaleFactor > 64:
672 scaleFactor = 1
673
674 raise BrokenError
675
676 -def get_mass_distribution(bestChirpmass,bestEta,bestSpin1z,bestSpin2z,scaleFactor,order,evecs,evals,evecsCV,maxmass1,minmass1,minmass2,maxmass2,maxNSspin,maxBHspin,f0,nsbh_flag = False):
677 chirpmass = bestChirpmass * ( 1 - (numpy.random.random(100) - 0.5) * 0.0001 * scaleFactor)
678 chirpmass[0] = bestChirpmass
679 eta = bestEta * ( 1 - (numpy.random.random(100) - 0.5) * 0.01 * scaleFactor)
680 eta[0] = bestEta
681
682 eta[eta > 0.25] = 0.25
683 eta[eta < 0.0001] = 0.0001
684 totmass = chirpmass / (eta**(3./5.))
685 spin1z = bestSpin1z + ( (numpy.random.random(100) - 0.5) * 0.01 * scaleFactor)
686 spin1z[0] = bestSpin1z
687 spin2z = bestSpin2z + ( (numpy.random.random(100) - 0.5) * 0.01 * scaleFactor)
688 spin2z[0] = bestSpin2z
689 beta,sigma,gamma,chis = get_beta_sigma_from_aligned_spins(totmass,eta,spin1z,spin2z)
690
691 diff = (totmass*totmass * (1-4*eta))**0.5
692 mass1 = (totmass + diff)/2.
693 mass2 = (totmass - diff)/2.
694
695 numploga1 = numpy.logical_and(abs(spin1z) < maxBHspin,mass1 > 2.99)
696 if nsbh_flag:
697 numploga = numploga1
698 else:
699 numploga2 = numpy.logical_and(abs(spin1z) < maxNSspin,mass1 < 3.01)
700 numploga = numpy.logical_or(numploga1,numploga2)
701 numplogb2 = numpy.logical_and(abs(spin2z) < maxNSspin,mass2 < 3.01)
702 if nsbh_flag:
703 numplogb = numplogb2
704 else:
705 numplogb1 = numpy.logical_and(abs(spin2z) < maxBHspin,mass2 > 2.99)
706 numplogb = numpy.logical_or(numplogb1,numplogb2)
707 numplog1 = numpy.logical_and(numploga,numplogb)
708 numplog = numpy.logical_not(numplog1)
709 beta[numplog] = 0
710 sigma[numplog] = 0
711 gamma[numplog] = 0
712 chis[numplog] = 0
713 spin1z[numplog] = 0
714 spin2z[numplog] = 0
715
716 totmass[mass1 < minmass1] = 0.0001
717 totmass[mass1 > maxmass1] = 0.0001
718 totmass[mass2 < minmass2] = 0.0001
719 totmass[mass2 > maxmass2] = 0.0001
720
721 new_xis = get_cov_params(totmass,eta,beta,sigma,gamma,chis,f0,evecs,evals,evecsCV,order)
722 return chirpmass,totmass,eta,spin1z,spin2z,diff,mass1,mass2,beta,sigma,gamma,chis,new_xis
723
724 -def stack_xi_direction_brute(xis,bestMasses,bestXis,f0,temp_number,direction_num,req_match,order,evecs,evals,evecsCV,maxmass1,minmass1,maxmass2,minmass2,maxNSspin,maxBHspin,nsbh_flag=False):
725
726 origScaleFactor = 0.8
727 numTestPoints = 3000
728
729
730 xi_size = len(xis)
731 origMasses = copy.deepcopy(bestMasses)
732 bestChirpmass = bestMasses[0] * (bestMasses[1])**(3./5.)
733 count = 0
734 unFixedCount = 0
735 xi3min = 10000000000
736 xi3max = -100000000000
737
738 for i in range(numTestPoints):
739
740 scaleFactor = origScaleFactor
741 chirpmass,totmass,eta,spin1z,spin2z,diff,mass1,mass2,beta,sigma,gamma,chis,new_xis = get_mass_distribution(bestChirpmass,bestMasses[1],bestMasses[2],bestMasses[3],scaleFactor,order,evecs,evals,evecsCV,maxmass1,minmass1,minmass2,maxmass2,maxNSspin,maxBHspin,f0,nsbh_flag = nsbh_flag)
742 cDist = (new_xis[0] - xis[0])**2
743 for j in range(1,xi_size):
744 cDist += (new_xis[j] - xis[j])**2
745 redCDist = cDist[cDist < req_match]
746 redXis = (new_xis[direction_num])[cDist < req_match]
747
748 if len(redCDist):
749 new_xis[direction_num][cDist > req_match] = -10000000
750 maxXi3 = (new_xis[direction_num]).max()
751 idx = (new_xis[direction_num]).argmax()
752 if maxXi3 > xi3max:
753 xi3max = maxXi3
754 bestMasses[0] = totmass[idx]
755 bestMasses[1] = eta[idx]
756 bestMasses[2] = spin1z[idx]
757 bestMasses[3] = spin2z[idx]
758 m1 = mass1[idx]
759 m2 = mass2[idx]
760 bestChirpmass = bestMasses[0] * (bestMasses[1])**(3./5.)
761
762 bestMasses = copy.deepcopy(origMasses)
763 bestChirpmass = bestMasses[0] * (bestMasses[1])**(3./5.)
764
765 for i in range(numTestPoints):
766
767 scaleFactor = origScaleFactor
768 chirpmass,totmass,eta,spin1z,spin2z,diff,mass1,mass2,beta,sigma,gamma,chis,new_xis = get_mass_distribution(bestChirpmass,bestMasses[1],bestMasses[2],bestMasses[3],scaleFactor,order,evecs,evals,evecsCV,maxmass1,minmass1,minmass2,maxmass2,maxNSspin,maxBHspin,f0,nsbh_flag = nsbh_flag)
769
770 cDist = (new_xis[0] - xis[0])**2
771 for j in range(1,xi_size):
772 cDist += (new_xis[j] - xis[j])**2
773 redCDist = cDist[cDist < req_match]
774 redXis = (new_xis[direction_num])[cDist < req_match]
775
776 if len(redCDist):
777 new_xis[direction_num][cDist > req_match] = 10000000
778 maxXi3 = (new_xis[direction_num]).min()
779 idx = (new_xis[direction_num]).argmin()
780 if maxXi3 < xi3min:
781 xi3min = maxXi3
782 bestMasses[0] = totmass[idx]
783 bestMasses[1] = eta[idx]
784 bestMasses[2] = spin1z[idx]
785 bestMasses[3] = spin2z[idx]
786 m1 = mass1[idx]
787 m2 = mass2[idx]
788 bestChirpmass = bestMasses[0] * (bestMasses[1])**(3./5.)
789
790 return xi3min,xi3max,xi3max-xi3min
791
793 diff = (mass*mass * (1-4*eta))**0.5
794 mass1 = (mass + diff)/2.
795 mass2 = (mass - diff)/2.
796 chiS = 0.5 * (spin1z + spin2z)
797 chiA = 0.5 * (spin1z - spin2z)
798 delta = (mass1 - mass2) / (mass1 + mass2)
799 spinspin = spin1z*spin2z
800
801 beta = (113. / 12. - 19./3. * eta) * chiS
802 beta += 113. / 12. * delta * chiA
803 sigma = eta / 48. * (474 * spinspin)
804 sigma += (1 - 2*eta) * (81./16. * (chiS*chiS + chiA*chiA))
805 sigma += delta * (81. / 8. * (chiS*chiA))
806 gamma = (732985./2268. - 24260./81.*eta - 340./9.*eta*eta)*chiS
807 gamma += (732985. / 2268. + 140./9. * eta) * delta * chiA
808 return beta,sigma,gamma,chiS
809
811
812 aMass1 = point1[0]
813 aMass2 = point1[1]
814 aSpin1 = point1[2]
815 aSpin2 = point1[3]
816 try:
817 leng = len(aMass1)
818 aArray = True
819 except:
820 aArray = False
821
822 bMass1 = point2[0]
823 bMass2 = point2[1]
824 bSpin1 = point2[2]
825 bSpin2 = point2[3]
826 try:
827 leng = len(bMass1)
828 bArray = True
829 except:
830 bArray = False
831
832 if aArray and bArray:
833 print "I cannot take point1 and point2 as arrays"
834
835 aTotMass = aMass1 + aMass2
836 aEta = (aMass1 * aMass2) / (aTotMass * aTotMass)
837 aCM = aTotMass * aEta**(3./5.)
838
839 bTotMass = bMass1 + bMass2
840 bEta = (bMass1 * bMass2) / (bTotMass * bTotMass)
841 bCM = bTotMass * bEta**(3./5.)
842
843 abeta,asigma,agamma,achis = get_beta_sigma_from_aligned_spins(aTotMass,aEta,aSpin1,aSpin2)
844 bbeta,bsigma,bgamma,bchis = get_beta_sigma_from_aligned_spins(bTotMass,bEta,bSpin1,bSpin2)
845
846 aXis = get_cov_params(aTotMass,aEta,abeta,asigma,agamma,achis,f0,evecs,evals,evecsCV,order)
847 if return_xis and not aArray:
848 xis1 = aXis
849
850 bXis = get_cov_params(bTotMass,bEta,bbeta,bsigma,bgamma,bchis,f0,evecs,evals,evecsCV,order)
851 if return_xis and not bArray:
852 xis2 = bXis
853
854 dist = (aXis[0] - bXis[0])**2
855 for i in range(1,len(aXis)):
856 dist += (aXis[i] - bXis[i])**2
857
858 if aArray and return_xis:
859 aXis = numpy.array(aXis)
860 xis1 = aXis[:,dist.argmin()]
861 if bArray and return_xis:
862 bXis = numpy.array(bXis)
863 xis2 = bXis[:,dist.argmin()]
864
865 if return_xis:
866 return xis1,xis2
867
868 return dist
869