#=================================================================================================== # for script description run : # python3 make_disnht_custom.py -h #=================================================================================================== import os import numpy as np import argparse import matplotlib.pyplot as plt #-------------------------------------------------------------------------------------------- # FUNCTIONS def rebin ( ene, spec, xLowLim, xUppLim, Nbins ) : binAx = np.logspace( np.log10(xLowLim), np.log10(xUppLim), Nbins ) rebinned = np.array([ np.nanmean( spec[ (ene > binAx[i]) * (ene < binAx[i+1]) ] ) for i in range( len(binAx)-1 ) ]) return rebinned #-------------------------------------------------------------------------------------------- # MAIN def main() : parser = argparse.ArgumentParser( description = 'This program computes the proper absorption spectrum for a user defined distribution of column densities. The input is a file including the arrray of column densitiey values. Other optional inputs are a cross-section .dat file (--cross_section) including the 2-d array [energy, cross-section] and boolean flags for input description and output rebin, plot or save. The output table format is a .dat tab (--save) and/or a etable readable by XSPEC (--xspec); NOTE: Rebinning of the spectrum (--rebin) is set by default as 500 log-spaced steps (--Nbins) from 0.1 (--BinLow_keV) to 10 (--BinUpp_keV) keV.' ) parser.add_argument("input_NH", type = str, help = 'str; .dat file name of N_H distribution, example : numpy.savetxt( "input_NH_example.dat", numpy.random.lognormal( 2.3 * 21, 0.3, 100 ) )' ) parser.add_argument( "-o", "--output", type = str, default = 'disnht.dat', help = 'str; output file path for output absorption tpectra; default = "./disnht.{dat,fits}"' ) parser.add_argument( "-c", "--cross_section", type = str, default = 'cross_angr.dat', help = 'str; file name containing absorption cross section as a function of energy, shape = ( 2, N_Ebins ) ); choose within { cross_aneb.dat cross_angr.dat cross_aspl.dat cross_feld.dat cross_grsa.dat cross_lodd.dat cross_wilm.dat } or set a custom path' ) parser.add_argument( "-l", "--logarithmic", type = bool, nargs = '?', default = False, const = True, help = 'bool; set flag if column densities in .dat are log_10 values' ) parser.add_argument( "-r", "--rebin", type = bool, nargs = '?', default = False, const = True, help = 'bool; set flag to decrease spectrum resolution in output file' ) parser.add_argument("--BinLow_keV", type = float, default = 0.1, help = 'float; lower energy bound [keV] for output spectrum; default = 0.1; setting a different value overrides --rebin as True' ) parser.add_argument("--BinUpp_keV", type = float, default = 10., help = 'float; upper energy bound [keV] for output spectrum; default = 10; setting a different value overrides --rebin as True' ) parser.add_argument("--Nbins", type = int, default = 500, help = 'int; number of output energy bins; default = 500; setting a different value overrides --rebin as True' ) parser.add_argument( "-p", "--plot", type = bool, nargs = '?', default = False, const = True, help = 'bool; set flag to plot the output absorption spectrum' ) parser.add_argument( "-s", "--save", type = bool, nargs = '?', default = False, const = True, help = 'bool; set flag to save spectrum in output file "disnht_.dat"' ) parser.add_argument( "-x", "--xspec", type = bool, nargs = '?', default = False, const = True, help = 'bool; set flag to convert the output into XSPEC etable model "disnht_.fits"; overrides --save as True' ) parser.add_argument("-v", "--verbose", help = "bool, increase output verbosity", action = "store_true" ) args = parser.parse_args() if args.BinLow_keV != 0.1 or args.BinUpp_keV != 10. or args.Nbins != 500 : args.rebin = True cross = np.loadtxt( args.cross_section ) if args.logarithmic : NH_distr = np.loadtxt( args.input_NH ) else : NH_distr = np.log10( np.loadtxt( args.input_NH ) ) NH_size = len( NH_distr ) if args.verbose : print( 'Reading cross section from ', args.cross_section ) E_keV = cross[0] # compute absorbed spectra with the logNH in the distribution if not args.rebin : length = len( E_keV ) else : length = args.Nbins - 1 if args.verbose : print( 'Rebinning spectra with ', args.Nbins, ' bins (log-spaced) from ', args.BinLow_keV, ' to ', args.BinUpp_keV, ' keV' ) SpcAbs = np.zeros([ length, NH_size ]) for n in range( NH_size ) : # compute single NH absorbed spectrum SpcAbs_n = np.exp( -cross[1] * 10 ** NH_distr[n] ) / NH_size if not args.rebin : SpcAbs[ :, n ] = SpcAbs_n else : SpcAbs[ :, n ] = rebin( E_keV, SpcAbs_n, args.BinLow_keV, args.BinUpp_keV, args.Nbins ) # compute the spectrum absorbed by (i.e. average) of the NH distribution SpcNHavg = np.exp( -cross[1] * 10 ** np.mean( NH_distr ) ) # rebin before printing output and plot if args.rebin : SpcNHavg = rebin( E_keV, SpcNHavg, args.BinLow_keV, args.BinUpp_keV, args.Nbins ) E_keV = rebin( E_keV, E_keV, args.BinLow_keV, args.BinUpp_keV, args.Nbins ) # sum the single spectra to get total emission SpcTot = np.sum( SpcAbs, axis=-1 ) # average the spectra (= SpcTot / NH_size) SpcAvg = np.mean(SpcAbs, axis=-1) SpcStd = np.std(SpcAbs, axis=-1) # build XSPEC readable table tab = np.zeros([ len(E_keV)-1, 3 ]) # left borders of energy bins tab[:,0] = E_keV[:-1] # right borders of energy bins tab[:,1] = E_keV[1:] # tabulate the spectrum in exp format (negative of the exponent) tab[:,2] = -np.log( SpcTot[:-1] ) tab[:,2][ tab[:,2] != tab[:,2] ] = 1e+30 tab[:,2][ tab[:,2] == float( 'inf' ) ] = 1e+29 if args.plot : fig, ( ax1, ax2 ) = plt.subplots( 2, 1 ) # plot the N spectra for n in range( NH_size ) : # mask out too small numbers (avoids numerical issues) mask_n = ( SpcAbs[:, n] > 1e-200 ) ax1.plot( E_keV[mask_n], np.log( SpcAbs[:,n][mask_n] ), 'b--', lw=0.3, alpha=0.3) # compute the ratio between the correct and the approximated spectra for the multiplicative model ratio = SpcTot / (SpcNHavg * NH_size) mask_avgNH = ( SpcNHavg > 1e-200 ) ax1.plot( E_keV[mask_avgNH], np.log( SpcNHavg[mask_avgNH] ) , 'k--', label=r'$\tilde{A}_\Theta = \exp(-\sigma \langle N_H \rangle)$' ) ax1.plot( E_keV, np.log( SpcTot ), 'b-.', lw=2, label=r'$A_\Theta = \sum_{i=1}^M \exp(-\sigma N_{H,i})$' ) ax2.plot( E_keV[mask_avgNH], np.log( ratio[mask_avgNH] ), 'k-' ) ax1.set_title( 'M={}, from '.format( NH_size ) + args.input_NH, fontsize=12 ) ax2.set_xlabel( 'E [keV]', fontsize=14 ) ax1.set_ylabel( 'ln R [ph/s]', fontsize=14 ) ax2.set_ylabel( 'ratio', fontsize=14 ) ax1.legend( fontsize=14 ) ax1.set_xscale('log') ax1.set_xlim([1e-1, 1e1]) ax2.set_xscale('log') ax2.set_xlim([1e-1, 1e1]) ax1.set_ylim([-4,0.1]) plt.tight_layout( True ) plt.show() if args.save or args.xspec : np.savetxt( args.output, tab ) print( '---> created file ', args.output ) if args.xspec : input_extension = args.input_NH.split( '.' )[-1] root_name = args.output.replace( '.' + input_extension, '' ) fits_name = root_name + '.fits' # to create the .fits XSPEC-readable table use : os.system( 'echo "ftflx2tab %s disnht %s additive=False" > convert_tab_to_etable.xcm' % ( args.output, fits_name ) ) os.system( 'xspec convert_tab_to_etable.xcm' ) print( '---> created file ' + fits_name ) os.system( 'rm convert_tab_to_etable.xcm' ) if args.verbose : print( 'cleaning temporary files' ) print( '' ) print( 'example : ' ) print( ' xspec> model pow * etable{%s}' % ( fits_name ) ) print( '' ) #-------------------------------------------------------------------------------------------- # RUN if __name__ == '__main__' : main()