Skip to content

Instantly share code, notes, and snippets.

@ZGainsforth
Created July 23, 2021 19:05
Show Gist options
  • Save ZGainsforth/c012bed5ee9ff97aa282e85863985958 to your computer and use it in GitHub Desktop.
Save ZGainsforth/c012bed5ee9ff97aa282e85863985958 to your computer and use it in GitHub Desktop.
Convert Newton/XMM RGS Fluxed spectrum into a text file which can be plotted easily.
from astropy.constants import c, h, e
from astropy.io import fits
import numpy as np
from sys import argv
def LoadSpectrum(SpectrumName='0090340201.fluxed'):
spectrumfits = fits.open(SpectrumName)
spec = spectrumfits[1]
angstrom = np.nan_to_num(spec.data.field('CHANNEL').copy())
eV = np.nan_to_num((h*c/(e.si*angstrom*1e-10)).value)
flux = np.nan_to_num(spec.data.field('FLUX').copy())
error = np.nan_to_num(spec.data.field('ERROR').copy())
spectrumfits.close()
# x-axis(lambda), x-axis(eV), y-axis, x-label, x-label, y-label
return angstrom.astype('float'), eV.astype('float'), flux.astype('float'), error.astype('float'), spec.header['TUNIT1'], 'eV', spec.header['TUNIT2'], spec.header['TUNIT3']
FileName = f'{argv[1]}.txt'
angstromsum, eVsum, fluxsum, errorsum, angstrom_label, eV_label, flux_label, error_label = LoadSpectrum(f'{argv[1]}.fluxed')
assert flux_label == '1/(s cm^2 A)', 'Flux units are not what we expect before converting to Chandra units of 1/(s cm^2 keV)'
assert error_label == '1/(s cm^2 A)', 'Flux error units are not what we expect before converting to Chandra units of 1/(s cm^2 keV)'
fluxchandra = fluxsum * 8.07e-2 * angstromsum**2
flux_label = '1/(s cm^2 keV)'
fluxerrorchandra = errorsum * 8.07e-2 * angstromsum**2
error_label = '1/(s cm^2 keV)'
# Now we have everything we need to write the file.
with open(FileName, 'w') as f:
f.write(f'# keV, flux in {flux_label:>12s}, flux error in {error_label:>12s}, Counts, {angstrom_label:>20s}, {eV_label:>19s}\n')
for i in reversed(range(len(eVsum))):
f.write(f'{eVsum[i]/1000:21.6f} {fluxchandra[i]:20.6f} {fluxerrorchandra[i]:21.6f} {1.000000:20.6f} {angstromsum[i]:20.6f} {eVsum[i]:21.6f}\n')
print(f'Wrote {FileName}')
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment