# # This is a script containing Sherpa Python commands for calculating # the K-corrections for a particular spectral model (in this case the thermal # plasma Mekal model) for a range of redshifts and parameter values (a range of # plasma temperatures). # # This script may be executed at the Sherpa prompt with the command # 'execfile("plot_kcorrs.py")', or at the CIAO prompt with # 'sherpa plot_kcorrs.py'. # # ================================================================== #Clear existing data and model definitions in the current Sherpa session: clean() #Define the range of redshifts in which you are interested, e.g. using the #Numpy 'arange' function: z = np.arange(0.0, 2.01, 0.01) #Define an energy grid which covers at least 0.5-2.0 keV for the #redshift range z = 0 to 2: elo = 0.5 ehi = 2.0 dataspace1d(0.1, 7.0, 0.01, id=1) #Define the model: set_source(xsmekal.plasma) plasma.abundanc = 0.3 #Calculate the &kcorr; for a 1 keV plasma model: plasma.kt = 1 k1 = calc_kcorr(z, elo, ehi, id=1) #Calculate the &kcorr; for a 2 keV plasma model: plasma.kt = 2 k2 = calc_kcorr(z, elo, ehi, id=1) #Calculate the &kcorr; for a 4 keV plasma model: plasma.kt = 4 k4 = calc_kcorr(z, elo, ehi, id=1) #Calculate the &kcorr; for a 10 keV plasma model: plasma.kt = 10 k10 = calc_kcorr(z, elo, ehi, id=1) #Calculate the analytic solution using a power law, #using the Numpy Python package (imported as 'np' or 'numpy' in Sherpa) #for mathematical operations on arrays: z = np.array(z) kpl = 1./np.sqrt(1+z) #Plot the model curves: clear() add_curve(z, k1) add_curve(z, k2) add_curve(z, k4) add_curve(z, k10) current_curve("all") set_curve(["symbol.style", "none"]) add_curve(z, kpl) set_curve(["symbol.style", "none"]) set_curve(["line.style", "shortdash"]) set_curve(["line.color","red"]) set_plot_xlabel("Redshift") set_plot_ylabel("K_{0.5-2.0 keV}") set_plot_title("xsmekal model, abundance = 0.3 solar") limits(X_AXIS, 0, 2) limits(Y_AXIS, 0, 3) add_label(1.2, 2.0, "kT = 1 keV") add_label(1.2, 1.3, "kT = 2 keV") add_label(1.5, 0.9, "kT = 4 keV") add_label(1.0, 0.5, "kT = 10 keV") print_window("plasmakT_kcorr.png") save("sherpa_plot_session.save") #finished quit()