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.

Regular grid (HRC)

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

(user.h is shipped with marx.)

#include <stdio.h>
#include <math.h>
#include "user.h"

static double Source_CosX;
static double Source_CosY;
static double Source_CosZ;

int user_open_source (char **argv, int argc, double area,
                   double cosx, double cosy, double cosz)
{
   return 0;
}

void user_close_source (void)
{
}

static double To_Radians = (M_PI / 180.0 / 3600.0);
#define ARC_SECONDS_PER_CELL 50
#define ANGULAR_STEPS 16

int user_create_ray (double *delta_t, double *energy,
                  double *cosx, double *cosy, double *cosz)
{
   static int last_i = 0;
   static int last_j = 0;
   double theta, phi;
   double cos_theta, sin_theta;

   if (last_j == ANGULAR_STEPS){
        last_j = 0;
        last_i++;
   }
   if (last_i == 20) last_i = 0;

   theta = To_Radians * last_i * ARC_SECONDS_PER_CELL;
   phi = (10. /180 * M_PI) + last_j * 2 * M_PI / ANGULAR_STEPS;

   sin_theta = sin(theta);

   *cosx = -cos (theta);
   *cosy = sin_theta * cos (phi);
   *cosz = sin_theta * sin (phi);

   *delta_t = -1.0;
   *energy = -1.0;

   if (last_i ==0){
     last_i++;
        }
   else {
     last_j++;
        }

   return 0;
}

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

Python : compile USER code

# 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'''
marxpath = self.conf.get('marx', 'path')
src = os.path.join(marxpath, 'share', 'doc', 'marx', 'examples',
                   'user-source', 'user.h')
shutil.copy(os.path.join(src),
            os.path.join(self.basepath, 'user.h'))

subprocess.call(['gcc', '-lm', '-fPIC',
                 '-shared', 'radialgrid.c', '-o', 'radialgrid.so'])

marx : run USER source

marx SourceType=USER OutputDir=points GratingType=NONE SourceRA=90.0 SourceDEC=0.0 RA_Nom=90.0 Dec_Nom=0 Roll_Nom=0 DetectorType=HRC-I UserSourceFile=radialgrid.so NumRays=-100000 ExposureTime=0

marx2fits : turn into fits file

marx2fits --pixadj=EDSER points points.fits

CIAO : ds9 image of the PSF

ds9 -width 500 -height 500 -log -cmap heat points.fits -pan to 16392 16392 physical -bin factor 16 -grid -saveimage /melkor/d1/guenther/marx/doc/source/tests/figures/RegularGridHRCI_ds9.png -exit

CIAO : Source detection

dmcopy "points.fits[EVENTS][bin x=8500:24500:8,y=8500:24500:8]" im.fits  option=image clobber=yes
mkpsfmap im.fits psf.map 1.4 ecf=0.5 clobber=yes
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
src = Table.read('src.fits')
# Find distance from input position.
src['RA_INPUT'] = src['RA'] - (src['RA'] // (360./16.)) * (360./16.) - 10.
# Problem: Might expect source at 1.0,
# but measure at 0.99. In this case distance to next lower source
# is 0.99. Thus shift input by 0.005 (about 50 arcsec / 2)
# before integer devision
src['DEC_INPUT'] = src['DEC'] - ((0.005 + src['DEC']) // (50./3600.)) * (50./3600.)

cen = SkyCoord(90., 0, unit='deg')
det = SkyCoord(src['RA'], src['DEC'], unit='deg')
d = cen.separation(det).arcsec
d_err = np.mod(d + 10, 50.) - 10

ang = cen.position_angle(det).degree
# Subtract offset that we placed in the C code to avoid 0./360. ambiguity
# Step width is 360./16 = 22.5 deg
# Offset is 10 deg. Complement we find here is 12.5 deg.
ang = ang - 12.5

ang_err = np.mod(ang + 2, 360. / 16.) - 2
ind = d > 10

fig = plt.figure(figsize=(8, 4))
ax1 = plt.subplot(121)
scat1 = ax1.scatter(d, d_err, c=src['NET_COUNTS'], lw=1)
ax1.set_xlabel('distance [arcsec]')
ax1.set_ylabel('distance error [arcsec]')
ax1.set_xlim([-10, 620])
ax1.set_ylim([-1, 0.5])
ax2 = plt.subplot(122)
scat2 = ax2.scatter(ang, ang_err, c=src['NET_COUNTS'], lw=1)
ax2.set_xlabel('pos ang [deg]')
ax2.set_ylabel('pos ang error [deg]')
ax2.set_xlim([-5, 350])
ax2.set_ylim([-0.3, 0.3])
cbar2 = fig.colorbar(scat2, ax=ax2)
cbar2.set_label('net counts per source')

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