Comment calculer le rayonnement X (CXB) de fond d'espace avec Python en spécifiant la plage de flux

Contexte

Lors de l'estimation du rayonnement X (CXB) de fond d'espace, non seulement nous voulons connaître une partie de la luminosité totale, mais également combien lors de l'imagerie et de l'observation, à l'exclusion des sources de points lumineux détectés par l'imagerie. Reste sans décomposition spatiale? Besoin d'être estimé.

Ici, je vais vous montrer comment calculer soft (1-2 keV) et hard (2-10 keV) avec python selon https://iopscience.iop.org/article/10.1086/374335.

Code source

Pour ceux qui peuvent lire le code, voir Colab of calc_cxb.

python

#!/bin/evn python 

import matplotlib.pyplot as plt
import numpy as np

# THE RESOLVED FRACTION OF THE COSMIC X-RAY BACKGROUND A. Moretti, 
# The Astrophysical Journal, Volume 588, Number 2 (2003)
# https://iopscience.iop.org/article/10.1086/374335

smax = 1e-1 # maxinum of flux range (erg cm^-2 c^-1)
smin = 1e-17 # minimum of flux range (erg cm^-2 c^-1)
n = int(1e5) # number of grid in flux 
fcut = 5e-14 # user-specified maxinum flux (erg cm^-2 c^-1)

def calc_cxb(s,fcut=5e-14, ftotal_walker=2.18e-11, plot=True, soft=True):

	if soft: # 1-2 keV
		print("..... soft X-ray 1-2 keV is chosen")
		a1=1.82; a2=0.6; s0=1.48e-14; ns=6150
		foutname="soft"
	else: # hard 2-10 keV 
		print("..... hard X-ray 2-10 keV is chosen")
		a1=1.57; a2=0.44; s0=4.5e-15; ns=5300
		foutname="hard"


	n_larger_s = ns * 2.0e-15 ** a1 / (s**a1 + s0**(a1-a2)*s**a2) # Moretti et al. 2013, eq(2)

	if plot:
		F = plt.figure(figsize=(6,6))
		ax = plt.subplot(1,1,1)
		plt.plot(s, n_larger_s, ".", label="test")
		plt.xscale('log')
		plt.ylabel(r"$N(>S)/deg^2$")
		plt.xlabel(r"$S (erg cm^{-2} s^{-1})$")
		plt.yscale('log')
		plt.savefig(foutname + "_ns.png ")
		plt.show()

	diff_n_larger_s = np.abs(np.diff(n_larger_s))
	diff_s = np.diff(s)
	div_ns = np.abs(diff_n_larger_s/diff_s)

	s_lastcut = s[:-1]
	totalflux = np.sum(div_ns*s_lastcut*diff_s) # Moretti et al. 2013, eq(4)

	fluxcut = np.where(s > fcut)[0][:-1]
	particalflux = np.sum(div_ns[fluxcut]*s_lastcut[fluxcut]*diff_s[fluxcut])
	fdiff = totalflux - particalflux 
	fcxb = ftotal_walker - particalflux # Walker et al. 2016, eq(1)

	print("totalflux    = ", "%.4e"%totalflux, " [erg cm^-2 c^-1]", " from ", "%.4e"%smin, " to ",  "%.4e"%smax)
	print("particalflux = ", "%.4e"%particalflux, " [erg cm^-2 c^-1]", " from ", "%.4e"%fcut, " to ",  "%.4e"%smax)
	print("fdiff        = totalflux - particalflux     = ", "%.4e"%fdiff, " [erg cm^-2 c^-1]", " from ", "%.4e"%smin, " to ",  "%.4e"%fcut)	
	print("fcxb(2-10keV)= ftotal_walker - particalflux = ", "%.4e"%fcxb, " [erg cm^-2 c^-1]")

	if plot:
		F = plt.figure(figsize=(6,6))
		ax = plt.subplot(1,1,1)
		plt.plot(s[:-1], div_ns, ".", label="test")
		plt.xscale('log')
		plt.ylabel(r"$dN/dS/deg^2$")
		plt.xlabel(r"$S (erg cm^{-2} s^{-1})$")
		plt.yscale('log')
		plt.savefig(foutname + "_nsdiff.png ")		
		plt.show()

	return s, n_larger_s, diff_n_larger_s, div_ns, diff_s

#s = np.linspace(sexcl,smax,n)
s = np.logspace(np.log10(smin),np.log10(smax),n)

# plot soft band
s, n_larger_s, diff_n_larger_s, div_ns, diff_s = calc_cxb(s,fcut=fcut)
# plot hard band
s, n_larger_s, diff_n_larger_s, div_ns, diff_s = calc_cxb(s,fcut=fcut,soft=False)

Résultat d'exécution

1-2 keV

2-10 keV

Méthode de calcul

Comme prévu, si elle est linéaire, la mémoire est insuffisante et elle ne peut pas être calculée, donc dans le log de l'espace s (erg cm ^ -2 s ^ -1), l'intégration a utilisé une simple approximation rectangulaire.

Cette fonction est calculée par n_larger_s = ns * 2.0e-15 ** a1 / (s ** a1 + s0 ** (a1-a2) * s ** a2).

Cette intégration est calculée par totalflux = np.sum (div_ns * s_lastcut * diff_s).

fcxb = ftotal_walker --particalflux # Walker et al.2016, eq (1) est un bonus, et le CXB de 2-10 keV dans cet article est estimé, et le calcul est basé sur celui-ci.