Note
The following shows the code used to run this marx test. You can inspect
it and adapt it to your needs, but you cannot copy and paste it directly
because it depends on local $PATH and other environment variables.
For example, we use a python function to manage the directory structure
for all the images generated by all the tests instead of giving the file
name directly to save images.
On-axis PSF at different energies¶
CIAO : Set up default marx.par file¶
cp /nfs/melkor/d1/guenther/marx/installed/dev//share/marx/pfiles/marx.par marx.par
pset marx.par GratingType=NONE
pset marx.par DetectorType=HRC-I
pset marx.par DitherModel=INTERNAL
pset marx.par ExposureTime=10000
pset marx.par SourceFlux=0.2
pset marx.par SourceRA=0.
pset marx.par SourceDEC=0.
pset marx.par RA_Nom=0.
pset marx.par Dec_Nom=0.
pset marx.par Roll_Nom=0.
marx : run marx for different energies¶
marx MinEnergy=0.25 MaxEnergy=0.25 OutputDir=marx0.25
marx MinEnergy=0.5 MaxEnergy=0.5 OutputDir=marx0.5
marx MinEnergy=1 MaxEnergy=1 OutputDir=marx1
marx MinEnergy=2 MaxEnergy=2 OutputDir=marx2
marx MinEnergy=3 MaxEnergy=3 OutputDir=marx3
marx MinEnergy=4 MaxEnergy=4 OutputDir=marx4
marx MinEnergy=6 MaxEnergy=6 OutputDir=marx6
marx MinEnergy=8 MaxEnergy=8 OutputDir=marx8
marxasp : Generate asol file¶
Since all simulations use the same pointing and exposure time, it is
enough to run marxasp once.
marxasp MarxDir=marx0.25
marx2fits : Fits files from marx runs¶
marx2fits --pixadj=NONE marx0.25 marx0.25.fits
marx2fits --pixadj=NONE marx0.5 marx0.5.fits
marx2fits --pixadj=NONE marx1 marx1.fits
marx2fits --pixadj=NONE marx2 marx2.fits
marx2fits --pixadj=NONE marx3 marx3.fits
marx2fits --pixadj=NONE marx4 marx4.fits
marx2fits --pixadj=NONE marx6 marx6.fits
marx2fits --pixadj=NONE marx8 marx8.fits
SAOTrace : Now use SAOTrace¶
trace-nest tag=sao0.25 srcpars="point{ position = { ra = 0., dec = 0., ra_aimpt=0., dec_aimpt=0. }, spectrum = { { 0.25, 0.2 } } } roll(0) dither_asol_marx{ file = 'marx_asol1.fits', ra = 0., dec = 0., roll = 0. }" tstart=362912400.01 limit=9999.98302428007 limit_type=sec
trace-nest tag=sao0.5 srcpars="point{ position = { ra = 0., dec = 0., ra_aimpt=0., dec_aimpt=0. }, spectrum = { { 0.5, 0.2 } } } roll(0) dither_asol_marx{ file = 'marx_asol1.fits', ra = 0., dec = 0., roll = 0. }" tstart=362912400.01 limit=9999.98302428007 limit_type=sec
trace-nest tag=sao1 srcpars="point{ position = { ra = 0., dec = 0., ra_aimpt=0., dec_aimpt=0. }, spectrum = { { 1, 0.2 } } } roll(0) dither_asol_marx{ file = 'marx_asol1.fits', ra = 0., dec = 0., roll = 0. }" tstart=362912400.01 limit=9999.98302428007 limit_type=sec
trace-nest tag=sao2 srcpars="point{ position = { ra = 0., dec = 0., ra_aimpt=0., dec_aimpt=0. }, spectrum = { { 2, 0.2 } } } roll(0) dither_asol_marx{ file = 'marx_asol1.fits', ra = 0., dec = 0., roll = 0. }" tstart=362912400.01 limit=9999.98302428007 limit_type=sec
trace-nest tag=sao3 srcpars="point{ position = { ra = 0., dec = 0., ra_aimpt=0., dec_aimpt=0. }, spectrum = { { 3, 0.2 } } } roll(0) dither_asol_marx{ file = 'marx_asol1.fits', ra = 0., dec = 0., roll = 0. }" tstart=362912400.01 limit=9999.98302428007 limit_type=sec
trace-nest tag=sao4 srcpars="point{ position = { ra = 0., dec = 0., ra_aimpt=0., dec_aimpt=0. }, spectrum = { { 4, 0.2 } } } roll(0) dither_asol_marx{ file = 'marx_asol1.fits', ra = 0., dec = 0., roll = 0. }" tstart=362912400.01 limit=9999.98302428007 limit_type=sec
trace-nest tag=sao6 srcpars="point{ position = { ra = 0., dec = 0., ra_aimpt=0., dec_aimpt=0. }, spectrum = { { 6, 0.2 } } } roll(0) dither_asol_marx{ file = 'marx_asol1.fits', ra = 0., dec = 0., roll = 0. }" tstart=362912400.01 limit=9999.98302428007 limit_type=sec
trace-nest tag=sao8 srcpars="point{ position = { ra = 0., dec = 0., ra_aimpt=0., dec_aimpt=0. }, spectrum = { { 8, 0.2 } } } roll(0) dither_asol_marx{ file = 'marx_asol1.fits', ra = 0., dec = 0., roll = 0. }" tstart=362912400.01 limit=9999.98302428007 limit_type=sec
marx : run marx for all SAOTrace runs¶
marx SourceType=SAOSAC SAOSACFile=sao0.25.fits OutputDir=saomarx0.25 DitherModel=FILE DitherFile=marx_asol1.fits
marx SourceType=SAOSAC SAOSACFile=sao0.5.fits OutputDir=saomarx0.5 DitherModel=FILE DitherFile=marx_asol1.fits
marx SourceType=SAOSAC SAOSACFile=sao1.fits OutputDir=saomarx1 DitherModel=FILE DitherFile=marx_asol1.fits
marx SourceType=SAOSAC SAOSACFile=sao2.fits OutputDir=saomarx2 DitherModel=FILE DitherFile=marx_asol1.fits
marx SourceType=SAOSAC SAOSACFile=sao3.fits OutputDir=saomarx3 DitherModel=FILE DitherFile=marx_asol1.fits
marx SourceType=SAOSAC SAOSACFile=sao4.fits OutputDir=saomarx4 DitherModel=FILE DitherFile=marx_asol1.fits
marx SourceType=SAOSAC SAOSACFile=sao6.fits OutputDir=saomarx6 DitherModel=FILE DitherFile=marx_asol1.fits
marx SourceType=SAOSAC SAOSACFile=sao8.fits OutputDir=saomarx8 DitherModel=FILE DitherFile=marx_asol1.fits
marx2fits : Fits files from marx + SAOTrace runs¶
marx2fits --pixadj=NONE saomarx0.25 saomarx0.25.fits
marx2fits --pixadj=NONE saomarx0.5 saomarx0.5.fits
marx2fits --pixadj=NONE saomarx1 saomarx1.fits
marx2fits --pixadj=NONE saomarx2 saomarx2.fits
marx2fits --pixadj=NONE saomarx3 saomarx3.fits
marx2fits --pixadj=NONE saomarx4 saomarx4.fits
marx2fits --pixadj=NONE saomarx6 saomarx6.fits
marx2fits --pixadj=NONE saomarx8 saomarx8.fits
Python : Image gallery of PSFs¶
# Note that this code might not run if you directly copy and paste it:
# - Not all import statements are shown here
# - `self` is a reference to a test instance, which allows access to
# parameters such as the directory where the test is run etc.
'''Image gallery of PSFs'''
import numpy as np
from matplotlib import pyplot as plt
from mpl_toolkits.axes_grid1 import AxesGrid
from matplotlib.ticker import StrMethodFormatter
from astropy.table import Table
fig = plt.figure(figsize=(10, 4))
grid = AxesGrid(fig, 111, # similar to subplot(142)
nrows_ncols=(2, len(self.parameter)),
axes_pad=0.0,
share_all=True,
label_mode="L",
cbar_location="top",
cbar_mode="single")
for i, e in enumerate(self.parameter):
for prog in ['marx', 'saomarx']:
f = os.path.join(self.basepath, '{0}{1}.fits'.format(prog, e))
with fits.open(f) as hdus:
header = hdus[1].header
tab = Table.read(f, hdu=1)
# The old question: Is an index the center of a pixel or the corner?
# Differs by 0.5...
d_ra = (tab['X'] - header['TCRPX9'] - 0.5) * header['TCDLT9'] * 3600.
# And do I count from the left or right (RA is reversed on the sky)?
d_dec = (tab['Y'] - header['TCRPX10'] + 0.5) * header['TCDLT10'] * 3600.
if prog == 'marx':
offset = 0
else:
offset = len(self.parameter)
im = grid[i + offset].hist2d(d_ra, d_dec,
range=[[-1, 1], [-1, 1]],
bins=[20, 20],
cmap=plt.get_cmap('OrRd'))
levels = np.max(im[0]) * np.array([0.3, 0.6, 0.9])
grid[i + offset].contour(0.5 * (im[1][1:] + im[1][:-1]),
0.5 * (im[2][1:] + im[2][:-1]),
# for some reason, the data is
# ordered for other orientation
im[0].T, levels=levels, colors='k') #,
#origin=im[3].origin)
grid[i + offset].grid()
if prog == 'marx':
grid[i].text(-.5, .5, '{0} keV'.format(e))
else:
grid[i].xaxis.set_major_formatter(StrMethodFormatter('{x:2.1g}'))
grid[0].yaxis.set_major_formatter(StrMethodFormatter('{x:2.1g}'))
grid[len(self.parameter)].yaxis.set_major_formatter(StrMethodFormatter('{x:2.1g}'))
grid.cbar_axes[0].colorbar(im[3])
for cax in grid.cbar_axes:
cax.toggle_label(False)
# This affects all axes as share_all = True.
grid.axes_llc.set_xticks([-.5, 0, .5])
grid.axes_llc.set_yticks([-.5, 0, .5])
grid[0].set_ylabel('Marx')
grid[offset].set_ylabel('Marx +\nSAOTrace')
grid[int(1.5 * offset)].set_xlabel('arcsec (measured from nominal source position)')
fig.subplots_adjust(top=1, bottom=0, left=0.1, right=0.99)
fig.savefig(self.figpath(list(self.figures.keys())[0]),
bbox_inches='tight')
Python : Plots of radial distribution¶
# Note that this code might not run if you directly copy and paste it:
# - Not all import statements are shown here
# - `self` is a reference to a test instance, which allows access to
# parameters such as the directory where the test is run etc.
'''Plots of radial distribution'''
import os
import numpy as np
from matplotlib import pyplot as plt
from astropy.table import Table
fig = plt.figure()
axecf = fig.add_subplot(111)
color = plt.cm.jet_r(np.linspace(0, 1, len(self.parameter)))
for i, e in enumerate(self.parameter):
for prog in ['marx', 'saomarx']:
f = os.path.join(self.basepath, '{0}{1}.fits'.format(prog, e))
with fits.open(f) as hdus:
header = hdus[1].header
tab = Table.read(f, hdu=1)
# The old question: Is an index the center of a pixel of the corner?
# Differs by 0.5...
d_ra = (tab['X'] - header['TCRPX9'] - 0.5) * header['TCDLT9'] * 3600.
# Count from the left or right (RA is reversed on the sky)?
d_dec = (tab['Y'] - header['TCRPX10'] + 0.5) * header['TCDLT10'] * 3600.
# The simulation is set up to make this simple: RA=DEC=0
# so cos(DEC) = 1 and we can approximate with Euklidian distance
r = np.linalg.norm(np.vstack([d_ra, d_dec]), axis=0)
val, edges = np.histogram(r, range=[0, 5], bins=50)
bin_mid_marx = 0.5 * (edges[:-1] + edges[1:])
ecf_marx = 1.0 * val.cumsum() / val.sum()
if prog == 'marx':
line, = axecf.plot(bin_mid_marx, ecf_marx, color=color[i],
label='{0} keV'.format(e))
else:
axecf.plot(bin_mid_marx, ecf_marx, color=line.get_color(), ls=':')
axecf.set_xscale('power', power=0.5)
axecf.legend(loc='lower right')
axecf.set_ylabel('encircled count fraction')
axecf.set_xlabel('radius [arcsec]')
axecf.set_xticks([0, .1, .2, .4, .6, .8, 1, 2, 3, 4, 5])
fig.savefig(self.figpath('ECF'), bbox_inches='tight')
