引力波数据分析第六章题目

发表于 更新于 分类于 Waline： 阅读次数：

本系列文章旨在充实博客内容记录引力波数据分析课程的作业成果，标记为引用的题目部分是直接复制自胡教授引力波数据处理讲义的内容，大部分公式也来自该讲义，不再一一标出。

如果真的有人搜索到了本系列文章并参考，最好能标记一下引用哦。

Exercise 6.1

作真值表，证明逻辑上 ，亦即逻辑与的分配律。

作真值表如下图，可知 ，即逻辑与的分配率。

A B C A,B A,C B+C A,B+A,C A,(B+C)
T T T T T T T T
T T F T F T T T
T F T F T T T T
T F F F F F F F
F T T F F T F F
F T F F F T F F
F F T F F T F F
F F F F F F F F

Exercise 6.2

给定特定引力波源以及对应的物理参数，可以计算出一个引力波探测器对此类源的极限探测距离，这一距离通常被定义为视界距离 horizon distance。对应的情形为，信号处于探测器坐标系中的天顶位置处，且轨道倾角取 ，极化角 。换言之，响应后的信号为 。通常情况下，选取单个 Michelson 干涉仪的信噪比 作为探测阈值。

请编写程序，利用非全天平均的天琴灵敏度曲线（见公式4.99），计算不同质量双黑
洞并合事件的视界距离。具体画图中，横坐标取总质量 ，纵坐标取视界距离 ，注意
请以 log-log 形式画图。计算过程中，假设质量比为 1 : 1，假设对所有的并合事件观测时长均为 3 个月（即，天琴名义连续工作时长，且双黑洞在天琴结束观测的瞬间刚好并合）。

天琴的探测器臂长 ，加速度噪声和位置噪声指标分别是 和 ，传递频率 。

其非全天平均灵敏度曲线，可以表示为

根据之前的知识，有

其中啁啾质量，根据题目假定质量比为 ，那么有。

考虑满足最优信噪比的条件，有

按照题目选定的单 Michelson 干涉仪信噪比 作为阈值，频率下限按照观测时间为3个月计算，对于给定的总质量 ，代入上式积分，把视界距离 提出来，即可求得此时对应的视界距离。

总质量范围取 ，可画出不同双黑洞并和事件的视界距离，如下图所示。绘图代码在最后给出。

不同双黑洞并和事件的视界距离

Exercise 6.3

误警率和探测率分别为

其中 服从自由度为 4 的中心 分布，而对于信噪比为 的信号而言， 服从自由度为 4，非中心参数为 的非中心 分布。可由单次测量误警率 ，代入计算出对应的 ，再代入探测率公式中，即可得到 时，信噪比 。

Exercise 6.4

根据题目参数空间范围是 ，，同时使用归一化条件 确定剩余参数距离 的取值。

使用数值差分计算出参数空间上定义的度规

注意到采取不同的差分步长会影响到最终的结果，这里在参数空间长度 范围内尝试步长，最后确定取变化较小的步长值，参数空间长度的 。

这里内积计算时的积分上限采用最内稳定圆轨道对应的频率

根据题目参数空间范围是 ，，可以使用数值差分计算出参数空间上定义的度规

之后使用蒙特卡洛方法，计算参数空间内的积分，即可估计波形模板的数量，其中 。

最后计算得到的 。

计算代码

Exercise 6.2

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
import numpy as np
import matplotlib.pyplot as plt
import scipy.constants as cons
import scipy.integrate as integrate

# constants
# sun mass
sun = 1.98855e30

mpc = 1e6 * cons.parsec

# length of tianqin
length = np.sqrt(3) * 1e8

# tianqin run time
run_time = 90 * cons.day

def chirp(mass):
    return np.power(1/2, 3/5) * mass

def frequency(t, mass):
    return np.power(5,3/8)/(8*cons.pi)*np.power(np.power(cons.c,3)/(cons.G*chirp(mass)),5/8)*np.power((-t),-3/8)

def min_f(mass):
    return frequency(-run_time, mass)

def hf(f, mass):
    return np.sqrt(5/384)*np.power(cons.G,5/6)*np.power(cons.c,-3/2)*np.power(chirp(mass),5/6)*np.power(np.pi,-2/3)*np.power(f,-7/6)

def sensitive(f):
    return (4*1e-30*(1+1e-4/f)/np.power((2*np.pi*f), 4)+1e-24)*(1+0.6*np.power((f/0.28),2))/np.power(length, 2)

def inner(f, mass):
    return np.power(hf(f, mass),2)/sensitive(f)

def inter(mass):
    v, _ = integrate.quad(inner, min_f(mass), np.inf, args=(mass))
    return np.sqrt(4 * v / 64)/ mpc

mass = np.logspace(0, 10, 5000)
distance = [dis for dis in map(inter, mass*sun)]

plt.figure(figsize = (10, 6))
plt.loglog(mass, distance)
plt.xlabel(r'Total Mass($M_\odot$)')
plt.ylabel(r'Horizon Distance($Mpc$)')
plt.show()

Exercise 6.4

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
import numpy as np
import scipy.constants as cons
import scipy.integrate as integrate
import multiprocessing as mp
import time

# constants
# sun mass
sun = 1.98855e30
mpc = 1e6 * cons.parsec

# length of tianqin
length = np.sqrt(3) * 1e8

def chirp(mass1, mass2):
    return np.power(mass1*mass2, 3/5) / np.power(mass1+mass2, 1/5)

def frequency(tc, mass1, mass2):
    return np.power(5,3/8)/(8*cons.pi)*np.power(np.power(cons.c,3)/(cons.G*chirp(mass1, mass2)),5/8)*np.power((-tc),-3/8)

def min_f(mass1, mass2, tc):
    return frequency(-tc, mass1, mass2)

def hf(f, mass1, mass2, tc, dis):
    phif = 2*np.pi*f*tc-cons.pi/4+3/4*np.power(8*cons.pi*cons.G*chirp(mass1, mass2)*f/np.power(cons.c,3),-5/3)
    return np.sqrt(5/384)*np.power(cons.G,5/6)*np.power(cons.c,-3/2)*np.power(chirp(mass1, mass2),5/6)*np.power(np.pi,-2/3)*np.power(f,-7/6)*np.exp(1j*phif)/dis

def sensitive(f):
    return (4*1e-30*(1+1e-4/f)/np.power((2*np.pi*f), 4)+1e-24)*(1+0.6*np.power((f/0.28),2))/np.power(length, 2)

def norm(f, mass1, mass2):
    freq = np.sqrt(5/384)*np.power(cons.G,5/6)*np.power(cons.c,-3/2)*np.power(chirp(mass1, mass2),5/6)*np.power(np.pi,-2/3)*np.power(f,-7/6)
    return np.power(freq, 2) / sensitive(f)

def partial(h1, h2, f, mass1, mass2, tc, dis):
    if h1 == 'm1' and h2 == 'm1':
        return (hf(f, mass1+step1, mass2, tc, dis) + hf(f, mass1-step1, mass2, tc, dis) - 2*hf(f, mass1, mass2, tc, dis))/np.power(step1, 2), step1, step1
    elif h1 == 'm2' and h2 == 'm2':
        return (hf(f, mass1, mass2+step1, tc, dis) + hf(f, mass1, mass2-step1, tc, dis) - 2*hf(f, mass1, mass2, tc, dis))/np.power(step1, 2), step1, step1
    elif h1 == 'tc' and h2 == 'tc':
        return (hf(f, mass1, mass2, tc+step2, dis) + hf(f, mass1, mass2, tc-step2, dis) - 2*hf(f, mass1, mass2, tc, dis))/np.power(step2, 2), step2, step2
    elif (h1 == 'm1' and h2 == 'm2') or (h1 == 'm2' and h2 == 'm1'):
        return (hf(f, mass1+step1, mass2+step1, tc, dis) + hf(f, mass1-step1, mass2-step1, tc, dis)
            - hf(f, mass1+step1, mass2-step1, tc, dis) - hf(f, mass1-step1, mass2+step1, tc, dis))/np.power(2*step1, 2), step1, step1
    elif (h1 == 'm1' and h2 == 'tc') or (h1 == 'tc' and h2 == 'm1'):
        return (hf(f, mass1+step1, mass2, tc+step2, dis) + hf(f, mass1-step1, mass2, tc-step2, dis)
            - hf(f, mass1+step1, mass2, tc-step2, dis) - hf(f, mass1-step1, mass2, tc+step2, dis))/(2*step1*2*step2), step1, step2
    elif (h1 == 'm2' and h2 == 'tc') or (h1 == 'tc' and h2 == 'm2'):
        return (hf(f, mass1, mass2+step1, tc+step2, dis) + hf(f, mass1, mass2-step1, tc-step2, dis)
            - hf(f, mass1, mass2+step1, tc-step2, dis) - hf(f, mass1, mass2-step1, tc+step2, dis))/(2*step1*2*step2), step1, step2

def inner(f, h1, h2, mass1, mass2, tc, dis):
    par = partial(h1, h2, f, mass1, mass2, tc, dis)
    return (hf(f, mass1, mass2, tc, dis).conjugate()*par[0]/sensitive(f)).real

def inter(h1, h2, mass1, mass2, tc):
    fisco = np.power(cons.c, 3) / (np.power(6, 2/3)*cons.pi*cons.G*(mass1+mass2))
    v, _ = integrate.quad(norm, min_f(mass1, mass2, tc), fisco, args=(mass1, mass2))
    dis = np.sqrt(4 * v)
    v, _ = integrate.quad(inner, min_f(mass1, mass2, tc), fisco, args=(h1, h2, mass1, mass2, tc, dis), limit=100)
    return -2 * v

def martrix(args):
    mass1 = args[0]
    mass2 = args[1]
    tc = args[2]
    gij = [inter(h1, h2, mass1, mass2, tc) for h1 in ['m1', 'm2', 'tc'] for h2 in ['m1', 'm2', 'tc']]
    gij = np.reshape(gij, (3, 3))
    return np.sqrt(np.abs(np.linalg.det(gij)))

# partial step
step1 = 1e-7*95*sun
step2 = 1e-7*5*cons.year

if __name__ == '__main__':
    num = 400
    time1 = time.time()
    mass1 = np.random.uniform(5, 100, num)*sun
    mass2 = np.random.uniform(5, 100, num)*sun
    tc = np.random.uniform(0, 5, num)*cons.year
    with mp.Pool() as p:
        a = p.map(martrix, [(m1, m2, tc) for m1, m2, tc in zip(mass1, mass2, tc)])
    ans = sum(a) * np.power(95*sun, 2) * 5*cons.year / num / np.power(0.03, 3/2)
    print(ans)
    with open('result.txt', 'a') as f:
        f.write(f'{step1:5e}:{ans:.5e}\n')
    print(time.time() - time1)