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.

Simulated off-axis PSF

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.
pset marx.par MinEnergy=4
pset marx.par MaxEnergy=4

marx : run marx for different energies

marx SourceRA=0.0 OutputDir=marx0
marx SourceRA=0.00833333333333 OutputDir=marx0.5
marx SourceRA=0.0166666666667 OutputDir=marx1
marx SourceRA=0.0333333333333 OutputDir=marx2
marx SourceRA=0.05 OutputDir=marx3
marx SourceRA=0.0833333333333 OutputDir=marx5
marx SourceRA=0.116666666667 OutputDir=marx7
marx SourceRA=0.166666666667 OutputDir=marx10
marx SourceRA=0.25 OutputDir=marx15

marxasp : Generate asol file

Since all simulations use the same pointing and exposure time, it is enough to run marxasp once.

marxasp MarxDir=marx0

marx2fits : Fits files from marx runs

marx2fits --pixadj=NONE marx0 marx0.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 marx5 marx5.fits
marx2fits --pixadj=NONE marx7 marx7.fits
marx2fits --pixadj=NONE marx10 marx10.fits
marx2fits --pixadj=NONE marx15 marx15.fits

SAOTrace : Now use SAOTrace

trace-nest tag=sao0 srcpars="point{ position = { ra = 0.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.98373971 limit_type=sec
trace-nest tag=sao0.5 srcpars="point{ position = { ra = 0.00833333333333, 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.98373971 limit_type=sec
trace-nest tag=sao1 srcpars="point{ position = { ra = 0.0166666666667, 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.98373971 limit_type=sec
trace-nest tag=sao2 srcpars="point{ position = { ra = 0.0333333333333, 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.98373971 limit_type=sec
trace-nest tag=sao3 srcpars="point{ position = { ra = 0.05, 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.98373971 limit_type=sec
trace-nest tag=sao5 srcpars="point{ position = { ra = 0.0833333333333, 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.98373971 limit_type=sec
trace-nest tag=sao7 srcpars="point{ position = { ra = 0.116666666667, 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.98373971 limit_type=sec
trace-nest tag=sao10 srcpars="point{ position = { ra = 0.166666666667, 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.98373971 limit_type=sec
trace-nest tag=sao15 srcpars="point{ position = { ra = 0.25, 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.98373971 limit_type=sec

marx : run marx for all SAOTrace runs

marx DitherFile=marx_asol1.fits SAOSACFile=sao0.fits DitherModel=FILE SourceType=SAOSAC OutputDir=saomarx0
marx DitherFile=marx_asol1.fits SAOSACFile=sao0.5.fits DitherModel=FILE SourceType=SAOSAC OutputDir=saomarx0.5
marx DitherFile=marx_asol1.fits SAOSACFile=sao1.fits DitherModel=FILE SourceType=SAOSAC OutputDir=saomarx1
marx DitherFile=marx_asol1.fits SAOSACFile=sao2.fits DitherModel=FILE SourceType=SAOSAC OutputDir=saomarx2
marx DitherFile=marx_asol1.fits SAOSACFile=sao3.fits DitherModel=FILE SourceType=SAOSAC OutputDir=saomarx3
marx DitherFile=marx_asol1.fits SAOSACFile=sao5.fits DitherModel=FILE SourceType=SAOSAC OutputDir=saomarx5
marx DitherFile=marx_asol1.fits SAOSACFile=sao7.fits DitherModel=FILE SourceType=SAOSAC OutputDir=saomarx7
marx DitherFile=marx_asol1.fits SAOSACFile=sao10.fits DitherModel=FILE SourceType=SAOSAC OutputDir=saomarx10
marx DitherFile=marx_asol1.fits SAOSACFile=sao15.fits DitherModel=FILE SourceType=SAOSAC OutputDir=saomarx15

marx2fits : Fits files from marx + SAOTrace runs

marx2fits --pixadj=NONE saomarx0 saomarx0.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 saomarx5 saomarx5.fits
marx2fits --pixadj=NONE saomarx7 saomarx7.fits
marx2fits --pixadj=NONE saomarx10 saomarx10.fits
marx2fits --pixadj=NONE saomarx15 saomarx15.fits

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.

def radial_distribution(x, y):
    '''Calculate the radial distance from the mean position for events.

    Parameters
    ----------
    x, y : np.array
        x and y coordinates of events

    Returns
    -------
    r : np.array
        radial distance from mean position of the events.
    '''
    centx = x - astropy.stats.sigma_clipped_stats(x, sigma=1.5, iters=20)[0]
    centy = y - astropy.stats.sigma_clipped_stats(y, sigma=1.5, iters=20)[0]

    xy = np.vstack([centx, centy])
    return np.linalg.norm(xy, axis=0)

'''Plots of radial distribution'''
import os
import numpy as np
from matplotlib import pyplot as plt
from astropy.table import Table

ecf = [0.5, 0.7, 0.9]
fig = plt.figure()
axecf = fig.add_subplot(111)
out = np.zeros((len(self.parameter), 2, len(ecf)))

for i, e in enumerate(self.parameter):

    for prog in ['marx', 'saomarx']:
        tab = Table.read(os.path.join(self.basepath,
                                      '{0}{1}.fits'.format(prog, e)),
                         hdu=1)
        # The old question: Is an index the center of a pixel of the corner?
        # Differs by 0.5...
        d_ra = (tab['X'] - tab.meta['TCRPX9'] - 0.5) * tab.meta['TCDLT9'] * 3600.
        # Count from the left or right (RA is reversed on the sky)?
        d_dec = (tab['Y'] - tab.meta['TCRPX10'] + 0.5) * tab.meta['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 = radial_distribution(d_ra, d_dec)
        if prog == 'marx':
            j = 0
        else:
            j = 1
        out[i, j, :] = np.percentile(r, np.array(ecf) * 100.)

for i, e in enumerate(ecf):
    line, = axecf.plot(self.parameter, out[:, 0, i],
                       label='{0} %'.format(e*100))
    axecf.plot(self.parameter, out[:, 1, i], color=line.get_color(),
               ls=':')
axecf.legend(loc='lower right')
axecf.set_ylabel('radius [arcsec]')
axecf.set_xlabel('off-axis angle [arcmin]')
axecf.grid()
fig.savefig(self.figpath('ECF'))