Skip to content

A
AMSSNCKU
Commit da4ca3ce authored by Chen-Kai Qiao's avatar Chen-Kai Qiao
Browse files
Options

Replace plot_xiaoqu.py

parent cb6ed6fe
Hide whitespace changes
Inline Side-by-side
Showing with 604 additions and 604 deletions
plot_xiaoqu.py
View file @ da4ca3ce
################################################# #################################################
## ##
## 这个文件对 AMSS-NCKU 数值相对论的结果进行画图 ## 这个文件对 AMSS-NCKU 数值相对论的结果进行画图
## 小曲 ## 小曲
## 2024/10/01 --- 2024/12/06 ## 2024/10/01 --- 2024/12/06
## ##
################################################# #################################################
import numpy ## 导入 numpy 包进行数组的操作 import numpy ## 导入 numpy 包进行数组的操作
import matplotlib.pyplot as plt ## 导入 matplotlib 包进行画图 import matplotlib.pyplot as plt ## 导入 matplotlib 包进行画图
from mpl_toolkits.mplot3d import Axes3D ## 画 3 维图需要 from mpl_toolkits.mplot3d import Axes3D ## 画 3 维图需要
import glob import glob
import os ## 导入 os 包进行系统操作 import os ## 导入 os 包进行系统操作
import plot_binary_data import plot_binary_data
import NR_xiaoqu_Input as input_data import AMSS_NCKU_Input as input_data
#################################################################################### ####################################################################################
## 该函数根据二进制数据画出所有二维图 ## 该函数根据二进制数据画出所有二维图
def generate_binary_data_plot(binary_outdir, figure_outdir): def generate_binary_data_plot(binary_outdir, figure_outdir):
# 生成若干文件夹存放图片 # 生成若干文件夹存放图片
surface_plot_outdir = os.path.join(figure_outdir, "surface plot") surface_plot_outdir = os.path.join(figure_outdir, "surface plot")
os.mkdir(surface_plot_outdir) os.mkdir(surface_plot_outdir)
density_plot_outdir = os.path.join(figure_outdir, "density plot") density_plot_outdir = os.path.join(figure_outdir, "density plot")
os.mkdir(density_plot_outdir) os.mkdir(density_plot_outdir)
contour_plot_outdir = os.path.join(figure_outdir, "contour plot") contour_plot_outdir = os.path.join(figure_outdir, "contour plot")
os.mkdir(contour_plot_outdir) os.mkdir(contour_plot_outdir)
print() print()
print(" 读取 AMSS-NCKU 程序的二进制数据 ") print(" 读取 AMSS-NCKU 程序的二进制数据 ")
print() print()
print(" 二进制数据列表 ") print(" 二进制数据列表 ")
## 设置对什么文件画图（这里设置对所有二进制文件画图） ## 设置对什么文件画图（这里设置对所有二进制文件画图）
globby = glob.glob( os.path.join(binary_outdir, '*.bin') ) globby = glob.glob( os.path.join(binary_outdir, '*.bin') )
file_list = [] file_list = []
for x in sorted(globby): for x in sorted(globby):
file_list.append(x) file_list.append(x)
print(x) print(x)
## 对列表中的所有文件画图 ## 对列表中的所有文件画图
for filename in file_list: for filename in file_list:
print(filename) print(filename)
plot_binary_data.plot_binary_data(filename, binary_outdir, figure_outdir) plot_binary_data.plot_binary_data(filename, binary_outdir, figure_outdir)
print() print()
print(" 二进制数据画图已完成 ") print(" 二进制数据画图已完成 ")
print() print()
return return
#################################################################################### ####################################################################################
#################################################################################### ####################################################################################
## 该函数对黑洞轨迹画图 ## 该函数对黑洞轨迹画图
def generate_puncture_orbit_plot(outdir, figure_outdir): def generate_puncture_orbit_plot(outdir, figure_outdir):
print() print()
print(" 对黑洞轨迹进行画图 ") print(" 对黑洞轨迹进行画图 ")
print() print()
# 打开文件路径 # 打开文件路径
file0 = os.path.join(outdir, "bssn_BH.dat") file0 = os.path.join(outdir, "bssn_BH.dat")
print(" 对应数据文件为 ", file0) print(" 对应数据文件为 ", file0)
# 读取整个文件数据，假设数据是以空格分隔的浮点数 # 读取整个文件数据，假设数据是以空格分隔的浮点数
data = numpy.loadtxt(file0) data = numpy.loadtxt(file0)
# print(data[:,0]) # print(data[:,0])
# print(data[:,2]) # print(data[:,2])
# -------------------------- # --------------------------
# 画出黑洞位移的轨迹图（XY图） # 画出黑洞位移的轨迹图（XY图）
plt.figure( figsize=(8,8) ) ## 这里 figsize 可以设定图形的大小 plt.figure( figsize=(8,8) ) ## 这里 figsize 可以设定图形的大小
plt.title( " Black Hole Trajectry ", fontsize=18 ) ## 这里 fontsize 可以设定文字大小 plt.title( " Black Hole Trajectry ", fontsize=18 ) ## 这里 fontsize 可以设定文字大小
for i in range(input_data.puncture_number): for i in range(input_data.puncture_number):
BH_x = data[:, 3*i+1] BH_x = data[:, 3*i+1]
BH_y = data[:, 3*i+2] BH_y = data[:, 3*i+2]
BH_z = data[:, 3*i+3] BH_z = data[:, 3*i+3]
if i==0: if i==0:
plt.plot( BH_x, BH_y, color='red', label="BH"+str(i+1), linewidth=2 ) plt.plot( BH_x, BH_y, color='red', label="BH"+str(i+1), linewidth=2 )
elif i==1: elif i==1:
plt.plot( BH_x, BH_y, color='green', label="BH"+str(i+1), linewidth=2 ) plt.plot( BH_x, BH_y, color='green', label="BH"+str(i+1), linewidth=2 )
elif i==2: elif i==2:
plt.plot( BH_x, BH_y, color='blue', label="BH"+str(i+1), linewidth=2 ) plt.plot( BH_x, BH_y, color='blue', label="BH"+str(i+1), linewidth=2 )
elif i==3: elif i==3:
plt.plot( BH_x, BH_y, color='gray', label="BH"+str(i+1), linewidth=2 ) plt.plot( BH_x, BH_y, color='gray', label="BH"+str(i+1), linewidth=2 )
plt.xlabel( "X [M]", fontsize=16 ) plt.xlabel( "X [M]", fontsize=16 )
plt.ylabel( "Y [M]", fontsize=16 ) plt.ylabel( "Y [M]", fontsize=16 )
plt.legend( loc='upper right' ) plt.legend( loc='upper right' )
# plt.show() # plt.show()
plt.savefig( os.path.join(figure_outdir, "BH_Trajectary_XY.pdf") ) plt.savefig( os.path.join(figure_outdir, "BH_Trajectary_XY.pdf") )
plt.close() plt.close()
# -------------------------- # --------------------------
# 画出黑洞位移的轨迹图（XZ图） # 画出黑洞位移的轨迹图（XZ图）
plt.figure( figsize=(8,8) ) ## 这里 figsize 可以设定图形的大小 plt.figure( figsize=(8,8) ) ## 这里 figsize 可以设定图形的大小
plt.title( " Black Hole Trajectry ", fontsize=18 ) ## 这里 fontsize 可以设定文字大小 plt.title( " Black Hole Trajectry ", fontsize=18 ) ## 这里 fontsize 可以设定文字大小
for i in range(input_data.puncture_number): for i in range(input_data.puncture_number):
BH_x = data[:, 3*i+1] BH_x = data[:, 3*i+1]
BH_y = data[:, 3*i+2] BH_y = data[:, 3*i+2]
BH_z = data[:, 3*i+3] BH_z = data[:, 3*i+3]
if i==0: if i==0:
plt.plot( BH_x, BH_z, color='red', label="BH"+str(i+1), linewidth=2 ) plt.plot( BH_x, BH_z, color='red', label="BH"+str(i+1), linewidth=2 )
elif i==1: elif i==1:
plt.plot( BH_x, BH_z, color='green', label="BH"+str(i+1), linewidth=2 ) plt.plot( BH_x, BH_z, color='green', label="BH"+str(i+1), linewidth=2 )
elif i==2: elif i==2:
plt.plot( BH_x, BH_z, color='blue', label="BH"+str(i+1), linewidth=2 ) plt.plot( BH_x, BH_z, color='blue', label="BH"+str(i+1), linewidth=2 )
elif i==3: elif i==3:
plt.plot( BH_x, BH_z, color='gray', label="BH"+str(i+1), linewidth=2 ) plt.plot( BH_x, BH_z, color='gray', label="BH"+str(i+1), linewidth=2 )
plt.xlabel( "X [M]", fontsize=16 ) plt.xlabel( "X [M]", fontsize=16 )
plt.ylabel( "Z [M]", fontsize=16 ) plt.ylabel( "Z [M]", fontsize=16 )
plt.legend( loc='upper right' ) plt.legend( loc='upper right' )
# plt.show() # plt.show()
plt.savefig( os.path.join(figure_outdir, "BH_Trajectary_XZ.pdf") ) plt.savefig( os.path.join(figure_outdir, "BH_Trajectary_XZ.pdf") )
plt.close() plt.close()
# -------------------------- # --------------------------
# 画出黑洞位移的轨迹图（YZ图） # 画出黑洞位移的轨迹图（YZ图）
plt.figure( figsize=(8,8) ) ## 这里 figsize 可以设定图形的大小 plt.figure( figsize=(8,8) ) ## 这里 figsize 可以设定图形的大小
plt.title( " Black Hole Trajectry ", fontsize=18 ) ## 这里 fontsize 可以设定文字大小 plt.title( " Black Hole Trajectry ", fontsize=18 ) ## 这里 fontsize 可以设定文字大小
for i in range(input_data.puncture_number): for i in range(input_data.puncture_number):
BH_x = data[:, 3*i+1] BH_x = data[:, 3*i+1]
BH_y = data[:, 3*i+2] BH_y = data[:, 3*i+2]
BH_z = data[:, 3*i+3] BH_z = data[:, 3*i+3]
if i==0: if i==0:
plt.plot( BH_y, BH_z, color='red', label="BH"+str(i+1), linewidth=2 ) plt.plot( BH_y, BH_z, color='red', label="BH"+str(i+1), linewidth=2 )
elif i==1: elif i==1:
plt.plot( BH_y, BH_z, color='green', label="BH"+str(i+1), linewidth=2 ) plt.plot( BH_y, BH_z, color='green', label="BH"+str(i+1), linewidth=2 )
elif i==2: elif i==2:
plt.plot( BH_y, BH_z, color='blue', label="BH"+str(i+1), linewidth=2 ) plt.plot( BH_y, BH_z, color='blue', label="BH"+str(i+1), linewidth=2 )
elif i==3: elif i==3:
plt.plot( BH_y, BH_z, color='gray', label="BH"+str(i+1), linewidth=2 ) plt.plot( BH_y, BH_z, color='gray', label="BH"+str(i+1), linewidth=2 )
plt.xlabel( "Y [M]", fontsize=16 ) plt.xlabel( "Y [M]", fontsize=16 )
plt.ylabel( "Z [M]", fontsize=16 ) plt.ylabel( "Z [M]", fontsize=16 )
plt.legend( loc='upper right' ) plt.legend( loc='upper right' )
# plt.show() # plt.show()
plt.savefig( os.path.join(figure_outdir, "BH_Trajectary_YZ.pdf") ) plt.savefig( os.path.join(figure_outdir, "BH_Trajectary_YZ.pdf") )
plt.close() plt.close()
# -------------------------- # --------------------------
# 得到黑洞 1 和黑洞 2 的坐标 # 得到黑洞 1 和黑洞 2 的坐标
BH_x1 = data[:, 1] BH_x1 = data[:, 1]
BH_y1 = data[:, 2] BH_y1 = data[:, 2]
BH_z1 = data[:, 3] BH_z1 = data[:, 3]
BH_x2 = data[:, 4] BH_x2 = data[:, 4]
BH_y2 = data[:, 5] BH_y2 = data[:, 5]
BH_z2 = data[:, 6] BH_z2 = data[:, 6]
# -------------------------- # --------------------------
# 画出黑洞位移的轨迹图（X2-X1 Y2-Y1） # 画出黑洞位移的轨迹图（X2-X1 Y2-Y1）
plt.figure( figsize=(8,8) ) plt.figure( figsize=(8,8) )
plt.title( " Black Hole Trajectry ", fontsize=18 ) plt.title( " Black Hole Trajectry ", fontsize=18 )
plt.plot( (BH_x2-BH_x1), (BH_y2-BH_y1), color='blue', linewidth=2 ) plt.plot( (BH_x2-BH_x1), (BH_y2-BH_y1), color='blue', linewidth=2 )
plt.xlabel( " $X_{2}$ - $X_{1}$ [M] ", fontsize=16 ) plt.xlabel( " $X_{2}$ - $X_{1}$ [M] ", fontsize=16 )
plt.ylabel( " $Y_{2}$ - $Y_{1}$ [M] ", fontsize=16 ) plt.ylabel( " $Y_{2}$ - $Y_{1}$ [M] ", fontsize=16 )
plt.legend( loc='upper right' ) plt.legend( loc='upper right' )
plt.savefig( os.path.join(figure_outdir, "BH_Trajectary_21_XY.pdf") ) plt.savefig( os.path.join(figure_outdir, "BH_Trajectary_21_XY.pdf") )
plt.close() plt.close()
# -------------------------- # --------------------------
# 画出黑洞位移的轨迹图（X2-X1 Z2-Z1） # 画出黑洞位移的轨迹图（X2-X1 Z2-Z1）
plt.figure( figsize=(8,8) ) plt.figure( figsize=(8,8) )
plt.title( " Black Hole Trajectry ", fontsize=18 ) plt.title( " Black Hole Trajectry ", fontsize=18 )
plt.plot( (BH_x2-BH_x1), (BH_z2-BH_z1), color='blue', linewidth=2 ) plt.plot( (BH_x2-BH_x1), (BH_z2-BH_z1), color='blue', linewidth=2 )
plt.xlabel( " $X_{2}$ - $X_{1}$ [M] ", fontsize=16 ) plt.xlabel( " $X_{2}$ - $X_{1}$ [M] ", fontsize=16 )
plt.ylabel( " $Z_{2}$ - $Z_{1}$ [M] ", fontsize=16 ) plt.ylabel( " $Z_{2}$ - $Z_{1}$ [M] ", fontsize=16 )
plt.legend( loc='upper right' ) plt.legend( loc='upper right' )
plt.savefig( os.path.join(figure_outdir, "BH_Trajectary_21_XZ.pdf") ) plt.savefig( os.path.join(figure_outdir, "BH_Trajectary_21_XZ.pdf") )
plt.close() plt.close()
# -------------------------- # --------------------------
# 画出黑洞位移的轨迹图（Y2-Y1 Z2-Z1） # 画出黑洞位移的轨迹图（Y2-Y1 Z2-Z1）
plt.figure( figsize=(8,8) ) plt.figure( figsize=(8,8) )
plt.title( " Black Hole Trajectry ", fontsize=18 ) plt.title( " Black Hole Trajectry ", fontsize=18 )
plt.plot( (BH_y2-BH_y1), (BH_z2-BH_z1), color='blue', linewidth=2 ) plt.plot( (BH_y2-BH_y1), (BH_z2-BH_z1), color='blue', linewidth=2 )
plt.xlabel( " $Y_{2}$ - $Y_{1}$ [M] ", fontsize=16 ) plt.xlabel( " $Y_{2}$ - $Y_{1}$ [M] ", fontsize=16 )
plt.ylabel( " $Z_{2}$ - $Z_{1}$ [M] ", fontsize=16 ) plt.ylabel( " $Z_{2}$ - $Z_{1}$ [M] ", fontsize=16 )
plt.legend( loc='upper right' ) plt.legend( loc='upper right' )
plt.savefig( os.path.join(figure_outdir, "BH_Trajectary_21_YZ.pdf") ) plt.savefig( os.path.join(figure_outdir, "BH_Trajectary_21_YZ.pdf") )
plt.close() plt.close()
# -------------------------- # --------------------------
# 报错 # 报错
# 这里 file0 只是个文件名，不涉及 file.open 操作 # 这里 file0 只是个文件名，不涉及 file.open 操作
# file0.close() # file0.close()
print() print()
print(" 对黑洞轨迹画图完成 ") print(" 对黑洞轨迹画图完成 ")
print() print()
return return
#################################################################################### ####################################################################################
#################################################################################### ####################################################################################
## 该函数对黑洞轨迹画 3 维图 ## 该函数对黑洞轨迹画 3 维图
def generate_puncture_orbit_plot3D(outdir, figure_outdir): def generate_puncture_orbit_plot3D(outdir, figure_outdir):
print() print()
print(" 对黑洞轨迹进行画 3 维图 ") print(" 对黑洞轨迹进行画 3 维图 ")
print() print()
# 打开文件路径 # 打开文件路径
file0 = os.path.join(outdir, "bssn_BH.dat") file0 = os.path.join(outdir, "bssn_BH.dat")
print(" 对应数据文件为 ", file0) print(" 对应数据文件为 ", file0)
# 读取整个文件数据，假设数据是以空格分隔的浮点数 # 读取整个文件数据，假设数据是以空格分隔的浮点数
data = numpy.loadtxt(file0) data = numpy.loadtxt(file0)
# 创建一个新的图 # 创建一个新的图
fig = plt.figure( figsize=(8,8) ) fig = plt.figure( figsize=(8,8) )
# 创建一个3D坐标轴 # 创建一个3D坐标轴
ax = fig.add_subplot(111, projection='3d') ax = fig.add_subplot(111, projection='3d')
# 添加标题 # 添加标题
ax.set_title( " Black Hole Trajectry ", fontsize=18 ) ax.set_title( " Black Hole Trajectry ", fontsize=18 )
for i in range(input_data.puncture_number): for i in range(input_data.puncture_number):
BH_x = data[:, 3*i+1] BH_x = data[:, 3*i+1]
BH_y = data[:, 3*i+2] BH_y = data[:, 3*i+2]
BH_z = data[:, 3*i+3] BH_z = data[:, 3*i+3]
if i==0: if i==0:
ax.plot( BH_x, BH_y, BH_z, color='red', label="BH"+str(i+1), linewidth=2 ) ax.plot( BH_x, BH_y, BH_z, color='red', label="BH"+str(i+1), linewidth=2 )
elif i==1: elif i==1:
ax.plot( BH_x, BH_y, BH_z, color='green', label="BH"+str(i+1), linewidth=2 ) ax.plot( BH_x, BH_y, BH_z, color='green', label="BH"+str(i+1), linewidth=2 )
elif i==2: elif i==2:
ax.plot( BH_x, BH_y, BH_z, color='blue', label="BH"+str(i+1), linewidth=2 ) ax.plot( BH_x, BH_y, BH_z, color='blue', label="BH"+str(i+1), linewidth=2 )
elif i==3: elif i==3:
ax.plot( BH_x, BH_y, BH_z, color='gray', label="BH"+str(i+1), linewidth=2 ) ax.plot( BH_x, BH_y, BH_z, color='gray', label="BH"+str(i+1), linewidth=2 )
# 设置坐标轴标签 # 设置坐标轴标签
ax.set_xlabel( "X [M]", fontsize=16 ) ax.set_xlabel( "X [M]", fontsize=16 )
ax.set_ylabel( "Y [M]", fontsize=16 ) ax.set_ylabel( "Y [M]", fontsize=16 )
ax.set_zlabel( "Z [M]", fontsize=16 ) ax.set_zlabel( "Z [M]", fontsize=16 )
plt.legend( loc='upper right' ) plt.legend( loc='upper right' )
# plt.show() # plt.show()
plt.savefig( os.path.join(figure_outdir, "BH_Trajectary_3D.pdf") ) plt.savefig( os.path.join(figure_outdir, "BH_Trajectary_3D.pdf") )
plt.close() plt.close()
return return
# 利用 Python 画 3 维图 # 利用 Python 画 3 维图
# import matplotlib.pyplot as plt # import matplotlib.pyplot as plt
# from mpl_toolkits.mplot3d import Axes3D # from mpl_toolkits.mplot3d import Axes3D
# 创建一个新的图 # 创建一个新的图
# fig = plt.figure() # fig = plt.figure()
# 创建一个3D坐标轴 # 创建一个3D坐标轴
# ax = fig.add_subplot(111, projection='3d') # ax = fig.add_subplot(111, projection='3d')
# 生成一些数据 # 生成一些数据
# x = [1, 2, 3, 4, 5] # x = [1, 2, 3, 4, 5]
# y = [5, 6, 2, 3, 13] # y = [5, 6, 2, 3, 13]
# z = [2, 3, 3, 3, 5] # z = [2, 3, 3, 3, 5]
# 在3D空间中绘制散点图 # 在3D空间中绘制散点图
# ax.scatter(x, y, z) # ax.scatter(x, y, z)
# 绘制三维曲线 # 绘制三维曲线
# ax.plot(x, y, z, label='3D curve') # ax.plot(x, y, z, label='3D curve')
# 设置坐标轴标签 # 设置坐标轴标签
# ax.set_xlabel('X Axis') # ax.set_xlabel('X Axis')
# ax.set_ylabel('Y Axis') # ax.set_ylabel('Y Axis')
# ax.set_zlabel('Z Axis') # ax.set_zlabel('Z Axis')
# 添加标题 # 添加标题
# ax.set_title('3D Curve Plot with Title') # ax.set_title('3D Curve Plot with Title')
#################################################################################### ####################################################################################
#################################################################################### ####################################################################################
## 该函数对引力波波形 Psi4 画图 ## 该函数对引力波波形 Psi4 画图
def generate_gravitational_waveform_plot(outdir, figure_outdir, detector_number_i): def generate_gravitational_waveform_plot(outdir, figure_outdir, detector_number_i):
print() print()
print(" 对引力波波形 Psi4 进行画图 ") print(" 对引力波波形 Psi4 进行画图 ")
print(" 对第 ", detector_number_i, " 个探测器半径数据进行画图 ") print(" 对第 ", detector_number_i, " 个探测器半径数据进行画图 ")
print() print()
# 打开文件路径 # 打开文件路径
file0 = os.path.join(outdir, "bssn_psi4.dat") file0 = os.path.join(outdir, "bssn_psi4.dat")
print(" 对应数据文件为 ", file0) print(" 对应数据文件为 ", file0)
# 读取整个文件数据，假设数据是以空格分隔的浮点数 # 读取整个文件数据，假设数据是以空格分隔的浮点数
data = numpy.loadtxt(file0) data = numpy.loadtxt(file0)
# 取出 phi4 文件中各列的数据 # 取出 phi4 文件中各列的数据
time = data[:,0] time = data[:,0]
psi4_l2m2m_real = data[:,1] psi4_l2m2m_real = data[:,1]
psi4_l2m2m_imaginary = data[:,2] psi4_l2m2m_imaginary = data[:,2]
psi4_l2m1m_real = data[:,3] psi4_l2m1m_real = data[:,3]
psi4_l2m1m_imaginary = data[:,4] psi4_l2m1m_imaginary = data[:,4]
psi4_l2m0_real = data[:,5] psi4_l2m0_real = data[:,5]
psi4_l2m0_imaginary = data[:,6] psi4_l2m0_imaginary = data[:,6]
psi4_l2m1_real = data[:,7] psi4_l2m1_real = data[:,7]
psi4_l2m1_imaginary = data[:,8] psi4_l2m1_imaginary = data[:,8]
psi4_l2m2_real = data[:,9] psi4_l2m2_real = data[:,9]
psi4_l2m2_imaginary = data[:,10] psi4_l2m2_imaginary = data[:,10]
# 报错 # 报错
# 这里 file0 只是个文件名，不涉及 file.open 操作 # 这里 file0 只是个文件名，不涉及 file.open 操作
# file0.close() # file0.close()
# python 中除法会返回浮点数，因此这里设置为整除 # python 中除法会返回浮点数，因此这里设置为整除
length = len(time) // input_data.Detector_Number length = len(time) // input_data.Detector_Number
time2 = numpy.zeros( (input_data.Detector_Number, length) ) time2 = numpy.zeros( (input_data.Detector_Number, length) )
psi4_l2m2m_real2 = numpy.zeros( (input_data.Detector_Number, length) ) psi4_l2m2m_real2 = numpy.zeros( (input_data.Detector_Number, length) )
psi4_l2m2m_imaginary2 = numpy.zeros( (input_data.Detector_Number, length) ) psi4_l2m2m_imaginary2 = numpy.zeros( (input_data.Detector_Number, length) )
psi4_l2m1m_real2 = numpy.zeros( (input_data.Detector_Number, length) ) psi4_l2m1m_real2 = numpy.zeros( (input_data.Detector_Number, length) )
psi4_l2m1m_imaginary2 = numpy.zeros( (input_data.Detector_Number, length) ) psi4_l2m1m_imaginary2 = numpy.zeros( (input_data.Detector_Number, length) )
psi4_l2m0_real2 = numpy.zeros( (input_data.Detector_Number, length) ) psi4_l2m0_real2 = numpy.zeros( (input_data.Detector_Number, length) )
psi4_l2m0_imaginary2 = numpy.zeros( (input_data.Detector_Number, length) ) psi4_l2m0_imaginary2 = numpy.zeros( (input_data.Detector_Number, length) )
psi4_l2m1_real2 = numpy.zeros( (input_data.Detector_Number, length) ) psi4_l2m1_real2 = numpy.zeros( (input_data.Detector_Number, length) )
psi4_l2m1_imaginary2 = numpy.zeros( (input_data.Detector_Number, length) ) psi4_l2m1_imaginary2 = numpy.zeros( (input_data.Detector_Number, length) )
psi4_l2m2_real2 = numpy.zeros( (input_data.Detector_Number, length) ) psi4_l2m2_real2 = numpy.zeros( (input_data.Detector_Number, length) )
psi4_l2m2_imaginary2 = numpy.zeros( (input_data.Detector_Number, length) ) psi4_l2m2_imaginary2 = numpy.zeros( (input_data.Detector_Number, length) )
# 将数据拆分为各探测器半径对应的数据 # 将数据拆分为各探测器半径对应的数据
for i in range(input_data.Detector_Number): for i in range(input_data.Detector_Number):
for j in range(length): for j in range(length):
time2[i,j] = time[ j*input_data.Detector_Number + i ] time2[i,j] = time[ j*input_data.Detector_Number + i ]
psi4_l2m2m_real2[i,j] = psi4_l2m2m_real[ j*input_data.Detector_Number + i ] psi4_l2m2m_real2[i,j] = psi4_l2m2m_real[ j*input_data.Detector_Number + i ]
psi4_l2m2m_imaginary2[i,j] = psi4_l2m2m_imaginary[ j*input_data.Detector_Number + i ] psi4_l2m2m_imaginary2[i,j] = psi4_l2m2m_imaginary[ j*input_data.Detector_Number + i ]
psi4_l2m1m_real2[i,j] = psi4_l2m1m_real[ j*input_data.Detector_Number + i ] psi4_l2m1m_real2[i,j] = psi4_l2m1m_real[ j*input_data.Detector_Number + i ]
psi4_l2m1m_imaginary2[i,j] = psi4_l2m1m_imaginary[ j*input_data.Detector_Number + i ] psi4_l2m1m_imaginary2[i,j] = psi4_l2m1m_imaginary[ j*input_data.Detector_Number + i ]
psi4_l2m0_real2[i,j] = psi4_l2m0_real[ j*input_data.Detector_Number + i ] psi4_l2m0_real2[i,j] = psi4_l2m0_real[ j*input_data.Detector_Number + i ]
psi4_l2m0_imaginary2[i,j] = psi4_l2m0_imaginary[ j*input_data.Detector_Number + i ] psi4_l2m0_imaginary2[i,j] = psi4_l2m0_imaginary[ j*input_data.Detector_Number + i ]
psi4_l2m1_real2[i,j] = psi4_l2m1_real[ j*input_data.Detector_Number + i ] psi4_l2m1_real2[i,j] = psi4_l2m1_real[ j*input_data.Detector_Number + i ]
psi4_l2m1_imaginary2[i,j] = psi4_l2m1_imaginary[ j*input_data.Detector_Number + i ] psi4_l2m1_imaginary2[i,j] = psi4_l2m1_imaginary[ j*input_data.Detector_Number + i ]
psi4_l2m2_real2[i,j] = psi4_l2m2_real[ j*input_data.Detector_Number + i ] psi4_l2m2_real2[i,j] = psi4_l2m2_real[ j*input_data.Detector_Number + i ]
psi4_l2m2_imaginary2[i,j] = psi4_l2m2_imaginary[ j*input_data.Detector_Number + i ] psi4_l2m2_imaginary2[i,j] = psi4_l2m2_imaginary[ j*input_data.Detector_Number + i ]
# 根据输入数据推算出探测器距离 # 根据输入数据推算出探测器距离
Detector_Interval = ( input_data.Detector_Rmax - input_data.Detector_Rmin ) / ( input_data.Detector_Number - 1 ) Detector_Interval = ( input_data.Detector_Rmax - input_data.Detector_Rmin ) / ( input_data.Detector_Number - 1 )
Detector_Distance_R = input_data.Detector_Rmax - Detector_Interval * detector_number_i Detector_Distance_R = input_data.Detector_Rmax - Detector_Interval * detector_number_i
plt.figure( figsize=(8,8) ) ## 这里 figsize 可以设定图形的大小 plt.figure( figsize=(8,8) ) ## 这里 figsize 可以设定图形的大小
plt.title( f" Gravitational Wave $\Psi_{4}$ Detector Distence = {Detector_Distance_R}", fontsize=18 ) ## 这里 fontsize 可以设定文字大小 plt.title( f" Gravitational Wave $\Psi_{4}$ Detector Distence = {Detector_Distance_R}", fontsize=18 ) ## 这里 fontsize 可以设定文字大小
plt.plot( time2[detector_number_i], psi4_l2m0_real2[detector_number_i], \ plt.plot( time2[detector_number_i], psi4_l2m0_real2[detector_number_i], \
color='red', label="l=2 m=0 real", linewidth=2 ) color='red', label="l=2 m=0 real", linewidth=2 )
plt.plot( time2[detector_number_i], psi4_l2m1_imaginary2[detector_number_i], \ plt.plot( time2[detector_number_i], psi4_l2m1_imaginary2[detector_number_i], \
color='orange', label="l=2 m=1 imaginary", linestyle='--', linewidth=2 ) color='orange', label="l=2 m=1 imaginary", linestyle='--', linewidth=2 )
plt.plot( time2[detector_number_i], psi4_l2m1_real2[detector_number_i], \ plt.plot( time2[detector_number_i], psi4_l2m1_real2[detector_number_i], \
color='green', label="l=2 m=1 real", linewidth=2 ) color='green', label="l=2 m=1 real", linewidth=2 )
plt.plot( time2[detector_number_i], psi4_l2m1_imaginary2[detector_number_i], \ plt.plot( time2[detector_number_i], psi4_l2m1_imaginary2[detector_number_i], \
color='cyan', label="l=2 m=1 imaginary", linestyle='--', linewidth=2 ) color='cyan', label="l=2 m=1 imaginary", linestyle='--', linewidth=2 )
plt.plot( time2[detector_number_i], psi4_l2m2_real2[detector_number_i], \ plt.plot( time2[detector_number_i], psi4_l2m2_real2[detector_number_i], \
color='black', label="l=2 m=2 real", linewidth=2 ) color='black', label="l=2 m=2 real", linewidth=2 )
plt.plot( time2[detector_number_i], psi4_l2m2_imaginary2[detector_number_i], \ plt.plot( time2[detector_number_i], psi4_l2m2_imaginary2[detector_number_i], \
color='gray', label="l=2 m=1 imaginary", linestyle='--', linewidth=2 ) color='gray', label="l=2 m=1 imaginary", linestyle='--', linewidth=2 )
plt.xlabel( "T [M]", fontsize=16 ) plt.xlabel( "T [M]", fontsize=16 )
plt.ylabel( "strength", fontsize=16 ) plt.ylabel( "strength", fontsize=16 )
plt.legend( loc='upper right' ) plt.legend( loc='upper right' )
plt.savefig( os.path.join(figure_outdir, "Gravitational_Psi4_Detector_" + str(detector_number_i) + ".pdf") ) plt.savefig( os.path.join(figure_outdir, "Gravitational_Psi4_Detector_" + str(detector_number_i) + ".pdf") )
print() print()
print(" 第 ", detector_number_i, " 个探测器半径数据画图完成 ") print(" 第 ", detector_number_i, " 个探测器半径数据画图完成 ")
print(" 对引力波波形 Psi4 的画图完成 ") print(" 对引力波波形 Psi4 的画图完成 ")
print() print()
return return
#################################################################################### ####################################################################################
#################################################################################### ####################################################################################
## 该函数对时空 ADM 质量画图 ## 该函数对时空 ADM 质量画图
def generate_ADMmass_plot(outdir, figure_outdir, detector_number_i): def generate_ADMmass_plot(outdir, figure_outdir, detector_number_i):
print() print()
print(" 对时空 ADM 质量进行画图 ") print(" 对时空 ADM 质量进行画图 ")
print(" 对第 ", detector_number_i, " 个探测器半径数据进行画图 ") print(" 对第 ", detector_number_i, " 个探测器半径数据进行画图 ")
print() print()
# 打开文件路径 # 打开文件路径
file0 = os.path.join(outdir, "bssn_ADMQs.dat") file0 = os.path.join(outdir, "bssn_ADMQs.dat")
print(" 对应数据文件为 ", file0) print(" 对应数据文件为 ", file0)
# 读取整个文件数据，假设数据是以空格分隔的浮点数 # 读取整个文件数据，假设数据是以空格分隔的浮点数
data = numpy.loadtxt(file0) data = numpy.loadtxt(file0)
# 取出 AMD 动量文件中各列的数据 # 取出 AMD 动量文件中各列的数据
time = data[:,0] time = data[:,0]
ADM_mass = data[:,1] ADM_mass = data[:,1]
ADM_Px = data[:,2] ADM_Px = data[:,2]
ADM_Py = data[:,3] ADM_Py = data[:,3]
ADM_Pz = data[:,4] ADM_Pz = data[:,4]
ADM_Jx = data[:,5] ADM_Jx = data[:,5]
ADM_Jy = data[:,6] ADM_Jy = data[:,6]
ADM_Jz = data[:,7] ADM_Jz = data[:,7]
# 报错 # 报错
# 这里 file0 只是个文件名，不涉及 file.open 操作 # 这里 file0 只是个文件名，不涉及 file.open 操作
# file0.close() # file0.close()
# python 中除法会返回浮点数，因此这里设置为整除 # python 中除法会返回浮点数，因此这里设置为整除
length = len(time) // input_data.Detector_Number length = len(time) // input_data.Detector_Number
''' '''
# 将数据拆分为各探测器半径对应的数据 # 将数据拆分为各探测器半径对应的数据
time2 = time.reshape( (input_data.Detector_Number, length) ) time2 = time.reshape( (input_data.Detector_Number, length) )
ADM_mass2 = ADM_mass.reshape( (input_data.Detector_Number, length) ) ADM_mass2 = ADM_mass.reshape( (input_data.Detector_Number, length) )
ADM_Px2 = ADM_Px.reshape( (input_data.Detector_Number, length) ) ADM_Px2 = ADM_Px.reshape( (input_data.Detector_Number, length) )
ADM_Py2 = ADM_Py.reshape( (input_data.Detector_Number, length) ) ADM_Py2 = ADM_Py.reshape( (input_data.Detector_Number, length) )
ADM_Pz2 = ADM_Pz.reshape( (input_data.Detector_Number, length) ) ADM_Pz2 = ADM_Pz.reshape( (input_data.Detector_Number, length) )
ADM_Jx2 = ADM_Jx.reshape( (input_data.Detector_Number, length) ) ADM_Jx2 = ADM_Jx.reshape( (input_data.Detector_Number, length) )
ADM_Jy2 = ADM_Jy.reshape( (input_data.Detector_Number, length) ) ADM_Jy2 = ADM_Jy.reshape( (input_data.Detector_Number, length) )
ADM_Jz2 = ADM_Jz.reshape( (input_data.Detector_Number, length) ) ADM_Jz2 = ADM_Jz.reshape( (input_data.Detector_Number, length) )
''' '''
# reshape 的行和列没有搞清楚，换成笨办法 # reshape 的行和列没有搞清楚，换成笨办法
time2 = numpy.zeros( (input_data.Detector_Number, length) ) time2 = numpy.zeros( (input_data.Detector_Number, length) )
ADM_mass2 = numpy.zeros( (input_data.Detector_Number, length) ) ADM_mass2 = numpy.zeros( (input_data.Detector_Number, length) )
ADM_Px2 = numpy.zeros( (input_data.Detector_Number, length) ) ADM_Px2 = numpy.zeros( (input_data.Detector_Number, length) )
ADM_Py2 = numpy.zeros( (input_data.Detector_Number, length) ) ADM_Py2 = numpy.zeros( (input_data.Detector_Number, length) )
ADM_Pz2 = numpy.zeros( (input_data.Detector_Number, length) ) ADM_Pz2 = numpy.zeros( (input_data.Detector_Number, length) )
ADM_Jx2 = numpy.zeros( (input_data.Detector_Number, length) ) ADM_Jx2 = numpy.zeros( (input_data.Detector_Number, length) )
ADM_Jy2 = numpy.zeros( (input_data.Detector_Number, length) ) ADM_Jy2 = numpy.zeros( (input_data.Detector_Number, length) )
ADM_Jz2 = numpy.zeros( (input_data.Detector_Number, length) ) ADM_Jz2 = numpy.zeros( (input_data.Detector_Number, length) )
# 将数据拆分为各探测器半径对应的数据 # 将数据拆分为各探测器半径对应的数据
for i in range(input_data.Detector_Number): for i in range(input_data.Detector_Number):
for j in range(length): for j in range(length):
time2[i,j] = time[ j*input_data.Detector_Number + i ] time2[i,j] = time[ j*input_data.Detector_Number + i ]
ADM_mass2[i,j] = ADM_mass[ j*input_data.Detector_Number + i ] ADM_mass2[i,j] = ADM_mass[ j*input_data.Detector_Number + i ]
ADM_Px2[i,j] = ADM_Px[ j*input_data.Detector_Number + i ] ADM_Px2[i,j] = ADM_Px[ j*input_data.Detector_Number + i ]
ADM_Py2[i,j] = ADM_Py[ j*input_data.Detector_Number + i ] ADM_Py2[i,j] = ADM_Py[ j*input_data.Detector_Number + i ]
ADM_Pz2[i,j] = ADM_Pz[ j*input_data.Detector_Number + i ] ADM_Pz2[i,j] = ADM_Pz[ j*input_data.Detector_Number + i ]
# 根据输入数据推算出探测器距离 # 根据输入数据推算出探测器距离
Detector_Interval = ( input_data.Detector_Rmax - input_data.Detector_Rmin ) / ( input_data.Detector_Number - 1 ) Detector_Interval = ( input_data.Detector_Rmax - input_data.Detector_Rmin ) / ( input_data.Detector_Number - 1 )
Detector_Distance_R = input_data.Detector_Rmax - Detector_Interval * detector_number_i Detector_Distance_R = input_data.Detector_Rmax - Detector_Interval * detector_number_i
# 画出最外 detector 半径的 AMD 动量 # 画出最外 detector 半径的 AMD 动量
plt.figure( figsize=(8,8) ) plt.figure( figsize=(8,8) )
plt.title(f" ADM Momentum Detector Distence = {Detector_Distance_R}", fontsize=18 ) plt.title(f" ADM Momentum Detector Distence = {Detector_Distance_R}", fontsize=18 )
plt.plot( time2[detector_number_i], ADM_mass2[detector_number_i], color='red', label="ADM Mass", linewidth=2 ) plt.plot( time2[detector_number_i], ADM_mass2[detector_number_i], color='red', label="ADM Mass", linewidth=2 )
plt.plot( time2[detector_number_i], ADM_Px2[detector_number_i], color='green', label="ADM Px", linewidth=2 ) plt.plot( time2[detector_number_i], ADM_Px2[detector_number_i], color='green', label="ADM Px", linewidth=2 )
plt.plot( time2[detector_number_i], ADM_Py2[detector_number_i], color='cyan', label="ADM Py", linewidth=2 ) plt.plot( time2[detector_number_i], ADM_Py2[detector_number_i], color='cyan', label="ADM Py", linewidth=2 )
plt.plot( time2[detector_number_i], ADM_Pz2[detector_number_i], color='blue', label="ADM Pz", linewidth=2 ) plt.plot( time2[detector_number_i], ADM_Pz2[detector_number_i], color='blue', label="ADM Pz", linewidth=2 )
plt.xlabel( "T [M]", fontsize=16 ) plt.xlabel( "T [M]", fontsize=16 )
plt.ylabel( "ADM Momentum [M]", fontsize=16 ) plt.ylabel( "ADM Momentum [M]", fontsize=16 )
plt.legend( loc='upper right' ) plt.legend( loc='upper right' )
plt.savefig( os.path.join(figure_outdir, "ADM_Mass_Dector_" + str(detector_number_i) + ".pdf") ) plt.savefig( os.path.join(figure_outdir, "ADM_Mass_Dector_" + str(detector_number_i) + ".pdf") )
print() print()
print(" 第 ", detector_number_i, " 个探测器半径数据画图完成 ") print(" 第 ", detector_number_i, " 个探测器半径数据画图完成 ")
print(" 对时空 ADM 质量的画图完成 ") print(" 对时空 ADM 质量的画图完成 ")
print() print()
return return
#################################################################################### ####################################################################################
#################################################################################### ####################################################################################
## 该函数对哈密顿约束违反性况画图 ## 该函数对哈密顿约束违反性况画图
def generate_constraint_check_plot(outdir, figure_outdir, input_level_number): def generate_constraint_check_plot(outdir, figure_outdir, input_level_number):
print() print()
print(" 对哈密顿约束违反性况进行画图 ") print(" 对哈密顿约束违反性况进行画图 ")
print(" 对第 ", input_level_number, " 层网格数据进行画图 ") print(" 对第 ", input_level_number, " 层网格数据进行画图 ")
print() print()
# 打开文件路径 # 打开文件路径
file0 = os.path.join(outdir, "bssn_constraint.dat") file0 = os.path.join(outdir, "bssn_constraint.dat")
print(" 对应数据文件为 ", file0) print(" 对应数据文件为 ", file0)
# 读取整个文件数据，假设数据是以空格分隔的浮点数 # 读取整个文件数据，假设数据是以空格分隔的浮点数
data = numpy.loadtxt(file0) data = numpy.loadtxt(file0)
# 取出约束数据文件中各列的数据 # 取出约束数据文件中各列的数据
time = data[:,0] time = data[:,0]
Constraint_H = data[:,1] Constraint_H = data[:,1]
Constraint_Px = data[:,2] Constraint_Px = data[:,2]
Constraint_Py = data[:,3] Constraint_Py = data[:,3]
Constraint_Pz = data[:,4] Constraint_Pz = data[:,4]
Constraint_Gx = data[:,5] Constraint_Gx = data[:,5]
Constraint_Gy = data[:,6] Constraint_Gy = data[:,6]
Constraint_Gz = data[:,7] Constraint_Gz = data[:,7]
# 报错 # 报错
# 这里 file0 只是个文件名，不涉及 file.open 操作 # 这里 file0 只是个文件名，不涉及 file.open 操作
# file0.close() # file0.close()
# 初始化各类数据 # 初始化各类数据
if (input_data.basic_grid_set == "Patch"): if (input_data.basic_grid_set == "Patch"):
level_number = input_level_number level_number = input_level_number
length0 = input_data.grid_level length0 = input_data.grid_level
# python 中除法会返回浮点数，因此这里设置为整除 # python 中除法会返回浮点数，因此这里设置为整除
length1 = len(time) // length0 length1 = len(time) // length0
elif (input_data.basic_grid_set == "Shell-Patch"): elif (input_data.basic_grid_set == "Shell-Patch"):
# 如果格点类型选择为 Shell-Patch，网格层数加 1 # 如果格点类型选择为 Shell-Patch，网格层数加 1
level_number = input_level_number + 1 level_number = input_level_number + 1
length0 = input_data.grid_level + 1 length0 = input_data.grid_level + 1
# python 中除法会返回浮点数，因此这里设置为整除 # python 中除法会返回浮点数，因此这里设置为整除
length1 = len(time) // length0 length1 = len(time) // length0
time2 = numpy.zeros( (length0, length1) ) time2 = numpy.zeros( (length0, length1) )
Constraint_H2 = numpy.zeros( (length0, length1) ) Constraint_H2 = numpy.zeros( (length0, length1) )
Constraint_Px2 = numpy.zeros( (length0, length1) ) Constraint_Px2 = numpy.zeros( (length0, length1) )
Constraint_Py2 = numpy.zeros( (length0, length1) ) Constraint_Py2 = numpy.zeros( (length0, length1) )
Constraint_Pz2 = numpy.zeros( (length0, length1) ) Constraint_Pz2 = numpy.zeros( (length0, length1) )
Constraint_Gx2 = numpy.zeros( (length0, length1) ) Constraint_Gx2 = numpy.zeros( (length0, length1) )
Constraint_Gy2 = numpy.zeros( (length0, length1) ) Constraint_Gy2 = numpy.zeros( (length0, length1) )
Constraint_Gz2 = numpy.zeros( (length0, length1) ) Constraint_Gz2 = numpy.zeros( (length0, length1) )
# 将数据拆分为各探测器半径对应的数据 # 将数据拆分为各探测器半径对应的数据
for i in range(length0): for i in range(length0):
for j in range(length1): for j in range(length1):
time2[i,j] = time[ j*length0 + i ] time2[i,j] = time[ j*length0 + i ]
Constraint_H2[i,j] = Constraint_H[ j*length0 + i ] Constraint_H2[i,j] = Constraint_H[ j*length0 + i ]
Constraint_Px2[i,j] = Constraint_Px[ j*length0 + i ] Constraint_Px2[i,j] = Constraint_Px[ j*length0 + i ]
Constraint_Py2[i,j] = Constraint_Py[ j*length0 + i ] Constraint_Py2[i,j] = Constraint_Py[ j*length0 + i ]
Constraint_Pz2[i,j] = Constraint_Pz[ j*length0 + i ] Constraint_Pz2[i,j] = Constraint_Pz[ j*length0 + i ]
# 画出最外 detector 半径的约束违反 # 画出最外 detector 半径的约束违反
plt.figure( figsize=(8,8) ) plt.figure( figsize=(8,8) )
plt.title( f" ADM Constraint Grid Level = {input_level_number}", fontsize=18 ) plt.title( f" ADM Constraint Grid Level = {input_level_number}", fontsize=18 )
plt.plot( time2[level_number], Constraint_H2[level_number], color='red', label="ADM Constraint H", linewidth=2 ) plt.plot( time2[level_number], Constraint_H2[level_number], color='red', label="ADM Constraint H", linewidth=2 )
plt.plot( time2[level_number], Constraint_Px2[level_number], color='green', label="ADM Constraint Px", linewidth=2 ) plt.plot( time2[level_number], Constraint_Px2[level_number], color='green', label="ADM Constraint Px", linewidth=2 )
plt.plot( time2[level_number], Constraint_Py2[level_number], color='cyan', label="ADM Constraint Py", linewidth=2 ) plt.plot( time2[level_number], Constraint_Py2[level_number], color='cyan', label="ADM Constraint Py", linewidth=2 )
plt.plot( time2[level_number], Constraint_Pz2[level_number], color='blue', label="ADM Constraint Pz", linewidth=2 ) plt.plot( time2[level_number], Constraint_Pz2[level_number], color='blue', label="ADM Constraint Pz", linewidth=2 )
plt.xlabel( "T [M]", fontsize=16 ) plt.xlabel( "T [M]", fontsize=16 )
plt.ylabel( "ADM Constraint", fontsize=16 ) plt.ylabel( "ADM Constraint", fontsize=16 )
plt.legend( loc='upper right' ) plt.legend( loc='upper right' )
plt.savefig( os.path.join(figure_outdir, "ADM_Constraint_Grid_Level_" + str(input_level_number) + ".pdf") ) plt.savefig( os.path.join(figure_outdir, "ADM_Constraint_Grid_Level_" + str(input_level_number) + ".pdf") )
print() print()
print(" 对哈密顿约束违反性况的画图完成 ") print(" 对哈密顿约束违反性况的画图完成 ")
print(" 第 ", input_level_number, " 层网格数据画图完成 ") print(" 第 ", input_level_number, " 层网格数据画图完成 ")
print() print()
return return
#################################################################################### ####################################################################################
#################################################################################### ####################################################################################
# 单独使用的例子 # 单独使用的例子
# outdir = "D:\AMSS-NCKU\AMSS-NCKU-Data\AMSSNCKU-xiaoqu-test5" # outdir = "D:\AMSS-NCKU\AMSS-NCKU-Data\AMSSNCKU-xiaoqu-test5"
# generate_puncture_orbit_plot(outdir, outdir) # generate_puncture_orbit_plot(outdir, outdir)
#################################################################################### ####################################################################################
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment