Problem4.py

import math

import numpy as np
import matplotlib.pyplot as plt
from math import factorial
import scipy.integrate as integrate

# 定义参数
from astropy.time import Time
from scipy.interpolate import interp1d

# 定义参数
T_obs = 10  # 观测时间 (s)
lambda_b = 1.54  # 背景光子流量密度 (ph/(s*cm^2))
lambda_s = 10 * lambda_b  # 脉冲星光子流量密度
A_eff = 250  # 探测器有效面积 (cm^2)

# 脉冲星自转参数
v = 29.647854750036593  # 自转频率 (s^-1)
dv = -368970.96e-15  # 自转频率一阶导数 (s^-2)

# 参考历元
pch_pm = 57715.000000295  # MJD
t0_mjd = Time(pch_pm, format='mjd').jd  # 转换为 Julian Date

# 当前时间，假设计算今天的相位
t_now = Time.now().jd  # 获取当前 Julian Date
delta_t = (t_now - t0_mjd) * 86400  # 时间差 (秒)

# 脉冲星位置变化参数
dt_ra = -14.7  # 赤经自行 (mas/year)
dt_dec = 2.0  # 赤纬自行 (mas/year)

# 转换自行为秒内的相位变化
# 假设每年的角度变化可以近似转换为相位变化
mas_to_rad = np.pi / (180 * 3600 * 1000)  # 角秒转换为弧度
dt_ra_rad_per_sec = dt_ra * mas_to_rad / (365.25 * 86400)  # 赤经自行变化率 (rad/s)
dt_dec_rad_per_sec = dt_dec * mas_to_rad / (365.25 * 86400)  # 赤纬自行变化率 (rad/s)


# 计算相位
def calculate_phase(t, t0=t0_mjd, phi_t0=0):
    delta_t = (t - t0) * 86400  # 转换为秒
    phase = phi_t0 + v * delta_t + 0.5 * dv * delta_t ** 2
    # 加入赤经和赤纬自行导致的相位漂移
    phase += dt_ra_rad_per_sec * delta_t + dt_dec_rad_per_sec * delta_t
    return phase % 1  # 返回周期相位


# 计算当前时间的相位
phase_now = calculate_phase(t_now)
print(f"在当前时间 {Time.now().iso} 时的脉冲相位为: {phase_now}")


def h1(_phi):
    # 数据点列表
    data_points = [
        (0.0019531250, 5.8764781107),
        (0.0058593750, 4.4038634533),
        (0.0097656250, 3.8322155694),
        (0.0136718750, 2.8939674122),
        (0.0175781250, 2.6229006846),
        (0.0214843750, 2.3829017768),
        (0.0253906250, 1.8050403290),
        (0.0292968750, 1.8291178893),
        (0.0332031250, 1.5471774249),
        (0.0371093750, 1.2846543478),
        (0.0410156250, 1.2271788812),
        (0.0449218750, 1.0733931733),
        (0.0488281250, 0.8893163411),
        (0.0527343750, 0.9498985897),
        (0.0566406250, 0.8683455627),
        (0.0605468750, 0.9211608564),
        (0.0644531250, 0.8823260816),
        (0.0683593750, 0.6283466550),
        (0.0722656250, 0.8652387808),
        (0.0761718750, 0.6174729181),
        (0.0800781250, 0.8225205286),
        (0.0839843750, 0.7324238512),
        (0.0878906250, 0.6625212567),
        (0.0917968750, 0.5506771055),
        (0.0957031250, 0.4629105147),
        (0.0996093750, 0.6516475198),
        (0.1035156250, 0.6563076928),
        (0.1074218750, 0.6664047342),
        (0.1113281250, 0.7440742836),
        (0.1152343750, 0.4706774696),
        (0.1191406250, 0.7619382800),
        (0.1230468750, 0.6050457901),
        (0.1269531250, 0.4201922625),
        (0.1308593750, 0.5895118803),
        (0.1347656250, 0.4007748751),
        (0.1386718750, 0.6811619486),
        (0.1425781250, 0.4465999093),
        (0.1464843750, 0.6477640423),
        (0.1503906250, 0.7728120169),
        (0.1542968750, 0.6710649072),
        (0.1582031250, 0.5242694587),
        (0.1621093750, 0.6578610838),
        (0.1660156250, 0.6174729181),
        (0.1699218750, 0.5933953577),
        (0.1738281250, 0.6322301324),
        (0.1777343750, 0.5374732821),
        (0.1816406250, 0.5281529362),
        (0.1855468750, 0.5910652712),
        (0.1894531250, 0.5017452894),
        (0.1933593750, 0.4372795634),
        (0.1972656250, 0.4660172966),
        (0.2011718750, 0.3712604463),
        (0.2050781250, 0.6594144748),
        (0.2089843750, 0.5662110154),
        (0.2128906250, 0.4652406012),
        (0.2167968750, 0.4869880750),
        (0.2207031250, 0.7269869828),
        (0.2246093750, 0.8147535736),
        (0.2285156250, 0.5048520714),
        (0.2324218750, 0.8605786078),
        (0.2363281250, 0.8629086943),
        (0.2402343750, 0.5009685939),
        (0.2441406250, 0.6306767415),
        (0.2480468750, 0.5421334551),
        (0.2519531250, 0.7184433323),
        (0.2558593750, 0.7409675017),
        (0.2597656250, 0.8442680024),
        (0.2636718750, 0.6811619486),
        (0.2675781250, 0.7518412386),
        (0.2714843750, 0.6959191630),
        (0.2753906250, 0.9452384167),
        (0.2792968750, 0.8333942655),
        (0.2832031250, 0.8085400097),
        (0.2871093750, 0.8318408745),
        (0.2910156250, 0.9281511158),
        (0.2949218750, 1.0710630868),
        (0.2988281250, 0.9708693680),
        (0.3027343750, 0.7930060998),
        (0.3066406250, 1.0104808382),
        (0.3105468750, 1.1961110614),
        (0.3144531250, 1.4042654539),
        (0.3183593750, 1.3304793819),
        (0.3222656250, 1.3242658180),
        (0.3261718750, 1.5386337744),
        (0.3300781250, 1.6077596734),
        (0.3339843750, 1.7164970426),
        (0.3378906250, 1.7164970426),
        (0.3417968750, 1.8601857091),
        (0.3457031250, 1.9145543937),
        (0.3496093750, 1.8594090136),
        (0.3535156250, 2.1017380079),
        (0.3574218750, 2.0908642709),
        (0.3613281250, 2.4411539389),
        (0.3652343750, 2.6073667747),
        (0.3691406250, 2.4729984542),
        (0.3730468750, 2.6663956323),
        (0.3769531250, 2.7417350953),
        (0.3808593750, 3.1137722371),
        (0.3847656250, 3.2784316819),
        (0.3886718750, 3.5091102437),
        (0.3925781250, 3.3467808854),
        (0.3964843750, 3.1836748316),
        (0.4003906250, 3.3809554872),
        (0.4042968750, 3.4174601754),
        (0.4082031250, 3.0042581724),
        (0.4121093750, 2.7906669114),
        (0.4160156250, 2.4062026417),
        (0.4199218750, 2.1483397375),
        (0.4238281250, 1.9526124729),
        (0.4277343750, 1.7126135652),
        (0.4316406250, 1.5557210753),
        (0.4355468750, 1.4190226683),
        (0.4394531250, 1.2341691406),
        (0.4433593750, 1.1611597641),
        (0.4472656250, 0.9763062365),
        (0.4511718750, 0.7844624494),
        (0.4550781250, 0.9483451987),
        (0.4589843750, 0.8170836601),
        (0.4628906250, 0.9118405104),
        (0.4667968750, 0.8947532096),
        (0.4707031250, 0.5965021397),
        (0.4746093750, 0.8900930366),
        (0.4785156250, 0.7363073287),
        (0.4824218750, 0.6858221216),
        (0.4863281250, 0.6011623127),
        (0.4902343750, 0.6772784711),
        (0.4941406250, 0.4015515706),
        (0.4980468750, 0.5770847523),
        (0.5019531250, 0.4955317254),
        (0.5058593750, 0.4442698228),
        (0.5097656250, 0.2695133366),
        (0.5136718750, 0.4714541651),
        (0.5175781250, 0.3968913976),
        (0.5214843750, 0.4768910336),
        (0.5253906250, 0.4116486120),
        (0.5292968750, 0.3316489761),
        (0.5332031250, 0.4023282661),
        (0.5371093750, 0.3891244427),
        (0.5410156250, 0.2935908969),
        (0.5449218750, 0.4497066913),
        (0.5488281250, 0.3805807923),
        (0.5527343750, 0.2073776970),
        (0.5566406250, 0.3634934914),
        (0.5605468750, 0.4132020030),
        (0.5644531250, 0.1918437871),
        (0.5683593750, 0.2726201185),
        (0.5722656250, 0.4372795634),
        (0.5761718750, 0.3844642697),
        (0.5800781250, 0.2361154303),
        (0.5839843750, 0.3122315888),
        (0.5878906250, 0.1545624034),
        (0.5917968750, 0.3417460175),
        (0.5957031250, 0.4093185256),
        (0.5996093750, 0.2858239419),
        (0.6035156250, 0.2213582159),
        (0.6074218750, 0.3013578518),
        (0.6113281250, 0.2174747384),
        (0.6152343750, 0.1801933547),
        (0.6191406250, 0.3549498409),
        (0.6230468750, 0.2998044608),
        (0.6269531250, 0.3829108787),
        (0.6308593750, 0.1561157944),
        (0.6347656250, 0.3751439238),
        (0.6386718750, 0.3013578518),
        (0.6425781250, 0.0000000000),
        (0.6464843750, 0.2687366411),
        (0.6503906250, 0.2788336825),
        (0.6542968750, 0.2260183889),
        (0.6582031250, 0.3533964500),
        (0.6621093750, 0.4396096498),
        (0.6660156250, 0.2928142014),
        (0.6699218750, 0.2974743744),
        (0.6738281250, 0.1988340466),
        (0.6777343750, 0.2027175240),
        (0.6816406250, 0.1304848431),
        (0.6855468750, 0.5118423308),
        (0.6894531250, 0.3914545292),
        (0.6933593750, 0.3324256716),
        (0.6972656250, 0.3743672283),
        (0.7011718750, 0.3099015023),
        (0.7050781250, 0.2966976789),
        (0.7089843750, 0.2928142014),
        (0.7128906250, 0.2205815204),
        (0.7167968750, 0.3401926266),
        (0.7207031250, 0.3106781978),
        (0.7246093750, 0.2780569870),
        (0.7285156250, 0.4248524354),
        (0.7324218750, 0.4077651346),
        (0.7363281250, 0.3067947203),
        (0.7402343750, 0.1599992718),
        (0.7441406250, 0.3510663635),
        (0.7480468750, 0.3751439238),
        (0.7519531250, 0.3580566229),
        (0.7558593750, 0.2011641331),
        (0.7597656250, 0.2493192537),
        (0.7636718750, 0.4349494769),
        (0.7675781250, 0.2198048249),
        (0.7714843750, 0.2897074194),
        (0.7753906250, 0.3953380067),
        (0.7792968750, 0.2337853438),
        (0.7832031250, 0.3324256716),
        (0.7871093750, 0.3215519347),
        (0.7910156250, 0.3114548933),
        (0.7949218750, 0.4132020030),
        (0.7988281250, 0.3813574878),
        (0.8027343750, 0.3332023671),
        (0.8066406250, 0.3782507058),
        (0.8105468750, 0.3277654986),
        (0.8144531250, 0.2998044608),
        (0.8183593750, 0.2772802915),
        (0.8222656250, 0.2143679565),
        (0.8261718750, 0.2780569870),
        (0.8300781250, 0.3891244427),
        (0.8339843750, 0.2563095132),
        (0.8378906250, 0.2330086483),
        (0.8417968750, 0.2524260357),
        (0.8457031250, 0.4209689580),
        (0.8496093750, 0.2827171600),
        (0.8535156250, 0.2733968140),
        (0.8574218750, 0.1638827493),
        (0.8613281250, 0.3370858446),
        (0.8652343750, 0.1615526628),
        (0.8691406250, 0.3697070553),
        (0.8730468750, 0.3976680931),
        (0.8769531250, 0.3751439238),
        (0.8808593750, 0.2213582159),
        (0.8847656250, 0.3091248068),
        (0.8886718750, 0.4706774696),
        (0.8925781250, 0.5436868461),
        (0.8964843750, 0.5902885757),
        (0.9003906250, 0.7231035053),
        (0.9042968750, 0.4256291309),
        (0.9082031250, 0.5576673650),
        (0.9121093750, 0.5739779704),
        (0.9160156250, 0.6532009108),
        (0.9199218750, 0.7992196637),
        (0.9238281250, 0.9871799734),
        (0.9277343750, 0.9398015482),
        (0.9316406250, 1.1634898506),
        (0.9355468750, 1.1029076021),
        (0.9394531250, 1.4834883943),
        (0.9433593750, 1.3693141566),
        (0.9472656250, 1.6256236698),
        (0.9511718750, 1.8702827505),
        (0.9550781250, 2.0217383719),
        (0.9589843750, 2.2780478851),
        (0.9628906250, 2.6834829332),
        (0.9667968750, 2.7595990916),
        (0.9707031250, 3.1331896244),
        (0.9746093750, 3.6621192561),
        (0.9785156250, 4.2283302716),
        (0.9824218750, 4.7176484330),
        (0.9863281250, 5.5595863490),
        (0.9902343750, 6.3424954073),
        (0.9941406250, 6.9560848479),
        (0.9980468750, 6.8139495724),
    ]
    # 将数据点转换为两个数组
    phi_values = np.array([point[0] for point in data_points])
    h_values = np.array([point[1] for point in data_points])
    # 计算边界点的斜率
    slope_left = (h_values[1] - h_values[0]) / (phi_values[1] - phi_values[0])
    slope_right = (h_values[-1] - h_values[-2]) / (phi_values[-1] - phi_values[-2])

    # 定义线性外推函数
    def linear_extrapolate(_x, x0, y0, slope):
        return y0 + slope * (_x - x0)

    # 计算边界值
    left_extrapolated_value = linear_extrapolate(0, phi_values[0], h_values[0], slope_left)
    right_extrapolated_value = linear_extrapolate(1, phi_values[-1], h_values[-1], slope_right)
    # 创建插值函数
    func_h = interp1d(phi_values, h_values, kind='cubic', bounds_error=False,
                      fill_value=(left_extrapolated_value, right_extrapolated_value))

    return func_h(_phi)


# 归一化脉冲轮廓函数
def h(_phi):
    return h1(_phi)


# 计算流量
def lambda_total(t):
    # _phi = (v * t) % 1  # 计算相位
    _phi = calculate_phase(t + t0_mjd)
    return lambda_b + lambda_s * h(_phi)


# 生成光子到达时间
def simulate_photons():
    arrival_times = []
    time = 0
    while time < T_obs:
        rate = lambda_total(time)
        dt = np.random.exponential(1 / (rate * A_eff))  # 使用指数分布生成光子到达时间间隔
        time += dt
        if time < T_obs:
            arrival_times.append(time)
    return np.array(arrival_times)


# 折叠脉冲轮廓
def fold_pulse_profile(_photon_times):
    phases = (v * _photon_times) % 1  # 计算到达时间的相位
    profile, _ = np.histogram(phases, bins=100, range=(0, 1))
    return profile


# 定义累计函数 Λ(t) 通过数值积分计算
def cumulative_intensity(t):
    return integrate.quad(lambda s: lambda_total(s), 0, t)[0]


# 计算光子到达概率
def poisson_probability(k, t_a, t_b):
    Lambda_diff = cumulative_intensity(t_b) - cumulative_intensity(t_a)
    return (np.exp(-Lambda_diff) * (Lambda_diff ** k)) / math.factorial(k)


# 小问1： 仿真光子到达时间
photon_times = simulate_photons()
# 小问2： 计算并绘制折叠脉冲轮廓
pulse_profile = fold_pulse_profile(photon_times)


# 时间范围和 k 值
t_values = np.linspace(0, 10, 500)  # 时间段 (t_b - t_a) 从 0 到 10
k_values = [0, 1, 2, 5]  # 不同 k 值

# 绘制 P(k; t_a, t_b) 随时间的变化曲线
plt.figure(figsize=(10, 6))

for k in k_values:
    P_values = [poisson_probability(k, 0, t_b) for t_b in t_values]
    plt.plot(t_values, P_values, label=f'k = {k}')

plt.xlabel('t_b - t_a')
plt.ylabel('P(k; t_a, t_b)')
plt.title('Poisson Probability vs Time Interval (t_b - t_a)')
plt.legend()
plt.grid(True)
plt.show()

"""

def gen_photon_sequence(T_obs, lambda_b, lambda_s, A):
    N_b = lambda_b * A * T_obs
    N_s = lambda_s * A * T_obs
    return np.random.poisson(N_b, size=int(N_b)) + np.random.poisson(N_s, size=int(N_s))


v1 = 29.647854750036593
dv1 = -368970.96e-15
# 设定 t_0 为参考历元时刻
t_0 = 57715.000000295
t_0_phi = t_0 * v1 % 1


def phi_t(t):
    return (t_0_phi + v1 * (t - 0) + dv1 * (t - 0) ** 2 / 2) % 1


# 在时间段(t_a,t_b )内，有k个光子到达探测器的概率为
def compute_probability(k, t_a, t_b, lambda_b, lambda_s):
    # 定义速率函数 λ(t)
    def lambda_func(t):
        return lambda_b + lambda_s * ContourCurvesFunc(phi_t(t))

    # 计算累积函数 Λ(t)
    def cumulative_lambda(t):
        # 使用数值积分方法计算 Λ(t)
        # 使用 np.linspace 生成时间间隔
        time_points = np.linspace(0, t, 1000)
        lambda_values = lambda_func(time_points)
        return np.trapz(lambda_values, dx=t / 1000)

    # 计算 Λ(t_b) 和 Λ(t_a)
    lambda_t_b = cumulative_lambda(t_b)
    lambda_t_a = cumulative_lambda(t_a)
    print(lambda_t_b, lambda_t_a)
    # 计算概率 P(k; t_a, t_b)
    numerator = np.exp(-(lambda_t_b - lambda_t_a)) * ((lambda_t_b - lambda_t_a) ** k)
    # K 的 阶乘
    denominator = math.factorial(k)
    print(numerator, denominator)
    return numerator / denominator
# 示例调用
k = 10
t_a = 0.0
t_b = 50.0
lambda_b = 1.54
lambda_s = 15.4

probability = compute_probability(k, t_a, t_b, lambda_b, lambda_s)
print(probability)

# 示例调用
gen = gen_photon_sequence(10, 1.54, 10 * 1.54, 250)

gen[0].shape, gen[1].shape


"""