diff options
Diffstat (limited to 'astro/marx')
| -rwxr-xr-x | astro/marx/marx_pntsrc.py | 121 | ||||
| -rwxr-xr-x | astro/marx/randpoints.py | 327 |
2 files changed, 448 insertions, 0 deletions
diff --git a/astro/marx/marx_pntsrc.py b/astro/marx/marx_pntsrc.py new file mode 100755 index 0000000..97b1575 --- /dev/null +++ b/astro/marx/marx_pntsrc.py @@ -0,0 +1,121 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +# +# Aaron LI +# 2015/06/16 + +""" +Run MARX simulation on a given list of point sources, merge the +output simulation results, and finally convert into FITS image. +""" + +__version__ = "0.1.0" +__date__ = "2015/06/16" + +import sys +import argparse +import subprocess +import re +import os + + +def marx_pntsrc(pfile, ra, dec, flux, outdir): + """ + Run MARX simulation for the provided point source. + """ + cmd = "marx @@%(pfile)s SourceRA=%(ra)s " % {"pfile": pfile, "ra": ra} + \ + "SourceDEC=%(dec)s SourceFlux=%(flux)s OutputDir=%(outdir)s" % \ + {"dec": dec, "flux": flux, "outdir": outdir} + print("CMD: %s" % cmd, file=sys.stderr) + subprocess.call(cmd, shell=True) + + +def marxcat(indirs, outdir): + """ + Concatenate a list of MARX simulation results. + + Note: the number of MARX results to be concatenated at *one* time + can not be to many, otherwise the 'marxcat' tool will failed. + """ + if isinstance(indirs, list): + pass + elif isinstance(indirs, str): + indirs = indirs.split() + else: + raise ValueError("invalid indirs type: %s" % indirs) + pid = os.getpid() + tempdir = "_marx_tmp%d" % pid + cmd = "cp -a %(marxdir)s %(tempdir)s" % \ + {"marxdir": indirs[0], "tempdir": tempdir} + print("CMD: %s" % cmd, file=sys.stderr) + subprocess.call(cmd, shell=True) + del indirs[0] + while len(indirs) > 0: + # concatenated 10 directories each time + catdirs = indirs[:9] + del indirs[:9] + catdirs = tempdir + " " + " ".join(catdirs) + # concatenate MARX results + cmd = "marxcat %(catdirs)s %(outdir)s" % \ + {"catdirs": catdirs, "outdir": outdir} + print("CMD: %s" % cmd, file=sys.stderr) + subprocess.call(cmd, shell=True) + # move output results to temporary directory + cmd = "rm -rf %(tempdir)s && mv %(outdir)s %(tempdir)s" % \ + {"tempdir": tempdir, "outdir": outdir} + print("CMD: %s" % cmd, file=sys.stderr) + subprocess.call(cmd, shell=True) + cmd = "mv %(tempdir)s %(outdir)s" % \ + {"tempdir": tempdir, "outdir": outdir} + print("CMD: %s" % cmd, file=sys.stderr) + subprocess.call(cmd, shell=True) + + +def marx2fits(indir, outfile, params=""): + """ + Convert the results of MARX simulation into FITS image. + """ + cmd = "marx2fits %(params)s %(indir)s %(outfile)s" % \ + {"params": params, "indir": indir, "outfile": outfile} + print("CMD: %s" % cmd, file=sys.stderr) + subprocess.call(cmd, shell=True) + + +def main(): + parser = argparse.ArgumentParser( + description="Run MARX on a given list of point sources") + parser.add_argument("-V", "--version", action="version", + version="%(prog)s " + "%s (%s)" % (__version__, __date__)) + parser.add_argument("pfile", help="marx paramter file") + parser.add_argument("srclist", help="point source list file") + parser.add_argument("outprefix", help="prefix of output directories") + args = parser.parse_args() + + outdirs = [] + i = 0 + for line in open(args.srclist, "r"): + if re.match(r"^\s*$", line): + # skip blank line + continue + elif re.match(r"^\s*#", line): + # skip comment line + continue + i += 1 + ra, dec, flux = map(float, line.split()) + print("INFO: ra = %g, dec = %g, flux = %g" % (ra, dec, flux), + file=sys.stderr) + outdir = "%sp%03d" % (args.outprefix, i) + print("INFO: outdir = %s" % outdir, file=sys.stderr) + outdirs.append(outdir) + marx_pntsrc(args.pfile, ra, dec, flux, outdir) + # merge results + merged = args.outprefix + "merged" + marxcat(outdirs, merged) + # convert to FITS image + merged_fits = merged + ".fits" + marx2fits(merged, merged_fits) + + +if __name__ == "__main__": + main() + diff --git a/astro/marx/randpoints.py b/astro/marx/randpoints.py new file mode 100755 index 0000000..da0bedc --- /dev/null +++ b/astro/marx/randpoints.py @@ -0,0 +1,327 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +# +# Aaron LI +# Created: 2015-12-03 +# Updated: 2015-12-05 +# +# ChangeLog: +# 2015-12-05: +# * Add class "RegionDS9" to parse region +# 2015-12-10: +# * Support both "erg/cm^2/s" and "photon/cm^2/s" flux units +# * Add argument "--outregion" to also save a region file +# + +""" +Generate random point source information for MARX simulation. +And make a point source data list for "marx_pntsrc.py" usage. +""" + +__version__ = "0.3.1" +__date__ = "2015-12-10" + + +import sys +import argparse +import re +import numpy as np + + +class RegionDS9: + """ + Process DS9 regions. + """ + def __init__(self, shape=None, xc=None, yc=None, + width=None, height=None, rotation=None): + self.shape = shape + self.xc = xc + self.yc = yc + self.width = width + self.height = height + self.rotation = rotation + + def parse(self, region): + """ + Parse a DS9 region string and update the instance. + + Region syntax: + box(xc,yc,width,height,rotation) + + Note: + "width", "height" may have a '"' suffix which means "arcsec" instead + of "degree". + """ + re_box = re.compile(r'^\s*(?P<shape>box)\(\s*(?P<xc>[\d.-]+)\s*,\s*(?P<yc>[\d.-]+)\s*,\s*(?P<width>[\d.-]+"?)\s*,\s*(?P<height>[\d.-]+"?)\s*,\s*(?P<rotation>[\d.-]+)\s*\).*$', re.I) + m_box = re_box.match(region) + if m_box is not None: + self.shape = "box" + self.xc = float(m_box.group("xc")) + self.yc = float(m_box.group("yc")) + self.width = self.parse_dms(m_box.group("width")) + self.height = self.parse_dms(m_box.group("height")) + self.rotation = float(m_box.group("rotation")) + else: + raise NotImplementedError("Only 'box' region supported") + + @staticmethod + def parse_dms(ms): + """ + Parse a value in format ?'?" into degree. + """ + re_arcmin = re.compile(r'^\s*(?P<arcmin>[\d.]+)\'.*') + re_arcsec = re.compile(r'^([^\']*\'|)\s*(?P<arcsec>[\d.]+)".*') + m_arcmin = re_arcmin.match(ms) + m_arcsec = re_arcsec.match(ms) + degree = 0.0 + if m_arcmin is not None: + degree += float(m_arcmin.group("arcmin")) / 60.0 + if m_arcsec is not None: + degree += float(m_arcsec.group("arcsec")) / 3600.0 + return degree + + +class RandCoord: + """ + Randomly generate the coordinates of point sources within a given box + region for MARX simulation. + + Arguments: + xc - central X position of the box (degree) + yc - central Y position of the box (degree) + width - width of the box (degree) + height - height of the box (degree) + mindist - minimum distance between each generated coordinate (degree) + """ + + def __init__(self, xc, yc, width, height, mindist=0): + self.xc = xc + self.yc = yc + self.width = width + self.height = height + self.mindist = mindist + # Record the generated coordinates: [(x1,y1), (x2,y2), ...] + self.xy = [] + + def clear(self): + """ + Clear previously generated coordinates. + """ + self.xy = [] + + def generate(self, n=1): + """ + Generate random coordinates. + """ + coord = [] + xmin = self.xc - 0.5 * self.width + xmax = self.xc + 0.5 * self.width + ymin = self.yc - 0.5 * self.height + ymax = self.yc + 0.5 * self.height + i = 0 + while i < n: + x = np.random.uniform(low=xmin, high=xmax) + y = np.random.uniform(low=ymin, high=ymax) + if self.checkDistance((x, y)): + i += 1 + coord.append((x, y)) + self.xy.append((x, y)) + return coord + + def checkDistance(self, coord): + """ + Check whether the given coordinate has a distance larger than + the specified "mindist" + """ + if len(self.xy) == 0: + return True + else: + xy = np.array(self.xy) # each row represents one coordinate + dist2 = (xy[:, 0] - coord[0])**2 + (xy[:, 1] - coord[1])**2 + if all(dist2 >= self.mindist**2): + return True + else: + return False + + +class RandFlux: + """ + Randomly generate the flux of point sources for MARX simulation. + + Arguments: + fmin - minimum flux + fmax - maximum flux + """ + + def __init__(self, fmin, fmax): + self.fmin = fmin + self.fmax = fmax + + @staticmethod + def fluxDensity(S): + """ + The *differential* number count - flux function: dN(>S)/dS + i.e., density function + + Broken power law: + dN/dS = (1) K (S/S_ref)^(-gamma_1); (S < S_b) + (2) K (S_b/S_ref)^(gamma_2-gamma_1) (S/S_ref)^(-gamma_2); (S >= S_b) + K: normalization constant + S_ref: normalization flux; [10^-15 erg/cm^2/s] + gamma_1: faint power-law index + gamma_2: bright power-law index + S_b: break flux; [10^-15 erg/cm^2/s] + + Reference: + [1] Kim et al. 2007, ApJ, 659, 29 + http://adsabs.harvard.edu/abs/2007ApJ...659...29K + http://hea-www.cfa.harvard.edu/CHAMP/PUBLICATIONS/ChaMP_ncounts.pdf + Table 4: ChaMP: 9.6 [deg^2]: 0.5-8 [keV]: 1.4 (photon index) + Differential number count; broken power law + K (normalization constant): 1557 (+28 / -50) + S_ref (normalization flux): 1.0 [10^-15 erg/cm^2/s] + gamma_1 (faint power-law index): 1.64 (+/- 0.01) + gamma_2 (bright power-law index): 2.48 (+/- 0.05) + S_b (break flux): 22.9 (+/- 1.6) [10^-15 erg/cm^2/s] + f_min (faint flux limit): 0.69 [10^-15 erg/cm^2/s] + f_max (bright flux limit): 6767.74 [10^-15 erg/cm^2/s] + """ + K = 1557 # normalization constant: 1557 (+28 / -50) + S_ref = 1.0 # normalization flux: 1.0 [10^-15 erg/cm^2/s] + gamma_1 = 1.64 # faint power-law index: 1.64 (+/- 0.01) + gamma_2 = 2.48 # bright power-law index: 2.48 (+/- 0.05) + S_b = 22.9 # break flux: 22.9 (+/- 1.6) [10^-15 erg/cm^2/s] + # Adjust unit/magnitude + S = S / 1e-15 # => unit: 10^-15 erg/cm^2/s + if isinstance(S, np.ndarray): + Np = np.zeros(S.shape) + Np[S<=0] = 0.0 + Np[S<=S_b] = K * (S[S<=S_b] / S_ref)**(-gamma_1) + Np[S>S_b] = K * (S_b/S_ref)**(gamma_2-gamma_1) * (S[S>S_b] / S_ref)**(-gamma_2) + else: + # "S" is a single number + if S <= 0.0: + Np = 0.0 + elif S <= S_b: + Np = K * (S/S_ref)**(-gamma_1) + else: + Np = K * (S_b/S_ref)**(gamma_2-gamma_1) * (S/S_ref)**(-gamma_2) + # + return Np + + def generate(self, n=1): + """ + Generate a sample of luminosity values within [min, max] from + the above luminosity distribution. + """ + results = [] + # Get the maximum value of the flux number density function, + # which is a monotonically decreasing. + M = self.fluxDensity(self.fmin) + for i in range(n): + while True: + u = np.random.uniform() * M + y = 10 ** np.random.uniform(low=np.log10(self.fmin), + high=np.log10(self.fmax)) + if u <= self.fluxDensity(y): + results.append(y) + break + return results + + +def main(): + parser = argparse.ArgumentParser( + description="Randomly generate point sources information for MARX") + parser.add_argument("-V", "--version", action="version", + version="%(prog)s " + "v%s (%s)" % (__version__, __date__)) + parser.add_argument("-n", "--number", dest="number", type=int, default=1, + help="number of point sources (default: 1)") + parser.add_argument("-m", "--fmin", dest="fmin", + type=float, default=1e-15, + help="minimum flux (default: 1e-15 erg/cm^2/s)") + parser.add_argument("-M", "--fmax", dest="fmax", + type=float, default=6000e-15, + help="maximum flux (default: 6000e-15 erg/cm^2/s)") + parser.add_argument("-r", "--region", dest="region", required=True, + help="region within which to generate coordinates ('box' only)") + parser.add_argument("-d", "--distance", dest="distance", default="0", + help="minimum distance between coordinates (default: 0) [unit: deg/arcmin]") + parser.add_argument("-u", "--unit", dest="unit", default="erg", + help="unit for input and output flux; 'erg' (default) / 'photon'") + parser.add_argument("-f", "--factor", dest="factor", type=float, + help="conversion factor from 'photon/cm^s/s' to 'erg/cm^2/s' (required if unit='photon')") + parser.add_argument("-o", "--outfile", dest="outfile", + help="output file to save the generate information list") + parser.add_argument("-O", "--outregion", dest="outregion", + help="write the generate information list as a DS9 region file") + + args = parser.parse_args() + + # Check flux unit + if args.unit == "erg": + unit = "erg/cm^2/s" + fmin = args.fmin + fmax = args.fmax + factor = 1.0 + elif args.unit == "photon": + unit = "photon/cm^2/s" + factor = args.factor + try: + fmin = args.fmin / factor + fmax = args.fmax / factor + except NameError: + raise ValueError("argument '--factor' required") + else: + raise ValueError("unsupported flux unit") + + region = RegionDS9() + region.parse(args.region) + # Check the box rotation + if not (abs(region.rotation) <= 1.0 or abs(region.rotation-360) <= 1.0): + raise NotImplementedError("rotated 'box' region not supported") + + # Minimum distance between generated coordinates + try: + mindist = float(args.distance) + except ValueError: + mindist = region.parse_dms(args.distance) + + randcoord = RandCoord(region.xc, region.yc, region.width, region.height, + mindist=mindist) + randflux = RandFlux(fmin, fmax) + coord = randcoord.generate(n=args.number) + flux = randflux.generate(n=args.number) + + if args.outfile: + outfile = open(args.outfile, "w") + else: + outfile = sys.stdout + + print("# region: %s" % args.region, file=outfile) + print("# mindist: %.9f [deg]" % mindist, file=outfile) + print("# f_min: %.9g; f_max: %.9g [%s]" % (fmin, fmax, unit), file=outfile) + print("# factor: %g [photon/cm^2/s] / [erg/cm^2/s]" % factor, file=outfile) + print("# R.A.[deg] Dec.[deg] Flux[%s]" % unit, file=outfile) + for ((ra, dec), f) in zip(coord, flux): + print("%.9f %.9f %.9g" % (ra, dec, f*factor), file=outfile) + + if args.outfile: + outfile.close() + + # Save the generated information as a DS9 region file if specified + if args.outregion: + reg_r = '3"' + reg_header = ["# Region file format: DS9 version 4.1", "fk5"] + regions = [ + "circle(%.9f,%.9f,%s) # text={%.9g}" % \ + (ra, dec, reg_r, f*factor) \ + for ((ra, dec), f) in zip(coord, flux) + ] + regfile = open(args.outregion, "w") + regfile.write("\n".join(reg_header + regions)) + regfile.close() + + +if __name__ == "__main__": + main() + |