sampFreq = 2.0**(-2.0)                #采样频率(Sampling frequency)，单位时间样本点个数，应大于 2f（即Nyquist频率)
duration = 2.0**16.0         #信号持续时间(duration of signal) 2^16, 0.75d, 2^25 ,1.06yr


n = int(duration * sampFreq)#采样点数(Sampling Number), 有时也称为信号长度(Length of Signal) 2^16
                            #为2的幂时，快速傅里叶变化效率最高
                            #n =  duration * sampFreqint = (duration / sampIntrvl)

sampIntrvl = 1.0 / sampFreq                   #采样周期(Sampling period)，隔多少时间取样一次，或步长
freqIntrvl = sampFreq / n                     #傅里叶变换 频率分辨率(Frequency Interval) 
                                              # freqIntrvl = 1 / duration = 1 / (n * sampIntrvl)
                                              #            = sampFreq / n  
        
f_min = 1e-4                         #可以设置
f_max = sampFreq/2.0             #信号模式的最大频率

print("采样频率为%fHz，信号持续时间%ds, 时域信号采样%d 个点"%(sampFreq,duration,n))
print("信号中可分析最大频率为%fHz"%f_max)
print("\n采样周期，即时域分辨率为%fs"%(sampIntrvl))
print("信号频域的频率间隔，即频域分辨率为%fHz"%freqIntrvl);

#注意： 会有除0警告，可以忽略
freqVec = np.linspace(0, 1, int(1/freqIntrvl))
#构建频率序列(必须等频率间隔)，这里我们开始于0Hz(必须开始于0，否则导入pyCBC的结果不对), 结束于1Hz，频率间隔为上面的freqIntrvl

#LISA PSD 模拟函数
psd = LISAcalPSD1(freqVec) 

'''
1. Babak, S., Fang, H., Gair, J. R., Glampedakis, K. & Hughes, S. A. Kludge’
gravitational waveforms for a test-body orbiting a Kerr black hole. Phys.Rev. D 75, 024005 (2007).

'''

#构建用于pyCBC的FrequencySeries
psd = FrequencySeries(psd, delta_f=freqIntrvl, epoch='', dtype=None, copy=True)

#绘制 频率 - sqrt(PSD) 图 
plt.figure()
plt.loglog(psd.sample_frequencies, np.sqrt(psd), label = 'LISA-psd')
plt.xlim(1e-5,1e0)
plt.ylim(1e-21,1e-15)
plt.xlabel("Frequency / Hz")
plt.ylabel("$\sqrt{S_{n}(f) \ / \ (Hz^{-1})}$")
plt.legend()
plt.show()

#   Parameters

atilde = 0.9  # 自旋参量
nu = 1e-05 # 质量比
M = (1e6 * const.M_sun).to(u.kg).value  #the mass of system, kilogram
# m = (1e1 * const.M_sun).to(u.kg).value  #the mass of small compact body, kilogram

l = np.arange(2,3)      
m = 2
k = np.arange(-2,12,1)                   #the nth harmonic


########################
R_unit = (G * M / C**2)#                #the unit length in G = C = 1 system, metre
t_unit = R_unit / C    #                #the unit time in G = C = 1 system, second
print(R_unit, t_unit)
########################

e = 0.5                                 #eccentricity of orbit
p = 12 * R_unit                         #semi-latus rectum

a = p / (1 - e**2)                      #the semi-major axis of orbit, metre

D = (1.00 * u.Gpc).to('m').value

# #   Time Vector                  
# tb = 17542 *  t_unit  #17542 * t_unit                 #end time, second
# step = t_unit                           #step

step = sampIntrvl

wave = teukolsky(atilde = atilde, nu = nu, M = M, l = l, k = k, e = e, a = a, D = D)
wave.ecalculate(duration = duration, delta_t = step)

# h_plus_Vec = wave.hplus                  # h_plus of gravitational-wave
# h_cross_Vec = wave.hcross                # h_cross of gravitational-wave


hp, hc = get_td_waveform(template=wave)
b1 = time.time()
print("用时%fs"%(b1-a1))

a1 = time.time()

lambda_g = (1.6e13 * u.km).to('m').value#the compton wavelength of graviton
wave.eadd_dispersion(lambda_g)
hpd, hcd = get_td_waveform(template=wave)

b1 = time.time()
print("用时%fs"%(b1-a1))

#绘制波形图 
plt.figure()
plt.plot(hp.sample_times, hp, label = '$h_{+}$')
plt.plot(hp.sample_times, hc, label = '$h_{\\times}$')
plt.xlim(0,8e3)
# plt.ylim(- 2e-22,2e-22)
plt.xlabel("Time / s")
plt.ylabel("$Strain$")
plt.legend()
plt.show()

#绘制 波形图 
plt.figure()
plt.plot(hpd.sample_times, hpd, label = '$h_{+}$')
plt.plot(hpd.sample_times, hcd, label = '$h_{\\times}$')
plt.xlim(0,8e3)
# plt.ylim(- 2e-22,2e-22)
plt.xlabel("Time / s")
plt.ylabel("$Strain$")
plt.legend()
plt.show()

overlap = overlap_func(hc, hcd, psd, sampIntrvl, f_min , f_max)
overlap

match = match_func(hc, hcd, psd, sampIntrvl, f_min , f_max)
match