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.

Chandra Orion Ultradeep project

Python : Make input coordinate table

Coordinates are relative to pointing direction in arcmin

# 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.

'''Make input coordinate table

Coordinates are relative to pointing direction in arcmin'''
import os
from astropy.table import Table
from astropy.io import fits

asolfile = self.get_data_file('asol')
asol = fits.getheader(asolfile, 1)
coup = Table.read(os.path.join(self.pkg_data, 'COUP.tsv'),
                  format='ascii.fast_tab')
tab = Table()
tab['RA'] = (coup['RAJ2000'] - asol['RA_NOM']) * 60
tab['DEC'] = (coup['DEJ2000'] - asol['DEC_NOM']) * 60
tab['weight'] = 10**(coup['Lt'] - 27)
tab['energy'] = coup['<E>']
tab.write('coup.marxin', format='ascii.no_header', overwrite=True)

C source code : C code for a grid of sources.

(user.h and jdmath.h are shipped with marx.)

#include <stdio.h>
#include <stdlib.h>
#include <jdmath.h>
#include "user.h"

/* This user source implements many point sources via a file that
 * specifies the source positions and energies.  The current implementation
 * assumes the format:
 *  RA  Dec weight energy
 * Here RA, Dec specifiy the source position, weight specifies the strength
 * of the source in relation to the others.
 */
typedef struct
{
   double cosx, cosy, cosz;
   double weight;
   double energy;
}
Point_Source_Type;

static unsigned int Num_Points;
static Point_Source_Type *Point_Sources;
static unsigned int Max_Num_Points;

static char *do_realloc (char *p, unsigned int len)
{
   if (p == NULL)
     p = malloc (len);
   else
     p = realloc (p, len);

   if (p == NULL)
     fprintf (stderr, "Not enough memory\n");

   return p;
}

static void free_sources (void)
{
   if (Point_Sources == NULL)
     return;

   free ((char *) Point_Sources);
   Point_Sources = NULL;
}

static int add_source (double ra, double dec, double weight, double energy)
{
   Point_Source_Type *p;
   double cosx, cosy, cosz;

   /* Convert to God's units from arc-min */
   ra = ra * (PI/(180.0 * 60.0));
   dec = dec * (PI/(180.0 * 60.0));

   if (Max_Num_Points == Num_Points)
     {
        Max_Num_Points += 32;
        p = (Point_Source_Type *)do_realloc ((char *)Point_Sources, Max_Num_Points * sizeof (Point_Source_Type));
        if (p == NULL)
          {
             free_sources ();
             return -1;
          }
        Point_Sources = p;
     }

   p = Point_Sources + Num_Points;
   /* Note the the minus sign is to generate a vector pointing from the
    * source to the origin
    */
   p->cosx = -cos (dec) * cos (ra);
   p->cosy = -cos (dec) * sin(ra);
   p->cosz = -sin (dec);

   p->weight = weight;
   p->energy = energy;
   Num_Points += 1;

   return 0;
}

static void normalize_sources (void)
{
   double total;
   unsigned int i;

   total = 0;
   for (i = 0; i < Num_Points; i++)
     {
     Point_Sources[i].weight += total;
     total = Point_Sources[i].weight;
     }

   for (i = 0; i < Num_Points; i++)
     Point_Sources[i].weight /= total;

   /* Make sure no round-off error affects the weight of the last point */
   Point_Sources[Num_Points - 1].weight = 1.0;
}

int user_open_source (char **argv, int argc, double area,
                   double cosx, double cosy, double cosz)
{
   FILE *fp;
   char line[1024];
   char *file;
   unsigned int linenum;

   file = argv[0];
   if (file == NULL)
     {
     fprintf (stderr, "UserSource Model requires FILE as argument\n");
     return -1;
     }

   fp = fopen (file, "r");
   if (fp == NULL)
     {
     fprintf (stderr, "Unable to open %s\n", file);
     return -1;
     }

   linenum = 0;
   while (NULL != fgets (line, sizeof (line), fp))
     {
     double ra, dec, weight, energy;

     linenum++;
     if (4 != sscanf (line, "%lf %lf %lf %lf", &ra, &dec, &weight, &energy))
       continue;

     if (weight <= 0.0)
       {
          fprintf (stderr, "weight on line %d of %s must be positive\n",
                   linenum, file);
          free_sources ();
          return -1;
       }

     if (-1 == add_source (ra, dec, weight, energy))
       {
          fclose (fp);
          return -1;
       }
     }

   fclose (fp);
   if (Num_Points == 0)
     {
     fprintf (stderr, "%s contains no sources\n", file);
     return -1;
     }

   normalize_sources ();
   return 0;
}

void user_close_source (void)
{
   free_sources ();
}


int user_create_ray (double *delta_t, double *energy,
                  double *cosx, double *cosy, double *cosz)
{
   double r;
   Point_Source_Type *p;

   p = Point_Sources;

   r = JDMrandom ();
   while (r > p->weight)
     p++;

   *delta_t = -1.0;
   *energy = p->energy;
   *cosx = p->cosx;
   *cosy = p->cosy;
   *cosz = p->cosz;

   return 0;
}

int main (int a, char **b)
{
   (void) a;
   (void) b;
   return 1;
}

Python : compile USER code

marx ships with a few examples of user sources. We pick one of them, copy them to the right directory and compile it with gcc.

# 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.

'''compile USER code

|marx| ships with a few examples of user sources. We pick one
of them, copy them to the right directory and compile it with gcc.
'''
marxpath = self.conf.get('marx', 'path')
src = os.path.join(marxpath,
                   'share', 'doc', 'marx', 'examples', 'user-source')
shutil.copy(os.path.join(src, 'user.h'),
            os.path.join(self.basepath, 'user.h'))

jdmath_h = os.path.join(marxpath, 'include')
jdmath_a = os.path.join(marxpath, 'lib', 'libjdmath.a')

subprocess.call(['gcc',
                 '-shared', 'pnts.c', '-o', 'pnts.so', '-fPIC',
                 '-I' + jdmath_h, jdmath_a])

shell : Unzip fits file.

MARX cannot read zipped fits files, so we need to unzip the .fits.gz asol files that we downloaded from the archive. On the other hand, CIAO tools work on both zipped or unzipped files, so there is no need to unzip all of them, just the files that MARX reads as input.

gunzip -f download/3744/primary/pcadf158603295N003_asol1.fits

marx : run marx USER source matching observation

marx RA_Nom=83.821008667032 Dec_Nom=-5.3944057451075 Roll_Nom=310.90882089387 GratingType=NONE ExposureTime=166703.44446998835 DitherModel=FILE DitherFile=download/3744/primary/pcadf158603295N003_asol1.fits TStart=158603295.25227 ACIS_Exposure_Time=3.1 SourceRA=83.819583 SourceDEC=-5.39 DetectorType=ACIS-I DetOffsetX=0.0014398546217030406 DetOffsetZ=0.0050286306017994775 OutputDir=COUP SourceType=USER UserSourceFile=pnts.so UserSourceArgs=coup.marxin

marx2fits : turn into fits file

marx2fits --pixadj=EDSER COUP COUP.fits

CIAO : ds9 image of the PSF

In the observation, the brightest sources are piled-up. We don’t bother simulating this here, so we just set the scaling limits to bring out the fainter details and ignore the bright peaks.

ds9 -log -cmap heat download/3744/primary/acisf03744N003_evt2.fits.gz COUP.fits -scale limits 0 2000 -frame 1 -regions command 'text 5:35:15 -5:22:09 # text=Observation font="helvetica 24"' -frame 2 -regions command 'text 5:35:15 -5:22:09 # text=MARX font="helvetica 24"' -region load src.fits -saveimage /melkor/d1/guenther/marx/doc/source/tests/figures/ONC_ds9.png -exit

CIAO : Source detection

dmcopy "COUP.fits[EVENTS][bin x=2500:5500:2,y=2500:5500:2]" im.fits  option=image clobber=yes
mkpsfmap im.fits psf.map 1.4 ecf=0.5
celldetect im.fits src.fits psffile=psf.map clobber=yes

Python : Check position of detected sources

# 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.

'''Check position of detected sources'''
import numpy as np
import matplotlib.pyplot as plt
from astropy.table import Table
from astropy.coordinates import SkyCoord
from astropy.io import fits
src = Table.read('src.fits')
srcin = Table.read(os.path.join(self.pkg_data, 'COUP.tsv'),
                   format='ascii.fast_tab')

src_co = SkyCoord(src['RA'], src['DEC'], unit='deg')
srcin_co = SkyCoord(srcin['RAJ2000'], srcin['DEJ2000'], unit='deg')
idx, d2d, d3d = src_co.match_to_catalog_sky(srcin_co)

asolfile = self.get_data_file('asol')
asol = fits.getheader(asolfile, 1)
cen = SkyCoord(asol['RA_NOM'], asol['DEC_NOM'], unit='deg')
d = cen.separation(src_co).arcsec

fig = plt.figure()
ax1 = plt.subplot(111)
scat1 = ax1.scatter(d, d2d.arcsec, c=np.log10(src['NET_COUNTS']), lw=1)
ax1.set_xlabel('distance from aimpoint [arcsec]')
ax1.set_ylabel('coordinate error [arcsec]')
ax1.set_xlim([0, 350])
ax1.set_ylim([0, 2])
cbar1 = fig.colorbar(scat1, ax=ax1)
cbar1.set_label('log(net counts per source)')

fig.savefig(self.figpath(list(self.figures.keys())[1]))