Source specification now uses astropy Quantities

CHANGES.rst

@@ -34,6 +34,11 @@ New Features
 
 API Changes
 ^^^^^^^^^^^
+
+- Sources now require input for flux, energy, and polarization as astropy
+  quantities. Any table-like input must not be a astropy `QTable` and flux and
+  polarization values are *densities* not per-bin. [#218]
+
 - The geometry of a class is now defined as a separate object from the new
   `marxs.math.geometry` module. Before, the geometry was backed into the 
   optical components itself which led to significant code duplication. The new

docs/examples.rst

@@ -53,6 +53,7 @@ We now turn the combination of X-ray source and first mirror around the y-axis (
 
 .. doctest-requires:: mayavi
 
+  >>> import astropy.units as u
   >>> # Make an object that keeps the photon position after every simulation step
   >>> # for later plotting of photon paths.
   >>> pos = KeepCol('pos')
@@ -66,7 +67,7 @@ We now turn the combination of X-ray source and first mirror around the y-axis (
   ...   light.position = np.dot(rotmat, light_pos)
   ...   light.dir = np.dot(rotmat, light_dir)
   ...   m1.geometry.pos4d = np.dot(rotmat, m1pos4d)
-  ...   rays = light(100)
+  ...   rays = light(100 * u.s)
   ...   pos.data = []
   ...   pos(rays)
   ...   rays = experiment(rays)

docs/newopticalelements.rst

@@ -10,8 +10,9 @@ First, let us generate some photons coming from some astrophysical source::
 
   >>> from astropy.coordinates import SkyCoord
   >>> from marxs import source, optics
+  >>> import astropy.units as u
   >>> mysource = source.PointSource(coords=SkyCoord(181., -23., unit='deg'))
-  >>> photons = mysource.generate_photons(10)
+  >>> photons = mysource.generate_photons(10 * u.s)
   >>> mypointing = source.FixedPointing(coords=SkyCoord('12h01m23s -22d59m12s'))
   >>> photons = mypointing(photons)
 

docs/pyplots/runexample.py

@@ -1,28 +1,31 @@
 import numpy as np
 from matplotlib import pyplot as plt
 from astropy.coordinates import SkyCoord
+import astropy.units as u
+from astropy.table import QTable
 
 from marxs.source import poisson_process
 from marxs.missions import chandra
 from marxs.source import PointSource
 
 ngc1313_X1 = SkyCoord("3h18m19.99s -66d29m10.97s")
 
-energies = np.arange(.3, 8., .01)
-spectrum1 = 6.3e-4 * energies**(-1.9)
-flux1 = (spectrum1[1:] * np.diff(energies)).sum()
-flux1 = poisson_process(flux1)
+energies = np.arange(.3, 8., .01) * u.keV
+fluxdensity = 6.3e-4 * energies.value**(-1.9) / u.s / u.cm**2 / u.keV
+fluxperbin = fluxdensity[1:] * np.diff(energies)
+flux = poisson_process(fluxperbin.sum())
+energytab = QTable({'energy': energies, 'fluxdensity': fluxdensity})
 
 aperture = chandra.Aperture()
 
-src1 = PointSource(coords=ngc1313_X1, energy={'energy': energies, 'flux': spectrum1},
-                   flux=flux1, geomarea=aperture.area)
+src1 = PointSource(coords=ngc1313_X1, energy=energytab,
+                   flux=flux, geomarea=aperture.area)
 
 pointing = chandra.LissajousDither(coords=ngc1313_X1)
 hrma = chandra.HRMA()
 acis = chandra.ACIS(chips=[4,5,6,7,8,9], aimpoint=chandra.AIMPOINTS['ACIS-S'])
 
-photons = src1.generate_photons(5e3)  # 5 ks exposure time
+photons = src1.generate_photons(5 * u.ks)
 photons = pointing(photons)
 photons = aperture(photons)
 photons = hrma(photons)

docs/pyplots/sourcespectrum.py

@@ -1,26 +1,20 @@
 import numpy as np
 import matplotlib.pyplot as plt
-from astropy.table import Table
+from astropy.table import QTable
+import astropy.units as u
 from marxs.source import Source
-en = np.arange(0.5, 7., .5)
-flux = en**(-2)
+en = np.arange(0.5, 7., .5) * u.keV
+fluxperbin = en.value**(-2) / u.s / u.cm**2 / u.keV
 
-# Input as dictionary
-dictspectrum = {'energy': en, 'flux': flux}
-s1 = Source(energy=dictspectrum, name='dict')
+# Input as astropy QTable
+tablespectrum = QTable([en, fluxperbin], names=['energy', 'fluxdensity'])
+s = Source(energy=tablespectrum, name='table')
 
-# Input as astropy Table
-tablespectrum = Table([en, flux], names=['energy', 'flux'])
-s2 = Source(energy=tablespectrum, name='table')
-
-# Input as numpy arrays
-numpyspectrum = np.vstack([en, flux])
-s3 = Source(energy=numpyspectrum, name='2d array')
 
 fig = plt.figure()
-for s in [s1, s2, s3]:
-    photons = s.generate_photons(1e3)
-    plt.hist(photons['energy'], histtype='step', label=s.name, bins=20)
+
+photons = s.generate_photons(1e3 * u.s)
+plt.hist(photons['energy'], histtype='step', label=s.name, bins=20)
 
 leg = plt.legend()
 lab = plt.xlabel('Energy [keV]')

docs/pyplots/vis_pol.py

@@ -1,6 +1,7 @@
 import numpy as np
 from transforms3d import euler
 from astropy.io import ascii
+import astropy.units as u
 import matplotlib.pyplot as plt
 
 from marxs.source import LabPointSourceCone
@@ -34,7 +35,7 @@
     light.position = np.dot(rotmat, light_pos)
     light.dir = np.dot(rotmat, light_dir)
     m1.geometry.pos4d = np.dot(rotmat, m1pos4d)
-    rays = light(10000)
+    rays = light(10000 * u.s)
     rays = experiment(rays)
     n_detected[i] = rays['probability'].sum()
 

docs/pyplots/vis_subaperturing.py

@@ -46,10 +46,10 @@ def specific_process_photons(self, *args, **kwargs):
 
 instrum = simulator.Sequence(elements=[aper, mirr, gas, det, projectfp])
 star = source.PointSource(coords=SkyCoord(30., 30., unit='deg'),
-                          energy=1., flux=1.)
+                          energy=1. * u.keV, flux=1. / u.s / u.cm**2)
 pointing = source.FixedPointing(coords=SkyCoord(30., 30., unit='deg'))
 
-photons = star.generate_photons(4000)
+photons = star.generate_photons(4000 * u.s)
 photons = pointing(photons)
 photons = instrum(photons)
 ind = (photons['probability'] > 0) & (photons['facet'] >=0)

docs/results.rst

@@ -14,11 +14,12 @@ The most important columns during the simultion are the homogeneous space positi
 Photons from astrophysical sources have the follwoing properties: Time, energy, true sky coordinates (RA and dec) and polarization angle (measures North through east). There is also a probability column (see below), which is initialized to one::
 
   >>> from astropy.coordinates import SkyCoord
+  >>> import astropy.units as u
   >>> from marxs import optics, source
   >>> coords = SkyCoord(345., -35., unit='deg')
-  >>> src = source.PointSource(coords=coords, energy=0.25)
+  >>> src = source.PointSource(coords=coords, energy=0.25 * u.keV)
   >>> point = source.FixedPointing (coords=coords)
-  >>> photons = src.generate_photons(5)
+  >>> photons = src.generate_photons(5 * u.s)
   >>> photons.colnames
   ['time', 'energy', 'polangle', 'probability', 'ra', 'dec']
   >>> photons['probability']

docs/runexample.rst

@@ -20,8 +20,10 @@ In many cases, the spectrum of the sources will come from a table that is genera
 Our source here has a simple power-low spectrum. We ignore extinction for the moment which would significantly reduce the flux we observe in the soft band. 
 
    >>> import numpy as np
-   >>> energies = np.arange(.3, 8., .01)
-   >>> spectrum1 = 6.3e-4 * energies**(-1.9)
+   >>> import astropy.units as u
+   >>> from astropy.table import QTable
+   >>> energies = np.arange(.3, 8., .01) * u.keV
+   >>> spectrum1 = 6.3e-4 * energies.value**(-1.9) / u.s / u.cm**2 / u.keV
 
 Marxs sources have separate parameters for the ``energy`` (shape of the spectrum) and the ``flux`` (normalization of the spectrum and temporal distribution). In many cases we want to apply the same spectrum with different normalizations to different sources, so it is useful to separate the shape of the spectrum and the total flux. We simulate a source that is constant in time, so the photon arrival times will just be Poisson distributed. In this case, the input spectrum is already properly normalized, so we just sum it all up to to get the total flux. Note that MARXS interprets the energy array as the *upper* energy bound for each bin. The lower energy bound is gives as the upper energy bound from the previous bin. This leaves the *lower* bound of the first bin undefined and thus we do not include it in the sum.
 
@@ -37,7 +39,7 @@ The last thing that the source needs to know is how large the geometric opening
 Now, can can finally create our source:
 
    >>> from marxs.source import PointSource
-   >>> src1 = PointSource(coords=ngc1313_X1, energy={'energy': energies, 'flux': spectrum1},
+   >>> src1 = PointSource(coords=ngc1313_X1, energy=QTable({'energy': energies, 'fluxdensity': spectrum1}),
    ...                    flux=flux1, geomarea=aperture.area)
 
 See :ref:`sources` for more details.
@@ -63,7 +65,7 @@ Run the simulation
 ==================
 In MARXS, a source generates a photon table and all following elements process this photon table by changing values in a column or adding new columns. Thus, we can stop and inspect the photons after every step (e.g. to look at their position right after they exit the mirror) or just execute all steps in a row::
 
-   >>> photons = src1.generate_photons(5e3)  # 5 ks exposure time
+   >>> photons = src1.generate_photons(5 * u.ks)
    >>> photons = pointing(photons)
    >>> photons = aperture(photons)
    >>> photons = hrma(photons)

docs/source.rst

@@ -19,7 +19,7 @@ MARXS offers many options to specify the flux, spectrum and polarization that ar
 In addition, we need to give the location of the source and its size and shape (most of the currently implemented sources are point sources, but additional shapes will be added in the future):
 
 - **Astrophysical source**: Needs coordiantes on the sky and a pointing model to translate from sky coordinates to the coordiante system of the satellite. See :ref:`sect-source-radec`. 
-- **Lab source**: Needs x,y,z coordianates in the coordinate system of the experiement as explained in :ref:`pos4d`. Lab source are described in more detail in :ref:`sect-source-lab`.
+- **Lab source**: Needs x,y,z coordianates in the coordinate system of the experiement as explained in :ref:`pos4d`. Lab sources are described in more detail in :ref:`sect-source-lab`.
 
 
  .. _sect-source-fluxenpol:
@@ -31,13 +31,14 @@ The source flux, the energy and the polarization of sources are specified in the
 
 Flux
 ----
-The source flux can just be a number, giving the total counts / second / mm^2 (if no number is given, the default is ``flux=1``).
+The source flux can just be a number with units::
 
      >>> from __future__ import print_function
      >>> from marxs.source import PointSource
      >>> from astropy.coordinates import SkyCoord
-     >>> star = PointSource(coords=SkyCoord("23h12m2.3s -3d4m12.3s"), flux=5.)
-     >>> photons = star.generate_photons(20)
+     >>> import astropy.units as u
+     >>> star = PointSource(coords=SkyCoord("23h12m2.3s -3d4m12.3s"), flux=5. / u.s / u.cm**2)
+     >>> photons = star.generate_photons(20 * u.s)
      >>> photons['time'].format='4.1f'
      >>> print(photons['time'][:6])
      time
@@ -50,27 +51,27 @@ The source flux can just be a number, giving the total counts / second / mm^2 (i
       0.8
       1.0
 
-This will generate 5 counts per second for 20 seconds with an absolutely constant (no Poisson) rate, so the photon list will contain 100 photons.
+This will generate 5 counts per second for 20 seconds with an absolutely constant (no Poisson) rate for a default effective area of :math:`1\;\mathrm{cm}^2`, so the photon list will contain 100 photons.
 
 ``flux`` can also be set to a function that takes the total exposure time as input and returns a list of times, one per photon. In the following example we show how to implement a `Poisson process <https://en.wikipedia.org/wiki/Poisson_process>`_ where the time intervals between two photons are `exponentially distributed <https://en.wikipedia.org/wiki/Poisson_process#Properties>`_ with an average rate of 100 events per second (so the average time difference between events is 0.01 s):
 
     >>> from scipy.stats import expon
     >>> def poisson_rate(exposuretime):
-    ...     times = expon.rvs(scale=0.01, size=exposuretime * 0.01 * 2.)
-    ...     return times[times < exposuretime]
+    ...     times = expon.rvs(scale=0.01, size=exposuretime.to(u.s).value * 0.01 * 2.)
+    ...     return times[times < exposuretime] * u.s
     >>> star = PointSource(coords=SkyCoord(0, 0, unit="deg"), flux=poisson_rate)
 
 Note that this simple implementation is incomplete (it can happen by chance that it does not generate enough photons). MARXS provides a better implementation called `~marxs.source.poisson_process` which will generate the appropriate function automatically given the expected rate:
 
     >>> from marxs.source.source import poisson_process
-    >>> star = PointSource(coords=SkyCoord("23h12m2.3s -3d4m12.3s"), flux=poisson_process(100.))
+    >>> star = PointSource(coords=SkyCoord("23h12m2.3s -3d4m12.3s"), flux=poisson_process(100. / u.s / u.cm**2))
 
 Energy
 ------
-Similarly to the flux, the input for ``energy`` can just be a number, which specifies the energy of a monochromatic source in keV (the default is ``energy=1``):
+Similarly to the flux, the input for ``energy`` can just be a number with a unit, which specifies the energy of a monochromatic source:
 
-    >>> FeKalphaline = PointSource(coords=SkyCoord(255., -33., unit="deg"), energy=6.7)
-    >>> photons = FeKalphaline.generate_photons(5)
+    >>> FeKalphaline = PointSource(coords=SkyCoord(255., -33., unit="deg"), energy=6.7 * u.keV)
+    >>> photons = FeKalphaline.generate_photons(5 * u.s)
     >>> print(photons['energy'])
     energy
      keV
@@ -81,22 +82,22 @@ Similarly to the flux, the input for ``energy`` can just be a number, which spec
        6.7
        6.7
 
-We can also specify a spectrum, by giving binned energy and flux density values. The energy values are taken as the *upper* egde of the bin; the first value of the flux density array is ignored since the lower bound for this bin is undefined. The spectrum can either be in the form of a (2, N) `numpy.ndarray` or it can be some type of table, e.g. an `astropy.table.Table` or a `dict <dict>` with columns named "energy" and "flux" (meaning: "flux density" in counts/s/unit area/keV). In the following exmaple, we specify the same spectrum in three differect ways (the plots look a little different because photon energies are randomly drawn from the spectrum, so there is a Poisson uncertainty):
+We can also specify a spectrum, by giving binned energy and flux density values. The energy values are taken as the *upper* egde of the bin; the first value of the flux density array is ignored since the lower bound for this bin is undefined. The format of the spectrum should be a `~astropy.table.QTable` with columns named "energy" and "fluxdensity" (photons/s/area/keV):
 
 .. plot:: pyplots/sourcespectrum.py
    :include-source:
 
 Two helpful hints:
 
-- If the input spectrum is in some type of file, e.g. fits or ascii, the `astropy.table.Table` `read/write interface <https://astropy.readthedocs.org/en/stable/io/unified.html>`_ offers a convenient way to read it into python::
+- If the input spectrum is in some type of file, e.g. fits or ascii, the `astropy.table.QTable` `read/write interface <https://astropy.readthedocs.org/en/stable/io/unified.html>`_ offers a convenient way to read it into python::
 
-      >>> from astropy.table import Table
-      >>> spectrum = Table.read('AGNspec.dat', format='ascii')  # doctest: +SKIP
+      >>> from astropy.table import QTable
+      >>> spectrum = QTable.read('AGNspec.dat', format='ascii')  # doctest: +SKIP
       >>> agn = PointSource(coords=SkyCoord("11h11m1s -2d3m2.3s", energy=spectrum)  # doctest: +SKIP
 
 - The normalization of the source flux is always is always given by the ``flux`` parameter, independent of the normalizion of ``spectrum['energy']``. If the input is a table and the flux density in that table is in units of photons/s/cm^2/keV, then it is easy to add all that up to set the flux parameter::
 
-    >>> flux = (spectrum['flux'][1:] * np.diff(spectrum['energy'])).sum()   # doctest: +SKIP
+    >>> flux = (spectrum['fluxdensity'][1:] * np.diff(spectrum['energy'])).sum()   # doctest: +SKIP
     >>> agn = PointSource(coords=SkyCoord("11h11m1s -2d3m2.3s", energy=spectrum,
     ...         flux=flux)  # doctest: +SKIP
 
@@ -106,12 +107,13 @@ Lastly, "energy" can be a function that assigns energy values based on the timin
     >>> from marxs.source.source import Source
     >>> import numpy as np
     >>> def time_dependent_energy(t):
+    ...     t = t.value  # convert quantity to plain numpy array
     ...     en = np.ones_like(t)
     ...     en[t <= 5.] = 0.5
     ...     en[t > 5.] = 2.
-    ...     return en
+    ...     return en * u.keV
     >>> mysource = Source(energy=time_dependent_energy)
-    >>> photons = mysource.generate_photons(7)
+    >>> photons = mysource.generate_photons(7 * u.s)
     >>> print(photons['time', 'energy'])
     time energy
      s    keV
@@ -129,11 +131,11 @@ Polarization
 An unpolarized source can be created with ``polarization=None`` (this is also
 the default). In this case, a random polarization is assigned to every
 photon. The other options are very similar to "energy": Allowed are a constant
-angle (in degrees) or a table of some form (see examples above) with two columns "angle" and "probability" (really "probability density") or a numpy array where the first column represents the angle and the second one the probability density. Here is an example where most polarizations are randomly oriented, but an orientation around :math:`35^{\circ}` (0.6 in radian) is a lot more likely.
+angle or a table with columns "angle" and "probabilitydensity". Here is an example where most polarizations are randomly oriented, but an orientation around :math:`35^{\circ}` (0.6 in radian) is a lot more likely.
 
-    >>> angles = np.array([0., 30., 40., 360])
-    >>> prob = np.array([1, 1., 8., 1.])
-    >>> polsource = PointSource(coords=SkyCoord(11., -5.123, unit='deg'), polarization={'angle': angles, 'probability': prob})
+    >>> angles = np.array([0., 30., 40., 360]) * u.degree
+    >>> prob = np.array([1, 1., 8., 1.]) / u.degree
+    >>> polsource = PointSource(coords=SkyCoord(11., -5.123, unit='deg'), polarization=QTable({'angle': angles, 'probabilitydensity': prob}))
 
 Lastly, if polarization is a function, it will be called with time and energy
 as parameters allowing for time and energy dependent polarization
@@ -142,10 +144,10 @@ the 6.4 keV Fe flourescence line after some polarized feature comes into view at
 
     >>> def polfunc(time, energy):
     ...     pol = np.random.uniform(high=2*np.pi, size=len(time))
-    ...     ind = (time > 1000.) & (energy > 6.3) & (energy < 6.5)
+    ...     ind = (time > 1000. * u.s) & (energy > 6.3 * u.keV) & (energy < 6.5 * u.keV)
     ...     # set half of all photons with these conditions to a specific polarization angle
     ...     pol[ind & (np.random.rand(len(time))> 0.5)] = 1.234
-    ...     return pol
+    ...     return pol * u.rad
     >>> polsource = Source(energy=tablespectrum, polarization=polfunc)   # doctest: +SKIP
 
 
@@ -154,7 +156,7 @@ the 6.4 keV Fe flourescence line after some polarized feature comes into view at
 Specify the sky position of an astrophysical source
 ===================================================
 
-An astrophysical source in Marxs must be followed by a pointing model as first optical element that translates the sky coordiantes into the coordinate system of the satellite (see :ref:`pos4d`) and an entrace aperture that selects an initial position for each ray (all rays from astrophysical sources are parallel, thus the position of the source on the sky only determines the direction of a photon but not if it hits the left or the right side of a mirror). See :ref:`sect-apertures` for more details.
+An astrophysical source in marxs must be followed by a pointing model as first optical element that translates the sky coordiantes into the coordinate system of the satellite (see :ref:`pos4d`) and an entrace aperture that selects an initial position for each ray (all rays from astrophysical sources are parallel, thus the position of the source on the sky only determines the direction of a photon but not if it hits the left or the right side of a mirror). See :ref:`sect-apertures` for more details.
 
 
 The following astropysical sources are included in marxs: