"""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 = relcfile+'.9' try: f = open(relcfile) lines = f.readlines() except: print('Error: Could not open %s' % relcfile) return() idep,itnu,delt,delp,delm = [],[],[],[],[] i = 3 while i < len(lines): xlin = lines[i] xlin = xlin.replace('D','E',5) xlin = xlin.strip() val = xlin.split() itnu.append(float(val[0])) idep.append(float(val[1])) delt.append(float(val[2])) delp.append(float(val[4])) delm.append(float(val[6])) i += 1 niter = int(itnu[-1]) nd = int(idep[0]) xtit1 = 'depth index' xtit2 = 'iteration' ytit1 = 'rel.change' ytit2 = 'log(rel.change)' ftit1 = 'temperature' ftit2 = 'max.population' ftit3 = 'max in state vector' if max(delt) == 0: delt=delp ftit1 = ftit2 print(delt[0],delt[-1]) fig, axes = plt.subplots(2,3) axes[0, 0].set_xlabel(xtit1) axes[0, 0].set_ylabel(ytit1) axes[0, 0].set_title(ftit1) axes[0, 1].set_xlabel(xtit1) axes[0, 1].set_ylabel(ytit2) # axes[0, 1].set_title(ftit1) axes[1, 0].set_xlabel(xtit1) axes[1, 0].set_ylabel(ytit1) axes[1, 0].set_title(ftit3) axes[1, 1].set_xlabel(xtit1) axes[1, 1].set_ylabel(ytit2) axes[1, 1].set_title(ftit3) i=0 tmaxt, tmaxm, tmaxi = [], [], [] while i < niter: i0, i1 = i*nd, (i+1)*nd-1 axes[0, 0].plot(idep[i0:i1],delt[i0:i1]) axes[0, 1].plot(idep[i0:i1],np.log10(np.abs(delt[i0:i1]))) axes[1, 0].plot(idep[i0:i1],delm[i0:i1]) axes[1, 1].plot(idep[i0:i1],np.log10(np.abs(delm[i0:i1]))) tmaxt.append(max(np.log10(np.abs(delt[i0:i1])))) tmaxm.append(max(np.log10(np.abs(delm[i0:i1])))) i+=1 tmaxi.append(i) axes[0, 2].set_xlabel(xtit2) axes[0, 2].set_ylabel('max(rel.changei)') axes[0, 2].set_title(ftit1) # axes[0, 2].plot(tmaxi,tmaxt,'ko',ls=':',ms=10) axes[0, 2].plot(tmaxi,tmaxt,'ko',ls='--') axes[1, 2].set_xlabel(xtit2) axes[1, 2].set_ylabel('max(rel.change') axes[1, 2].set_title(ftit3) # axes[1, 2].plot(tmaxi,tmaxm,'ko',ls=':',ms=10) axes[1, 2].plot(tmaxi,tmaxm,'ko',ls='--') # overall title # contains the cores name of the model, and also shows the total # execution time used by Tlusty to generate the model form file # [core name].69 (if it exists) timfile = relcfile[0:-2]+'.69' try: ft = open(timfile) linti = ft.readlines() titlt = 'time = '+ linti[-1].split()[2] +' s' except: print('Error: Could not open %s' % relcfile) titlt=' ' titl=relcfile + ': ' + titlt fig.suptitle(titl) fig.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)