"""Utility programs for examining output from Tlusty At the moment, it contains the following routines: 1) pconv - for ploting the convergence log contasined in the output file from Tlusty (fort.9); 2) pmodel - extracting values of column mass and a selected state parameter from the Tlusty output file fort.7m (coindensed model) or fort.12 (b-factors for NLTE models); 3) pmodels - similar to pmodel, but with an option to plot selected state parameter for one or more Tlusty models, and a possibility of plotting differences of between individual models 4) pmods - a plot of a range pf state parameters for one or more Tlusty models Examples: -------- To plot a convergence log of the last computed model: >>> import matplotlib.pyplot as plt >>> import tlusty as tl >>> plt.ion() >>> tl.pconv() a convergence log for a previously computed model, say, hhe35lt: >>> tl.pconv('hhe35lt') To extract log(column mass) and temperature of the same model >>> (m35,t35) = tl.pmodel('hhe35lt') which can then be plotted or used for some other purpose. To plot the temperature as a function of log(m) for the last computed model, do >>> tl.pmodels() and to plot the temperature for three models, hhe35lt','hhe35nc, and hhe35nl, and then differences from the first model, do >>> tl.pmodels(['hhe35lt','hhe35nc','hhe35nl']) >>> tl.pmodels(['hhe35lt','hhe35nc','hhe35nl'],diff=True) To plot level populations for the first five levels of hydrogen for the model hhe35nl, do >>> tl.pmods('hhe35nl',[3,4,5,6,7]) """ # Tlusty输出分析工具程序 # 当前包含以下功能： # 1) pconv - 绘制Tlusty输出文件（fort.9）的收敛日志 # 2) pmodel - 从Tlusty输出文件（fort.7m或fort.12）提取列质量及指定状态参数 # 3) pmodels - 类似pmodel，支持多模型对比及差异分析 # 4) pmods - 绘制多个Tlusty模型的多状态参数 # 5) pflux - 绘制辐射通量或通量差异 import matplotlib.pyplot as plt # 导入绘图库 import numpy as np # 导入数值计算库 import sys # 系统相关功能 ################################################### def pconv(relcfile='fort'): """" 绘制Tlusty输出文件（单位9）的收敛日志 类似旧IDL程序PCONV.PRO，生成2x3子图： 第一行显示温度（或电离态）的相对变化 第二行显示所有状态参数的最大相对变化 左列显示各深度的相对变化值 中列显示相对变化的对数 右列显示每迭代步的最大变化值 右侧子图最能反映整体收敛情况 标题会显示运行时间（若存在fort.69文件） 参数： relcfile: str, 可选 - 模型核心文件名（自动添加.9后缀） 默认值：'fort'（对应fort.9） 返回：无 """ relcfile += '.9' # 补充文件后缀 try: f = open(relcfile) # 打开文件 lines = f.readlines() # 读取所有行 except: print(f'错误：无法打开{relcfile}') # 异常处理 return() idep,itnu,delt,delp,delm = [],[],[],[],[] # 初始化变量 i = 3 # 从第4行开始读取数据 while i < len(lines): # 循环处理每一行 xlin = lines[i].replace('D','E') # 替换指数符号 val = list(map(float, xlin.split())) # 解析数值 itnu.append(val[0]) # 迭代次数 idep.append(val[1]) # 深度索引 delt.append(val[2]) # 温度变化 delp.append(val[4]) # 电离变化 delm.append(val[6]) # 状态向量最大变化 i += 1 # 下一行 niter = int(itnu[-1]) # 总迭代次数 nd = int(idep[0]) # 深度层数 # 设置坐标轴标题 xtit1 = '深度索引' # x轴标题（左列） xtit2 = '迭代次数' # x轴标题（右列） ytit1 = '相对变化值' # y轴标题（左列） ytit2 = '相对变化对数' # y轴标题（中列） ftit1 = '温度' # 第一行标题 ftit2 = '最大电离变化' # 第二行左列标题 ftit3 = '状态向量最大变化' # 第二行右列标题 if max(delt) == 0: # 若温度未变化则使用电离变化 delt = delp ftit1 = ftit2 fig, axes = plt.subplots(2,3) # 创建2x3子图 # 设置各子图标题和坐标轴 axes[0,0].set(xlabel=xtit1, ylabel=ytit1, title=ftit1) axes[0,1].set(xlabel=xtit1, ylabel=ytit2) axes[1,0].set(xlabel=xtit1, ylabel=ytit1, title=ftit3) axes[1,1].set(xlabel=xtit1, ylabel=ytit2) i = 0 tmaxt, tmaxm, tmaxi = [], [], [] while i < niter: # 循环绘制各迭代步数据 i0, i1 = i*nd, (i+1)*nd # 绘制温度变化 axes[0,0].plot(idep[i0:i1], delt[i0:i1]) # 左上 axes[0,1].plot(idep[i0:i1], np.log10(abs(delt[i0:i1]))) # 中上 # 绘制状态向量变化 axes[1,0].plot(idep[i0:i1], delm[i0:i1]) # 左下 axes[1,1].plot(idep[i0:i1], np.log10(abs(delm[i0:i1]))) # 中下 # 记录最大值 tmaxt.append(max(np.log10(abs(delt[i0:i1])))) tmaxm.append(max(np.log10(abs(delm[i0:i1])))) tmaxi.append(i+1) # 迭代次数 i += 1 # 绘制右侧最大值变化曲线 axes[0,2].set(xlabel=xtit2, ylabel='最大相对变化', title=ftit1) axes[0,2].plot(tmaxi, tmaxt, 'ko--') # 右上 axes[1,2].set(xlabel=xtit2, ylabel='最大相对变化', title=ftit3) axes[1,2].plot(tmaxi, tmaxm, 'ko--') # 右下 # 添加总标题（包含运行时间） timfile = relcfile[:-2]+'.69' try: with open(timfile) as ft: linti = ft.readlines() titlt = f"耗时：{linti[-1].split()[2]} 秒" except: titlt = '' fig.suptitle(f"{relcfile}：{titlt}") plt.tight_layout() plt.show() # Next lines permit one to run the routine from the command line #if __name__ == "__main__": # import sys # if len(sys.argv)==4: # doit(sys.argv[1],sys.argv[2],sys.argv[3]) # else: # print (__doc__) # return() #################################################### def pmodel(modfile='fort', ind=0): """ 从Tlusty的压缩模型大气文件（单元7）中提取列质量和一个状态参数，也可从类似输出文件（单元12）提取b因子。 功能类似旧IDL程序PMODEL.PRO（或MODELS11.PRO）。 参数： ----------- modfile: str, 可选 模型大气文件的核心名称。若名称不含小数点"."，则自动添加".7" 默认值：'fort'，即文件名默认为'fort.7' ind: int, 可选 - 状态参数的索引 默认0，即温度 状态参数索引说明： 0 - 温度（K） 1 - 电子密度（cm⁻³） 2 - 质量密度（g/cm³） 3 - 总粒子数密度（仅含分子模型） 后续索引 - 能级布居数 返回： ----------- dm: array 列质量的对数（log10） par: array 所选状态参数的值（温度为原始值，其他为对数） """ # 判断文件名格式并补充后缀 if modfile.count('.') == 0: modfile += '.7' # 补充默认后缀".7" # 尝试打开文件 try: f = open(modfile) lines = f.readlines() # 读取所有行 except: print(f'错误：无法打开 {modfile}') # 异常处理 return () # 读取输入参数的深度数量和参数数量 x = lines[0] # 第一行数据 nd = int(x.split()[0]) # 总深度层数 npa = abs(int(x.split()[1])) # 状态参数数量 # 设置每行读取的深度数量（根据文件类型） ndrow = 6 if modfile.count('.12') > 0: ndrow = 8 # 若为b因子文件则调整 # 计算需要读取的行数 nr = nd // ndrow if nd % ndrow != 0: nr += 1 # 处理余数 # 读取深度数据 dm = [] # 初始化列质量列表 i = 0 while i < nr: x = lines[i+1] # 从第二行开始读取深度数据 x = x.replace('D','E',6) # 替换指数符号 val = x.split() # 分割数值 j = 0 while j < ndrow: k = ndrow * i + j if k < nd: dm.append(float(val[j])) # 存储深度值 j += 1 i += 1 # 读取状态参数数据 nrx = npa // 5 # 每行参数数量 if npa % 5 != 0: nrx += 1 # 处理余数 par = [] # 初始化状态参数列表 i = nr + 1 # 从深度数据后开始读取参数 id = 0 while id < nd: xpar = [] # 临时参数列表 j = 0 while j < nrx: x = lines[i] # 读取当前行 x = x.replace('D','E',6) # 替换指数符号 val = x.split() # 分割数值 k = 0 while k < 5: if 5*j + k < npa: xpar.append(float(val[k])) # 存储参数值 k += 1 j += 1 i += 1 par.append(xpar[ind]) # 保存指定索引的参数值 id += 1 # 数据处理与绘图 dm = np.log10(dm) # 转换为列质量对数 par = np.array(par) # 转换为数组 if ind > 0: par = np.log10(par) # 非温度参数取对数 # 创建绘图 fig, ax = plt.subplots() ax.set_xlabel('log mass') # X轴标签 ax.set_title(modfile) # 标题 # 根据参数类型设置Y轴标签和绘图 if ind == 0: ax.set_ylabel('temperature') # 温度标签 ax.plot(dm, par) # 绘制温度曲线 else: ax.set_ylabel('log n') # 密度/布居数标签 ax.plot(dm, np.log10(par)) # 绘制对数曲线 plt.show() # 显示图形 # return (dm, par) # 返回数据数组 #################################################### def pmodels(modfiles=['fort'], ind=0, diff=False): """ 类似pmodel函数，但支持多模型对比和差异分析 可绘制一个或多个Tlusty模型的某个状态参数，并可显示与首个模型的差异 参数： ----------- modfiles: list of str, 可选 模型文件的核心名称列表。若名称不含小数点"."，则自动添加".7" 默认值：['fort']（对应fort.7） ind: int, 可选 - 状态参数的索引 默认0（温度） diff: bool, 可选 若为True，绘制指定状态参数与首个模型的差异 默认False（绘制实际参数值） 状态参数索引说明： 0 - 温度（K） 1 - 电子密度（cm⁻³） 2 - 质量密度（g/cm³） 3 - 总粒子数密度（仅含分子模型） 后续索引 - 能级布居数 返回： ----------- 无（直接绘图） """ imod = 0 # 初始化模型计数器 for modfile in modfiles: # 遍历所有模型文件 if modfile.count('.') == 0: # 自动补充后缀".7" modfile += '.7' # 尝试打开文件 try: f = open(modfile) lines = f.readlines() except: print(f'错误：无法打开 {modfile}') return() # 读取输入参数的深度数量和参数数量 x = lines[0] # 第一行数据 nd = int(x.split()[0]) # 总深度层数 npa = abs(int(x.split()[1])) # 状态参数数量 # 设置每行读取的深度数量（根据文件类型） ndrow = 6 if modfile.count('.12') > 0: ndrow = 8 # 若为b因子文件则调整 # 计算需要读取的行数 nr = nd // ndrow if nd % ndrow != 0: nr += 1 # 处理余数 # 读取深度数据 dm = [] # 初始化列质量列表 i = 0 while i < nr: x = lines[i+1] # 从第二行开始读取深度数据 x = x.replace('D','E',6) # 替换指数符号 val = x.split() # 分割数值 j = 0 while j < ndrow: k = ndrow * i + j if k < nd: dm.append(float(val[j])) # 存储深度值 j += 1 i += 1 # 读取状态参数数据 nrx = npa // 5 # 每行参数数量 if npa % 5 != 0: nrx += 1 # 处理余数 par = [] # 初始化状态参数列表 i = nr + 1 # 从深度数据后开始读取参数 id = 0 while id < nd: xpar = [] # 临时参数列表 j = 0 while j < nrx: x = lines[i] # 读取当前行 x = x.replace('D','E',6) # 替换指数符号 val = x.split() # 分割数值 k = 0 while k < 5: if 5*j + k < npa: xpar.append(float(val[k])) # 存储参数值 k += 1 j += 1 i += 1 par.append(xpar[ind]) # 保存指定索引的参数值 id += 1 # 数据处理与绘图 dm = np.log10(dm) # 转换为列质量对数 if ind != 0: par = np.log10(par) # 非温度参数取对数 # 初始化绘图（首个模型） if imod == 0: fig, ax = plt.subplots() # 创建图形 ax.set_xlabel('log mass') # X轴标签 titl = modfile[:-2] # 标题初始化 dm0 = dm.copy() # 存储首个模型的深度数据 par0 = par.copy() # 存储首个模型的参数值 if diff: # 差异模式 par = np.array(par) - np.array(par0) # 计算差异 else: # 后续模型 titl += f', {modfile[:-2]}' # 更新标题 if diff: # 差异模式处理 # 线性插值到首个模型的深度网格 pari = np.interp(dm0, dm, par) par = pari - par0 # 计算与首个模型的差异 # 设置Y轴标签 if ind == 0: ax.set_ylabel('temperature') # 温度标签 else: ax.set_ylabel('log n') # 密度/布居数标签 # 绘制数据 if not diff: ax.plot(dm, par) # 绘制实际参数曲线 else: ax.plot(dm0, par) # 绘制差异曲线 imod += 1 # 增加模型计数器 ax.set_title(titl) # 设置最终标题 plt.show() # 显示图形 ##################################################################### def pmods(modfiles, ind): """ 类似pmodels函数，但支持同时绘制多个状态参数 可绘制一个或多个Tlusty模型的多个状态参数曲线 参数： ----------- modfiles: list of str 模型文件的核心名称列表。若名称不含小数点"."，则自动添加".7" ind: list of int 需要绘制的状态参数索引列表 状态参数索引说明： 0 - 温度（K） 1 - 电子密度（cm⁻³） 2 - 质量密度（g/cm³） 3 - 总粒子数密度（仅含分子模型） 后续索引 - 能级布居数 返回： ----------- 无（直接绘图） """ imod = 0 # 初始化模型计数器 for modfile in modfiles: # 遍历所有模型文件 if modfile.count('.') == 0: # 自动补充后缀".7" modfile += '.7' # 尝试打开文件 try: f = open(modfile) lines = f.readlines() except: print(f'错误：无法打开 {modfile}') return() # 读取输入参数的深度数量和参数数量 x = lines[0] # 第一行数据 nd = int(x.split()[0]) # 总深度层数 npa = abs(int(x.split()[1])) # 状态参数数量 # 设置每行读取的深度数量（根据文件类型） ndrow = 6 if modfile.count('.12') > 0: ndrow = 8 # 若为b因子文件则调整 # 计算需要读取的行数 nr = nd // ndrow if nd % ndrow != 0: nr += 1 # 处理余数 # 读取深度数据 dm = [] # 初始化列质量列表 i = 0 while i < nr: x = lines[i+1] # 从第二行开始读取深度数据 x = x.replace('D','E',6) # 替换指数符号 val = x.split() # 分割数值 j = 0 while j < ndrow: k = ndrow * i + j if k < nd: dm.append(float(val[j])) # 存储深度值 j += 1 i += 1 # 读取所有状态参数数据 nrx = npa // 5 # 每行参数数量 if npa % 5 != 0: nrx += 1 # 处理余数 totpar = [] # 存储所有参数的二维数组 i = nr + 1 # 从深度数据后开始读取参数 id = 0 while id < nd: xpar = [] # 临时参数列表 j = 0 while j < nrx: x = lines[i] # 读取当前行 x = x.replace('D','E',6) # 替换指数符号 val = x.split() # 分割数值 k = 0 while k < 5: if 5*j + k < npa: xpar.append(float(val[k])) # 存储参数值 k += 1 j += 1 i += 1 totpar.append(xpar) # 存储所有参数的当前层数据 id += 1 # 数据处理 pararr = np.log10(np.array(totpar)) # 转换为对数形式 dm = np.log10(dm) # 转换为列质量对数 # 初始化绘图（首个模型） if imod == 0: fig, ax = plt.subplots() # 创建图形 ax.set_xlabel('log mass') # X轴标签 ax.set_ylabel('log n') # Y轴标签 titl = modfile[:-2] # 标题初始化 # 绘制指定状态参数 for indx in ind: # 遍历所有需要绘制的参数索引 ax.plot(dm, pararr[:, indx]) # 绘制对应参数曲线 # 更新标题 if imod != 0: titl += f', {modfile[:-2]}' # 添加模型名称 imod += 1 # 增加模型计数器 ax.set_title(titl) # 设置最终标题 plt.show() # 显示图形 def pflux(modfiles=['fort'], rel=False, ilogx=False): """ 绘制Tlusty输出文件（单元14）的出射通量，或多个模型的通量相对差异 也可用于绘制Synspec生成的合成光谱 参数： ----------- modfiles: list of str, 可选 模型文件的核心名称列表。若名称不含小数点"."，则自动添加".14" 默认值：['fort']（对应fort.14） rel: bool, 可选 若为True，绘制指定模型的通量相对于首个模型的百分比差异： rel(i) = (flux(i)-flux(0))/flux(0)*100% 默认False（绘制实际通量值） ilogx: bool, 可选 若为True，设置X轴（波长）为对数坐标 默认False（线性坐标） 返回： ----------- 无（直接绘图） """ imod = 0 # 初始化模型计数器 for modfile in modfiles: # 遍历所有模型文件 if modfile.count('.') == 0: # 自动补充后缀".14" modfile += '.14' # 读取波长和通量数据 try: wl, fl = np.loadtxt(modfile, unpack=True) # 加载文件数据 except: print(f'错误：无法打开 {modfile}') return() # 初始化绘图（首个模型） if imod == 0: fig, ax = plt.subplots() # 创建图形 ax.set_xlabel('波长') # X轴标签 titl = modfile[:-3] # 标题初始化 wl0 = wl.copy() # 存储首个模型的波长数据 fl0 = fl.copy() # 存储首个模型的通量数据 if rel: # 相对差异模式 fl = (fl - fl0)/fl0 * 100 # 计算初始差异 else: # 后续模型 titl += f', {modfile[:-3]}' # 更新标题 if rel: # 相对差异模式处理 # 线性插值到首个模型的波长网格 fli = np.interp(wl0, wl, fl) fl = (fli - fl0)/fl0 * 100 # 计算相对差异 # 绘制数据 if not rel: if ilogx: ax.semilogx(wl, fl) # 对数坐标X轴 else: ax.plot(wl, fl) # 线性坐标X轴 else: if ilogx: ax.semilogx(wl0, fl) # 对数坐标X轴 + 相对差异 else: ax.plot(wl0, fl) # 线性坐标X轴 + 相对差异 imod += 1 # 增加模型计数器 ax.set_title(titl) # 设置最终标题 plt.show() # 显示图形 if __name__ == "__main__": file = sys.argv[1] pconv(file)