Source code for xpsi.Photosphere

from __future__ import division, print_function

from .global_imports import *
from . import global_imports
from . import _warning, make_verbose, verbose

from os.path import join as _join

from .Spacetime import Spacetime
from .HotRegion import HotRegion
from .Elsewhere import Elsewhere
from .Everywhere import Everywhere

from .Parameter import Parameter
from .ParameterSubspace import ParameterSubspace

from .pixelmesh.integrator import integrate as _integrate
from .tools.energy_integrator import energy_integrator
from .tools.phase_integrator import phase_integrator

if _mpl is not None:
    import matplotlib
    from matplotlib import pyplot as plt
    from matplotlib.figure import Figure
    from matplotlib.backends.backend_agg import FigureCanvasAgg as FigureCanvas
    from matplotlib import rcParams
    from matplotlib.ticker import MultipleLocator, AutoLocator, AutoMinorLocator
    from matplotlib import gridspec
    from matplotlib import cm
    from matplotlib import animation
    import matplotlib.image as mgimg
else:
    pass

[docs]class Photosphere(ParameterSubspace): """ A photosphere embedded in an ambient Schwarzschild spacetime. :param obj hot: An instance of :class:`~.HotRegion.HotRegion` (or a derived class). This objects represents the hot regions of the surface that in most use-cases will be assumed to contain radiating material that is hotter than that *elsewhere*. :param obj elsewhere: An instance of :class:`~.Elsewhere.Elsewhere` (or a derived class). :param obj everywhere: An instance of :class:`~.Everywhere.Everywhere` (or a derived class). Note that if you use construct the surface radiation field in this way, you should use the :attr:`~.Photosphere.Photosphere.hot_atmosphere` property to pass a buffer of numerical data to the integrator routines. You then need to ensure that the extension modules ``xpsi/surface_radiation_field/hot_radiation_field.pyx`` and ``xpsi/surface_radiation_field/elsewhere_radiation_field.pyx`` match. .. note:: You cannot specify the surface radiation field *everywhere* if you use hot regions (the latter usage may also include specification of the radiation field *elsewhere*). :param dict bounds: Bounds are supplied for instantiation of a frequency parameter. The parameter name ``'mode_frequency'`` must be a key in the dictionary unless the parameter is *fixed* or *derived*. If a bound is ``None`` that bound is set equal to a strict hard-coded bound. If ``None``, lock the coordinate rotation frequency of a mode of asymmetry in the photosphere to a fixed frequency, e.g., the stellar rotation frequency. If bounds are passed, the frequency is interpreted as a free parameter. :param dict values: Either the fixed value of the mode frequency, a callable if the frequency is *derived*, or a value upon initialisation if the frequency is free. The dictionary must have a key with name ``'mode_frequency'`` if it is *fixed* or *derived*. If the asymmetry is locked to the stellar spin, then you need to pass the spin frequency. If fixed but different to the spin frequency, this value needs to be passed instead. In the hot region base class this mode frequency is applied to normalise the ray lags instead of the stellar rotation frequency. :param iterable custom: A :class:`~.Parameter.Parameter` instance or iterable over such instances. Might be useful for calling image plane extensions and passing global variables, without having to instantiate surface-discretisation classes and without having to handle global variable values at compile time or from disk for runtime access. .. note:: In basic modelling patterns the frequency is the spin frequency, and thus you only need to explicitly pass the spin as ``value`` whilst leaving ``bounds`` to default. If the spin frequency happens to be a free parameter (perhaps with informative prior information), then pass a callable instead that can be used to get the spin frequency dynamically when the derived mode frequency variable is called for. """ required_names = ['mode_frequency'] def __init__(self, hot = None, elsewhere = None, everywhere = None, bounds = None, values = None, custom = None, **kwargs): if everywhere is not None: if hot or elsewhere is not None: raise ValueError('Cannot use hot region nor elsewhere ' 'functionality if constructing the ' 'radiation field everywhere.') if not isinstance(everywhere, Everywhere): raise TypeError('Invalid type for everywhere object.') elif hot is None and elsewhere is None: pass # can call image-plane extensions else: if elsewhere is not None: if not isinstance(elsewhere, Elsewhere): raise TypeError('Invalid type for an elsewhere object.') if hot is None: raise ValueError('Hot region object(s) must be used in ' 'conjuction with an elsewhere object.') self._elsewhere_atmosphere = () # including derived classes if hot is not None and hot is not isinstance(hot, HotRegion): if hasattr(hot, 'objects'): for obj in getattr(hot, 'objects'): if not isinstance(obj, HotRegion): raise TypeError('Invalid object for the hot ' 'region(s).') else: raise TypeError('Invalid object for the hot region(s).') self._hot = hot self._hot_atmosphere = () self._elsewhere = elsewhere self._everywhere = everywhere if bounds is None: bounds = {} if values is None: values = {} doc = """ Coordinate frequency of the mode of radiative asymmetry in the photosphere that is assumed to generate the pulsed signal [Hz]. """ mode_frequency = Parameter('mode_frequency', strict_bounds = (0.0, 2000.0), bounds = bounds.get('mode_frequency', None), doc = doc, symbol = r'$f_{\rm mode}$', value = values.get('mode_frequency', None)) super(Photosphere, self).__init__(mode_frequency, hot, elsewhere, everywhere, custom, **kwargs) @property def hot_atmosphere(self): """ Get the numerical atmosphere buffers for hot regions if used. To preload a numerical atmosphere into a buffer, subclass and overwrite the setter. The underscore attribute set by the setter must be an :math:`n`-tuple whose :math:`n^{th}` element is an :math:`(n-1)`-dimensional array flattened into a one-dimensional :class:`numpy.ndarray`. The first :math:`n-1` elements of the :math:`n`-tuple must each be an ordered one-dimensional :class:`numpy.ndarray` of parameter values for the purpose of multi-dimensional interpolation in the :math:`n^{th}` buffer. The first :math:`n-1` elements must be ordered to match the index arithmetic applied to the :math:`n^{th}` buffer. An example would be ``self._hot_atmosphere = (logT, logg, mu, logE, buf)``, where: ``logT`` is a logarithm of local comoving effective temperature; ``logg`` is a logarithm of effective surface gravity; ``mu`` is the cosine of the angle from the local surface normal; ``logE`` is a logarithm of the photon energy; and ``buf`` is a one-dimensional buffer of intensities of size given by the product of sizes of the first :math:`n-1` tuple elements. It is highly recommended that buffer preloading is used, instead of loading from disk in the customisable radiation field extension module, to avoid reading from disk for every signal (likelihood) evaluation. This can be a non-negligible waste of compute resources. By preloading in Python, the memory is allocated and references to that memory are not in general deleted until a sampling script exits and the kernel stops. The likelihood callback accesses the same memory upon each call without I/O. """ return self._hot_atmosphere @hot_atmosphere.setter def hot_atmosphere(self, path): """ Implement if required. """ raise NotImplementedError('Implement setter if required.') @property def elsewhere_atmosphere(self): """ Get the numerical atmosphere buffers for elsewhere if used. To preload a numerical atmosphere into a buffer, subclass and overwrite the setter. The underscore attribute set by the setter must be an :math:`n`-tuple whose :math:`n^{th}` element is an :math:`(n-1)`-dimensional array flattened into a one-dimensional :class:`numpy.ndarray`. The first :math:`n-1` elements of the :math:`n`-tuple must each be an ordered one-dimensional :class:`numpy.ndarray` of parameter values for the purpose of multi-dimensional interpolation in the :math:`n^{th}` buffer. The first :math:`n-1` elements must be ordered to match the index arithmetic applied to the :math:`n^{th}` buffer. An example would be ``self._hot_atmosphere = (logT, logg, mu, logE, buf)``, where: ``logT`` is a logarithm of local comoving effective temperature; ``logg`` is a logarithm of effective surface gravity; ``mu`` is the cosine of the angle from the local surface normal; ``logE`` is a logarithm of the photon energy; and ``buf`` is a one-dimensional buffer of intensities of size given by the product of sizes of the first :math:`n-1` tuple elements. It is highly recommended that buffer preloading is used, instead of loading from disk in the customisable radiation field extension module, to avoid reading from disk for every signal (likelihood) evaluation. This can be a non-negligible waste of compute resources. By preloading in Python, the memory is allocated and references to that memory are not in general deleted until a sampling script exits and the kernel stops. The likelihood callback accesses the same memory upon each call without I/O. """ return self._elsewhere_atmosphere @elsewhere_atmosphere.setter def elsewhere_atmosphere(self, path): """ Implement if required. """ raise NotImplementedError('Implement setter if required.') @property def hot(self): """ Get the instance of :class:`~.HotRegion.HotRegion`. """ return self._hot @property def elsewhere(self): """ Get the instance of :class:`~.Elsewhere.Elsewhere`. """ return self._elsewhere @property def everywhere(self): """ Get the instance of :class:`~.Everywhere.Everywhere`. """ return self._everywhere @property def spacetime(self): """ Return instance of :class:`~.Spacetime.Spacetime`. """ return self._spacetime @spacetime.setter def spacetime(self, obj): if not isinstance(obj, Spacetime): raise TypeError('Invalid type for spacetime object.') # otherwise store a reference to the spacetime object self._spacetime = obj
[docs] def embed(self, fast_total_counts, threads): """ Embed the photosphere in an ambient Schwarzschild spacetime. In other words, generate a discrete representation of the photospheric radiation field and the null mapping from the photosphere to infinity, for use in flux integrators called by distant observers. """ if self._everywhere is not None: self._everywhere.embed(self._spacetime, self, threads) else: if self._elsewhere is not None: self._elsewhere.embed(self._spacetime, threads) if self._hot is not None: self._hot.embed(self._spacetime, self, fast_total_counts, threads, self._elsewhere._compute_cellParamVecs) elif self._hot is not None: self._hot.embed(self._spacetime, self, fast_total_counts, threads)
[docs] def integrate(self, energies, threads): """ Integrate over the photospheric radiation field. :param energies: A one-dimensional :class:`numpy.ndarray` of energies in keV. :param int threads: Number of ``OpenMP`` threads to spawn for signal integration. """ if self._everywhere is not None: spectrum = self._everywhere.integrate(self._spacetime, energies, threads, self._hot_atmosphere) if spectrum.ndim == 1: self._signal = ((spectrum.reshape(-1,1),),) else: self._signal = ((spectrum,),) else: if self._elsewhere is not None: spectrum = self._elsewhere.integrate(self._spacetime, energies, threads, *self._elsewhere_atmosphere) if self._hot is not None: self._signal = self._hot.integrate(self._spacetime, energies, threads, self._hot_atmosphere, self._elsewhere_atmosphere) if not isinstance(self._signal[0], tuple): self._signal = (self._signal,) # add time-invariant component to first time-dependent component if self._elsewhere is not None: for i in range(self._signal[0][0].shape[1]): self._signal[0][0][:,i] += spectrum
@property def signal(self): """ Get the stored signal. :returns: A tuple of tuples of *ndarray[m,n]*. Here :math:`m` is the number of energies, and :math:`n` is the number of phases. Units are photon/s/keV; the distance is a fast parameter so the areal units are not yet factored in. If the signal is a spectrum because the signal is time-invariant, then :math:`n=1`. """ return self._signal @property def global_variables(self): """ Get a vector of global surface radiation field variables. :returns: An *ndarray[n]* of scalars required to evaluate variables that control the radiation field w.r.t local comoving frames across the stellar surface. The following code block is how one would pass the properties of a single-temperature circular ``HotRegion`` to the extension modules. If you have more than one ``HotRegion`` object merged into the subspace associated with the ``Photosphere`` object, they may each be prefixed, meaning that the set of parameter names below would need to be prefixed at the least, and unless you only want to image one ``HotRegion``, the parameters of the ``HotRegions`` object are required. .. highlight:: python .. code-block:: python return _np.array([self['super_colatitude'], self['phase_shift'] * _2pi, self['super_radius'], self['super_temperature']]) The phase shift controls the initial rotational phase of the ``HotRegion`` when imaging commences. """ try: return _np.array([self['temperature']]) except KeyError: raise NotImplementedError('Subclass and provide an implementation.') @property def global_to_local_file(self): try: return self._global_to_local_file except AttributeError: return None @global_to_local_file.setter def global_to_local_file(self, filepath): if not isinstance(filepath, _six.string_types): raise TypeError('File path must be a string.') elif not _os.path.isfile(filepath): raise IOError('File does not exist.') self._global_to_local_file = filepath @property def images(self): """ Get the precomputed image information. """ return self._images @images.setter def images(self, images): """ Store an *ndarray[i,j,k]* of images and associated information. """ try: for i, obj in enumerate(images): if not isinstance(obj, _np.ndarray): if i < len(images) - 3: raise TypeError('Image information was expected to be ' 'contained in an ndarray.') elif obj is not None and not isinstance(obj, float): raise TypeError('Unexpected type for image information.') except TypeError: raise TypeError('An iterable of objects containing image ' 'information must be supplied.') if len(images) != 13: raise ValueError('There must be six ndarray objects specifing ' 'image information.') msg = 'Image information element %i must have %i dimensions.' # tuple elements: # energy-phase resolved signal (2D array) # x coordinate on image plane (1D array) # y coordinate on image plane (1D array) # colatitude mapped to point (x,y) on image plane (1D array) # azimuth mapped to point (x,y) on image plane (1D array) # radial coord mapped to point (x,y) on image plane (1D array) # phase lag # redshift # aberrated ray angle to local surface normal # elliptical image-plane radial array # elliptical image-plane semi-major axis # elliptical image-plane semi-minor axis # energy-phase resolved specific intensity sky maps (3D array) # the last element is None if intensities not cached assert images[0].ndim == 2, msg % (0, 2) assert images[1].ndim == 1, msg % (1, 1) assert images[2].ndim == 1, msg % (2, 1) assert images[3].ndim == 1, msg % (3, 1) assert images[4].ndim == 1, msg % (4, 1) assert images[5].ndim == 1, msg % (5, 1) assert images[6].ndim == 1, msg % (6, 1) assert images[7].ndim == 1, msg % (7, 1) assert images[8].ndim == 1, msg % (8, 1) assert images[9].ndim == 1, msg % (9, 1) if images[12] is not None: assert images[12].ndim == 3, msg % (12, 3) _num_rays = len(images[1]) for i in range(2,9): assert len(images[i] == _num_rays),\ ('Ray map: array length mismatch (array at tuple index %i is ' 'not equal in length to array at tuple index 1).' % i) assert int( _m.sqrt( _num_rays - 1 ) ) == len(images[9]),\ ('Ray map: array length mismatch for image-plane radial ' 'coordinate array (array at tuple index 9).') if images[12] is not None: assert images[12].shape[0] == images[0].shape[1],\ ('Intensity cache dimension 0 does not match the length of ' 'dimension 1 of the specific flux array ' '(at tuple index 1), meaning the number of phases ' 'is mismatched.') assert images[12].shape[2] == _num_rays,\ ('Intensity cache dimension 2 does not match the length of ' 'ray map arrays (e.g., array at tuple index 1), meaning ' 'the number of rays is mismatched.') self._images = images @images.deleter def images(self): del self._images
[docs] def load_image_data(self, directory): """ Load imaging data from disk. :param str directory: Path to directory to load files from. Should contain files written to disk by :meth:`write_image_data`. """ _d = directory photon_specific_flux = _np.load(_join(_d, 'photon_specific_flux.npy')) x_coordinate = _np.load(_join(_d, 'x_coordinate.npy')) y_coordinate = _np.load(_join(_d, 'y_coordinate.npy')) colatitude = _np.load(_join(_d, 'colatitude.npy')) azimuth = _np.load(_join(_d, 'azimuth.npy')) radial_coord = _np.load(_join(_d, 'radial_coord.npy')) phase_lag = _np.load(_join(_d, 'phase_lag.npy')) redshift = _np.load(_join(_d, 'redshift.npy')) abberated_angle = _np.load(_join(_d, 'abberated_angle.npy')) IP_radial_array = _np.load(_join(_d, 'IP_radial_array.npy')) IP_ellipse_axes = _np.load(_join(_d, 'IP_ellipse_axes.npy')) intensity = _np.load(_join(_d, 'intensity.npy')) self.images = [photon_specific_flux, x_coordinate, y_coordinate, colatitude, azimuth, radial_coord, phase_lag, redshift, abberated_angle, IP_radial_array, IP_ellipse_axes[0], IP_ellipse_axes[1], intensity]
[docs] def write_image_data(self, directory): """ Write imaging data to disk. :param str directory: Path to directory to write to. Must exist. """ _d = directory _np.save(_join(_d, 'photon_specific_flux.npy'), self.images[0]) _np.save(_join(_d, 'x_coordinate.npy'), self.images[1]) _np.save(_join(_d, 'y_coordinate.npy'), self.images[2]) _np.save(_join(_d, 'colatitude.npy'), self.images[3]) _np.save(_join(_d, 'azimuth.npy'), self.images[4]) _np.save(_join(_d, 'radial_coord.npy'), self.images[5]) _np.save(_join(_d, 'phase_lag.npy'), self.images[6]) _np.save(_join(_d, 'redshift.npy'), self.images[7]) _np.save(_join(_d, 'abberated_angle.npy'), self.images[8]) _np.save(_join(_d, 'IP_radial_array.npy'), self.images[9]) _np.save(_join(_d, 'IP_ellipse_axes.npy'), _np.array(self.images[10:12], dtype=_np.double)) _np.save(_join(_d, 'intensity.npy'), self.images[12])
@property def photon_specific_flux(self): """ Get the photon specific flux as a function of phase and energy. :return: A two-dimensional :class:`numpy.ndarray`, where photon energy varies with row number, and phase varies with column number. """ return self._images[0] @property def photon_specific_intensity(self): """ Get the photon specific intensity. Function of phase, energy and sky direction. :return: A three-dimensional :class:`numpy.ndarray`, where the first dimension is phase, the second dimension is photon energy, and the third dimension is sky direction (flattened from two-dimensional sky coordinates to one dimension). """ return self._images[12]
[docs] @make_verbose('Imaging the star', 'Star imaged') def image(self, reimage = False, reuse_ray_map = True, energies = None, num_phases = None, phases = None, phases_in_cycles = False, sqrt_num_rays = 100, epsabs_ray = 1.0e-12, epsrel_ray = 1.0e-12, max_steps = 100000, init_step = 0.1, image_plane_radial_increment_power = 1.0 / 2.0, threads = 1, cache_intensities = False, cache_energy_indices = None, cache_phase_indices = None, single_precision_intensities = True, plot_sky_maps = False, sky_map_kwargs = None, animate_sky_maps = False, free_memory = True, animate_kwargs = None, **kwargs): """ Image the star as a function of phase and energy. :param bool reimage: (Re)image the star. If ``False``, but the spacetime configuration has been updated or the photosphere parameters have been updated, a warning will be generated. In principle, one might want to plot sky maps using cached imaging information, or animate sky maps using images on disk, so reimaging is not forced if (non-fixed) parameters have been changed. :param bool reuse_ray_map: Reuse a precomputed ray map from the stellar surface to the image plane. If the spacetime configuration has changed (non-fixed parameters have changed), a cached ray map will *not* be reused. If the spacetime configuration is unchanged, but resolution settings have changed for ray tracing, pass ``False`` to adhere to the new resolution settings. :param ndarray[n] energies: Energies in keV to evaluate incident specific intensities at. :param int num_phases: The number of phases spanning the unit interval (zero and unity inclusive) to image at. :param ndarray[m] phases: Phases in *radians* or *cycles* at which to evaluate incident specific intensities at. If not ``None``, takes precedence over :obj:`num_phases`. The units need to be specified with the :obj:`phases_in_cycles` keyword argument: if ``False``, give the phase array in *radians*. :param bool phases_in_cycles: Is the phase array, if not ``None``, in units of rotational cycles? :param int sqrt_num_rays: Square-root of the number of rays. This is the level of discretisation in both a radial coordinate and a polar coordinate on an elliptical image plane. .. note:: When the spacetime is static or extremely close to being static in a numerical context, at the resolutions we are interested in, we need to mitigate problems with rays that graze the pole infinitesimally close to the polar axis. In the vicinity of the polar coordinate singularity the ODE system is stiff and the solution is unstable. The most straightforward way to mitigate this is to perform a fallback forward Euler step for a ray that passes exactly through the pole, and use that ray as an approximation for the grazing ray that started very nearby on the image plane. Internally, if a ray intersects the image plane at :math:`x`-coordinate that is numerically very close to, but not exactly, zero (which would mean alignment to the rotational axis), it is approximated by a ray that intersects :math:`x=0`. Image-plane interpolation of quantities (such as intensity) for the purpose of visualisation will then smooth out any such artefacts. Moreover, as an additional measure against artefacts in the sky maps in the vicinity of the rotational pole, rays are distributed accordingingly. For example, if we request :math:`n=400` rays per dimension, a maximal spacing of the rays from the rotational axis is achieved by rotating the *spokes* of rays (by up to :math:`\pm\pi/n`) so that no spoke is aligned (or anti-aligned) with the :math:`y`-direction. :param float epsabs_ray: Absolute error tolerance per ray to adhere to during numerical integration. :param float epsrel_ray: Relative error tolerance per ray to adhere to during numerical integration. :param int max_steps: Maximum number of steps to permit per ray before forced termination of integration. :param float init_step: The initial *suggested* step size at the image plane for the affine parameter for each ray. :param float image_plane_radial_increment_power: Controls the behaviour of the radial discretisation. Higher values towards unity result in linear spacing of rays with the radial coordinate. Lower values towards zero squeeze the rays towards the visible stellar limb, which is necessary for resolving images of extended radiating elements. Values above unity are not recommended, and would squeeze rays towards the image-plane origin, compromising resolution at the stellar limb. :param int threads: Number of OpenMP threads to spawn for parallel blocks of code. Parallel blocks include ray integration to generate a global ray map from image plane to surface; and image calculation at a sequence of rotational phases. :param float cache_intensities: Cache the photon specific intensity sky maps in memory, as a function of phase and energy? The type must be a float (greater than or equal to zero) or ``False``. The value represents the limiting size in GB that can be allocated for the intensity cache. Defaults to zero because this dominates memory consumption. You need to activate this option if you want to plot the sky maps (see below). To activate, supply a limit. A hard limit of 2 GB is imposed for safety. To override, use the secret :obj:`_OVERRIDE_MEM_LIM` keyword argument to supply a positive limit in GB. :param ndarray[m] cache_phase_indices: A one-dimensional :class:`numpy.ndarray` of ``dtype=numpy.int32``, specifying the phase-array indices to cache intensities at. This is useful to save memory when you want to plot specific intensity skymaps but also compute the specific flux at many more phases. If ``None``, intensities will be cached at all phases subject to memory constraints. Note that the order of the list matters for plotting order, so the indices should generally increase, as should the phases themselves. If plotting the pulse-profile and spectrum, then this is a case where many more phases are useful for the resolving specific flux pulse-profile than are needed to plot specific intensity skymaps and specific flux spectra at three representative phases. :param ndarray[m] cache_energy_indices: A one-dimensional :class:`numpy.ndarray` of ``dtype=numpy.int32``, specifying the energy-array indices to cache intensities at. This is useful to save memory when you want to plot specific intensity skymaps but also compute the specific flux at many more energies. If ``None``, intensities will be cached at all energies subject to memory constraints. Note that the order of the list matters for plotting order, so the indices should generally increase, as should the energies themselves. If plotting the pulse-profile and spectrum, then this is a case where many more energies are useful for the resolving specific flux spectrum than are needed to plot specific intensity skymaps and specific flux pulse-profiles at three representative energies. :param bool single_precision_intensities: Cache the intensities in single precision? In most use cases, double precision is simply unnecessary, and because memory consumption can be high, choosing single precision can reduce memory requirements by a factor of two. Note that this only applies to the caching of intensities, not the calculation of intensities, which is done safely in double precision; only the final caching operation is a demotion cast to single precision. The default is single precision caching. Option ignored if intensities are not cached. :param bool plot_sky_maps: Plot (specific) intenity sky maps at a sequence of phases, or by averaging over phase. Maps can be made at one more energies or energy intervals. The images will be written to disk and can be used as frames in an animated sequence. :param dict sky_map_kwargs: Dictionary of keyword arguments passed to :meth:`~Photosphere._plot_sky_maps`. Refer to the associated method docstring for available options. :param bool animate_sky_maps: Compile images from disk into an animated sequence. :param bool free_memory: Try to free the imaging information before animating a sequence of sky maps written to disk, to try to avoid high memory usage. For safety the default is to free the memory, so deactivate this at your own risk. If there are other non-weak references created to the underlying objects, the memory may fail to be freed. In the methods below, the aim is that the native garbage collection cleans up the references because they only exist in the method local scope (no closures or globals). .. note:: Memory used for plotting the sky maps and loading the images from disk to animate a phase sequence might not be straightforwardly freed despite efforts to do so, because of non-weak references covertly held by the matplotlib module. :param dict animate_kwargs: Dictionary of keyword arguments passed to :meth:`~Photosphere._animate`. Refer to the associated method docstring for available options. :param bool deactivate_all_verbosity: Deactivate the verbose output? Note that despite this keyword argument not appearing in the method signature, it is a valid switch. """ ref = self._spacetime # geometry shortcut saves characters try: _DV = deactivate_verbosity except NameError: _DV = False _exc = ValueError('You need to cache intensity sky maps if you ' 'want to plot them.') try: self.images except AttributeError: if not reimage: if plot_sky_maps: raise _exc else: yield ('Warning: star will not be reimaged... assuming ' 'images exist on disk.') else: if not reimage and plot_sky_maps and self.images[-1] is None: raise _exc if phases is not None and not isinstance(phases, _np.ndarray): raise TypeError('Imaging phases must be in a 1D ndarray.') elif isinstance(phases, _np.ndarray): if phases_in_cycles: if phases[0] != 0.0 or phases[-1] != 1.0: _warning('Phase array does not span the unit interval.') phases *= _2pi elif phases is None: if num_phases is None or not isinstance(num_phases, int): raise TypeError('Integer number of phases required.') phases = _np.linspace(0.0, 1.0, num_phases) * _2pi if not isinstance(energies, _np.ndarray): raise TypeError('Imaging energies must be in a 1D ndarray.') time_is_space = sky_map_kwargs.get('time_is_space', False) if reimage: if plot_sky_maps and not cache_intensities: raise _exc if cache_intensities: _override_mem_lim = kwargs.get('_OVERRIDE_MEM_LIM', 1.0) if not isinstance(_override_mem_lim, float): raise TypeError('Intensity cache limit override must be a ' 'float.') elif _override_mem_lim < 0.0: raise ValueError('Intensity cache limit override must be ' 'positive or zero.') if not isinstance(cache_intensities, float): raise TypeError('Intensity cache limit must be a float.') elif not 0.0 <= cache_intensities <= _override_mem_lim: raise ValueError('Intensity cache limit must be positive ' 'and less than the safety limit, which ' 'in turn can be overridden as described ' 'in the method docstring.') if cache_energy_indices is None: cache_energy_indices = _np.arange(len(energies), dtype=_np.int32) elif not isinstance(cache_energy_indices, _np.ndarray): raise TypeError('Energy indices for intensity caching ' 'must be supplied in a 1D numpy.ndarray.') elif cache_energy_indices.dtype != _np.int32: raise TypeError('Energy indices for intensity caching ' 'must be integers.') elif time_is_space and len(cache_energy_indices) != len(energies): raise TypeError('Sky maps must be cached at all energies.') if cache_phase_indices is None: cache_phase_indices = _np.arange(len(phases), dtype=_np.int32) elif not isinstance(cache_phases_indices, _np.ndarray): raise TypeError('Phase indices for intensity caching ' 'must be supplied in a 1D numpy.ndarray.') elif cache_phase_indices.dtype != _np.int32: raise TypeError('Phase indices for intensity caching ' 'must be integers.') elif not time_is_space and len(cache_phases_indices) != len(phases): raise TypeError('Sky maps must be cached at all phases.') _req_size = 4.0 if single_precision_intensities else 8.0 _req_size *= len(cache_phase_indices) * len(cache_energy_indices) # bytes _req_size *= sqrt_num_rays**2.0 # + 1.0 # origin ray negligible if _req_size/1.0e9 >= cache_intensities: raise MemoryError('Too much memory would be required to ' 'cache the intensities at this ' 'resolution. Try decreasing the number ' 'of rays, energies, and/or phases, or ' 'override the cache size limit if ' 'safe.') cache_intensities = True else: cache_intensities = False try: self.images except AttributeError: if reuse_ray_map: yield ('Warning: a ray map has not been cached... ' 'tracing new ray set') else: # if spacetime configuration was updated if ref.needs_update or not reuse_ray_map: # try to free up memory; CPython reference counting means # this should have immediate effect del self.images else: # del self.images[0] # doesn't require much memory del self.images[-1] # requires far more memory try: _ray_map = tuple(self.images[1:]) yield 'Cached ray set to be reused... commencing imaging' except AttributeError: _ray_map = None yield 'Commencing ray tracing and imaging' images = _integrate(threads, ref.r_s, ref.R, ref.Omega, self['mode_frequency'], ref.zeta, ref.epsilon, ref.a, # dimensionless spin ref.q, # mass quadrupole ref.d, ref.i, sqrt_num_rays, epsabs_ray, epsrel_ray, max_steps, init_step, image_plane_radial_increment_power, self.global_variables, energies, phases, cache_intensities, cache_energy_indices, cache_phase_indices, single_precision_intensities, _ray_map, self.global_to_local_file, self._hot_atmosphere) if images[0] == 1: raise Exception('A numerical error arose during imaging ' 'computation... terminating simulation.') elif _ray_map is not None: # only recalculated info is returned # tuple elements: # energy-phase resolved signal (2D array) # energy-phase resolved specific intensity sky maps (3D array) # the last element is None if intensities not cached # transpose so signal phase increments along columns self.images[0] = images[1].T self.images.append(images[2]) else: # the ray map is also returned # tuple elements: # energy-phase resolved signal (2D array) # x coordinate on image plane (1D array) # y coordinate on image plane (1D array) # colatitude mapped to point (x,y) on image plane (1D array) # azimuth mapped to point (x,y) on image plane (1D array) # radial coord mapped to point (x,y) on image plane (1D array) # phase lag # redshift # aberrated ray angle to local surface normal # elliptical image-plane radial array # elliptical image-plane semi-major axis # elliptical image-plane semi-minor axis # energy-phase resolved specific intensity sky maps (3D array) # the last element is None if intensities not cached # transpose so signal phase increments along columns self.images = [images[1].T] + list(images[2:]) yield 'Ray tracing complete.' yield 'Ray set cached.' if cache_intensities: yield 'Intensity caching complete.' else: if len(phases) > 1: yield 'Phase-resolved specific flux integration complete.' else: yield 'Specific flux integration complete.' # memoization self._spacetime([param.value for param in self._spacetime]) if sky_map_kwargs is None: sky_map_kwargs = {} if animate_kwargs is None: animate_kwargs = {} if plot_sky_maps or animate_sky_maps: root_dir = sky_map_kwargs.pop('root_dir', './images') file_root = sky_map_kwargs.pop('file_root', 'skymap') file_root = _os.path.join(root_dir, file_root) phase_average = sky_map_kwargs.get('phase_average', False) if phase_average and time_is_space: raise ValueError('Cannot phase average sky maps when spatial ' 'dimensions are used to render time.') if phase_average and animate_sky_maps: raise ValueError('Phase averaged sky maps cannot be animated.') bolometric = sky_map_kwargs.get('bolometric', False) if bolometric and not time_is_space: raise ValueError('Cannot energy-integrate sky maps when spatial ' 'dimensions are used to render energy.') if bolometric and animate_sky_maps: raise ValueError('Bolometric sky maps cannot be animated.') if plot_sky_maps: if not _os.path.isdir(root_dir): _os.mkdir(root_dir) elif _os.path.isfile(file_root + '_0.png'): yield ('\nWarning: at least one image file exists ' 'in ``%s``.' % root_dir) yield ('Attempting to move image files to a subdirectory ' 'of ``%s``.' % root_dir) try: # to archive the existing image files from datetime import datetime obj = datetime.now() temp = '__datetime__%i.%i.%i__%i.%i.%i' % (obj.day, obj.month, obj.year, obj.hour, obj.minute, obj.second) temp = _os.path.join(root_dir, 'archived_%s' % temp) _os.mkdir(temp) image_files = _os.listdir(root_dir) for image in image_files: if '.png' in image: _os.rename(_os.path.join(root_dir, image), _os.path.join(temp, image)) except Exception as e: raise Exception('Aborting: image files would be ' 'overwritten. %s' % str(e)) else: yield 'Image files archived in subdirectory ``%s``.' % temp figsize, dpi, num_frames = self._plot_sky_maps(file_root, _phases = phases, _energies = energies, _c_idxs = cache_energy_indices, _c_pidxs = cache_phase_indices, _redraw = True, deactivate_verbosity = _DV, **sky_map_kwargs) elif animate_sky_maps: if reimage: raise ValueError('Star was reimaged but sky maps were not ' 'plotted... aborting animation.') figsize, dpi, num_frames = self._plot_sky_maps(file_root, _phases = phases, _energies = energies, _c_idxs = cache_energy_indices, _c_pidxs = cache_phase_indices, _redraw = False, deactivate_verbosity = _DV, **sky_map_kwargs) if animate_sky_maps: if not _os.path.isfile(file_root + '_0.png'): raise IOError('No images located for animation.') if num_frames is None and reimage: if not time_is_space: num_frames = self.images[-1].shape[0] else: num_frames = self.images[-1].shape[1] elif num_frames is None: if not time_is_space: try: num_frames = len(phases) except TypeError: raise TypeError('You need to declare the image phases ' 'in order to include all images from disk.') else: try: num_frames = len(energies) except TypeError: raise TypeError('You need to declare the image energies ' 'in order to include all images from disk.') if free_memory: try: del self.images # try to free up memory except AttributeError: pass self._animate(file_root, num_frames, figsize, dpi, deactivate_verbosity = _DV, **animate_kwargs) yield None
[docs] @make_verbose('Plotting intensity sky maps', 'Intensity sky maps plotted') def _plot_sky_maps(self, _file_root, _phases, _energies, _c_idxs, _c_pidxs, _redraw, threads = 1, with_pulse_profile_and_spectrum = False, time_is_space = False, panel_layout = None, panel_indices = None, cycles = 1, phase_average = False, bolometric = False, energy_bounds = None, phase_bounds = None, num_levels = 100, add_zero_intensity_level = True, normalise_each_panel = True, invert = False, annotate_energies=False, annotate_phases=False, energy_annotation_format='[%.1f keV]', phase_annotation_format='[%.1f cycles]', annotate_location=(0.05,0.05), colormap = None, figsize = (10,10), usetex = False, fontsize_scale = 1.0, tick_spacing = (0.2,1.0), tick_length_scaling = 1.0, dpi_scale = 1.0, **kwargs): """ Helper method for specific intensity sky map visualization. Uses Delaunay triangulation to create an irregular sky mesh and calculate photon (specific) intensity contours at a sequence of phases. Each figure generated contains a sequence of panels arranged in one or two spatial dimensions. Each panel is an intensity sky map, either at a particular energy or integrated over a finite energy interval. Panels cannot mix specific and integrated intensities. Only a some sequence of energy (intervals), in any order, can be identified as labelling panels in one or two spatial dimensions. Time (rotational phase), whilst it could be defined to label a sequence of panels, is only identified as a labelling a sequence of *figures*. Similarly, sequence of energies could be identified as a variable labelling a sequence of figures, but is not. Moreover, energy and time could label panels in to spatial dimensions, but such mixing is not permitted. Finally, variables controlling the source-receiver system could be identified as labels of panels and/or figures, but this is also not supported. More complicated rendering patterns may be supported in future versions, but for now can be achieved via custom extensions building off of the current functionality. :param str _file_root: Relative or absolute path to parent directory for images, extended with the root name for the image files. E.g., the default is ``./images/skymap``. You do not need to change this unless you wish to, and is otherwise reserved for internal use. You may supply a custom file path via keywords ``root_dir`` and ``file_root`` upon calling :meth:`~Photosphere.image`, which are concatenated appropriately. :param ndarray[n] _phases: The phases at which the star was imaged. This is handled internally, so do *not* pass a keyword argument. If phase averaging, the minimum and maximum phases must be zero and :math:`2\pi` radians (i.e., zero and one cycles). :param ndarray[n] _energies: The energies at which the star was imaged. This is handled internally, so do *not* pass a keyword argument. :param ndarray[n] _c_idxs: The energy indices for which the intensity maps were cached for memory-efficieny plotting. This is handled internally, so do *not* pass a keyword argument. :param ndarray[n] _c_pidxs: The energy indices for which the intensity maps were cached for memory-efficieny plotting. This is handled internally, so do *not* pass a keyword argument. :param bool _redraw: Redraw the sky maps? This is handled internally, so do *not* pass a keyword argument. :param int threads: Number of OpenMP threads to spawn. :param bool with_pulse_profile_and_spectrum: A setting that fundamentally changes some behaviours. If deactivated (the default), only photon (specific) intensity skymaps are plotted. The following frame is an example: .. image:: _static/_skymap_plot.png :param bool with_pulse_profile_and_spectrum: If deactivated, the following keyword arguments do not have a use: :obj:`cycles` and :obj:`colormap`. If *activated*, photon specific intensity skymaps at three energies are plotted in each frame, together with their associated photon specific flux pulse-profiles, and also the photon specific flux spectrum at a finer array of energies. Use :obj:`panel_indices` to select the energies. The pulse-profiles are each normalised to their respective maxima, and the spectrum shows the relative orders of magnitude of the specific flux signals. The following frame is an example: .. image:: _static/_skymap_with_pulse_profile_and_spectrum_plot.png :param bool with_pulse_profile_and_spectrum: If activated, a subset of other keyword arguments are ignored: :obj:`energy_bounds`, :obj:`phase_average`, :obj:`panel_layout`, and :obj:`invert`. The panel layout is rigid (not customisable) in order to focus on the plot quality. If more energies were added, the information density in the plot-space might become too high without adding much more new information. :param bool time_is_space: Each image is at constant energy (or is a spectral trace up to that energy) instead of being at constant phase (or a pulse-profile trace up to that phase). :param tuple[int,int] panel_layout: Two elements: the number of rows and columns of panels. If ``None``, a layout is automatically determined based on the number of images to be plotted. :param iterable panel_indices: These ordered integers will be used to select intensity information by indexing the energy dimension of a 3D intensity array. If specific intensites are plotted, these integers should index a subset of energies at which the star was imaged. If intensities are plotted, these integers should index a subset of the energy intervals over which specific intensities are integrated. See the :obj:`energy_bounds` keyword argument. If the flux is calculated at more energies than specific intensities are cached at, then these integers need to index the :obj:`cache_energy_indices` array appropriately. :param int cycles: Nuber of cycles to generate images for. Only relevant if one cycle is different to the next in terms of the frame, e.g., most commonly if plotting the pulse-profile traces over more than one cycle. If frames separated by one cycle are identical, declare the number of cycles to the animator instead. :param bool phase_average: Average each sky map over one revolution (cycle) in phase? Note that the resulting image is incompatible with the currently supported animation mode. The following image is an example: .. image:: _static/_skymap_phaseaveraged.png :param bool bolometric: Integrate each sky map over energy? Note that the resulting image is incompatible with the currently supported animation mode. :param iterable energy_bounds: A set of two-element containers. Each container has an ordered pair of energies which delimit an integration domain. Specific intensity is integrated along each sky direction, at each phase, between these energy bounds. The bounds must be between the minimum and maximum energies at which the star was imaged. If ``None``, specific intensity sky maps will be plotted (the default). This option is ignored if energy is defined as the time dimension, meaning the sky maps are animated with respect to energy. :param iterable phase_bounds: A set of two-element containers. Each container has an ordered pair of energies which delimit an integration domain. Specific intensity is integrated along each sky direction, at each phase, between these energy bounds. The bounds must be between the minimum and maximum energies at which the star was imaged. If ``None``, specific intensity sky maps will be plotted (the default). This option is ignored if phase is defined as the time dimension, meaning the sky maps are animated with respect to phase. .. note:: To use this functionality, the specific intensities must have been cached at all energies the specific flux is calculated at. :param int num_levels: Number of contour levels in (specific) intensity, distributed between minimum finite, and maximum values per panel, or over all panels. See :obj:`normalise_each_panel` keyword argument. :param bool add_zero_intensity_level: Add a contour level at zero intensity such that the colormap minimum corresponds to zero intensity? If ``True`` (the default), then the background sky, where there is by definition zero model intensity, has the same colour only as the subset of the image of the surface that is not radiating in the model. The disadvantage of this choice is that the intensity structure of the image as a function of phase and sky direction is generally not as well- resolved by the colour and greyscale variation. In the limit that the minimum finite intensity of the image is far smaller than the maximum, then the intensity resolution by colour and greyscale values is highest. If ``False``, then the minimum colour is assigned to the minimum finite intensity as a function of phase and sky direction. This also maximally resolves the intensity by colour and greyscale values, which is useful for models wherein the surface radiation field is constructed, for instance, from uniform-temperature localised hot regions. However, in this case the background sky colour is undefined; the background sky colour is thus set to the minimum colour in the colormap, meaning that the fainest subset of the image over phase and sky direction merges with the background sky in terms of colour and greyscale values. :param bool normalise_each_panel: Normalise the contour colormap to each skymap panel uniquely, or globally over all panels? The former yields relative intensity as function of phase and sky direction for an energy or energy interval, whilst the latter offers more spectral information but emission in some panels may not be discernable. :param bool invert: Invert the greyscale to show bright pixels as dark on a white plot background. If a colormap is manually supplied, this just controls the plot background colour. Inversion is recommended for printed format, whilst a black background is more intuitive when in digital format. :param obj colormap: Usage dependent on other settings. If not plotting the pulse-profile and spectrum, then this is simply a (matplotlib) colormap object. Choose something appropriate and *accessible* for a non-negative scalar field (sky intensities). If plotting the pulse-profile and spectrum too, then :obj:`colormap` can be the string ``'RedGreyBlue'`` to invoke the default colour scheme which is reds for the lowest energy intensity skymap and pulse-profile; pure greyscale for the intermediate energy; and blues for the highest energy. Alternatively, if :obj:`colormap` is simply ``None``, the default greyscale will be used for all energies, with all pulse-profiles in black. Lastly, you can supply a three-element list or tuple of colormap objects, ordered from lowest to highest energy; the pulse-profile line colours will be retrieved as the midpoint of the colormap. Note that the background sky colour will be set to the lowest colour in each colourmap. :param tuple(int,int) figsize: The figure size (width, height) in *inches*. If the dimensions are inconsistent with the aspect ratio suggested by the :obj`panel_layout` settings, the height of the figure will be automatically rescaled to achieve congruence, meaning each panel is approximately square. :param bool usetex: Use TeX backend for figure text. :param float fontsize_scale: Use this argument to scale the font size of figure text relative to the default font size that is automatically determined based on the approximate panel size, which is in turn based on the figure size and the panel layout. :param tuple[float,float] tick_spacing: A two-element container. The first element is the minor tick spacing, and the second is the major tick spacing, for both the :math:`x` and :math:`y` directions on the image plane. The units are both the maximum possible angular size of the image of the surface in an ambient Schwarzschild spacetime, :math:`R_{\\rm eq}/\\sqrt{1-r_{\\rm s}/R_{\\rm eq}}`. :param float tick_length_scaling: Use this argument to scale the axis tick lengths relative to the default lengths that are automatically determined based on the panel size. :param float dpi_scale: Use this argument to scale the dots per inch of the figure, relative to the default that is automatically determined based on the panel size. """ file_root = _file_root phases = _phases energies = _energies redraw = _redraw if with_pulse_profile_and_spectrum: if len(panel_indices) != 3: raise ValueError('Selected plot type designed for showcasing ' 'the specific photon intensity skymaps and ' 'their associated specific photon flux ' 'pulse-profiles or spectra specifically at ' 'at three energies or phases to avoid ' 'excessive information density.') panel_layout = (2,3) if not isinstance(cycles, int): raise TypeError('Declare the number of cycles with an integer.') elif cycles < 1: cycles = 1 # quietly ignore input if not time_is_space: num_frames = len(phases) else: num_frames = len(energies) elif cycles > 1: if not time_is_space: num_frames = len(phases) + (cycles - 1) * (len(phases) - 1) else: num_frames = len(energies) elif panel_layout is None: x = int(_m.ceil(_m.sqrt(len(panel_indices)))) if x * (x - 1) >= len(panel_indices): panel_layout = (x, x - 1) else: panel_layout = (x, x) if not time_is_space: num_frames = len(phases) else: num_frames = len(energies) # try to improve the aspect ratio so that each panel is # approximately square width = panel_layout[1] + (panel_layout[1] - 1) * 0.2 _hspace = 0.25 if with_pulse_profile_and_spectrum else 0.2 height = panel_layout[0] + (panel_layout[0] - 1) * _hspace aspect_ratio = height/float(width) if aspect_ratio != 1.0: if _np.abs(figsize[1]/float(figsize[0]) - aspect_ratio)/aspect_ratio > 0.1: figsize = (figsize[0], figsize[0] * aspect_ratio) # calculate an appropriate dpi to resolve each panel adequately dpi = (max(panel_layout) / 2.0) * 150.0 * dpi_scale if redraw: rcParams['text.usetex'] = usetex try: iter(panel_indices) except TypeError: raise TypeError('Panel indices object must be iterable.') if _np.product(panel_layout) < len(panel_indices): raise ValueError('There are too few panels for the requested ' 'number of intensity sky maps.') # some scaling for appropriate fontsize panel_size = max(figsize[1]/float(panel_layout[0]), figsize[0]/float(panel_layout[1])) fontsize = (panel_size/5.0) * 14.0 * fontsize_scale rcParams['font.size'] = int(fontsize) tick_length = int((panel_size/5.0) * 8 * tick_length_scaling) # get coordinates of irregular set of points for triangulation X = self.images[1] Y = self.images[2] #if not isinstance(energies, _np.ndarray): # raise TypeError('Imaging energies must be in an ndarray.') #if not isinstance(phases, _np.ndarray): # raise TypeError('Imaging phases must be in an ndarray.') images = self.images[-1] if time_is_space: # transpose dimensions _images_ = _np.zeros((images.shape[1], images.shape[0], images.shape[-1]), dtype=_np.double) for i in range(images.shape[-1]): _images_[:,:,i] = images[:,:,i].T images = _images_ if with_pulse_profile_and_spectrum: if not time_is_space: flux = self.images[0] else: flux = self.images[0].T if not with_pulse_profile_and_spectrum and not time_is_space and energy_bounds: with verbose(True, 'Integrating specific intensity over energy intervals', 'Integrated specific intensity over energy intervals'): for bounds in energy_bounds: if bounds[0] > bounds[1]: raise ValueError('Energy bounds in a tuple must be ' 'ordered.') for bound in bounds: if not energies[0] <= bound <= energies[-1]: raise ValueError('Extrapolation would be required.') if len(panel_indices) < len(energy_bounds): yield 'Warning: fewer panels than energy intervals.' integrated = _np.zeros((images.shape[0], len(energy_bounds), images.shape[2]), dtype=_np.double) intensities = _np.zeros((images.shape[1], images.shape[2]), dtype=_np.double) for i in range(images.shape[0]): # phases intensities[...] = images[i,...] # sky directions for k in range(len(energy_bounds)): bounds = _np.log10( _np.array(energy_bounds[k]) ) _integrated = energy_integrator(threads, intensities, _np.log10(energies), bounds) integrated[i,k,:] = _integrated[0,:] images = integrated elif not with_pulse_profile_and_spectrum and not time_is_space: if len(panel_indices) != len(_c_idxs): yield ('Warning: number of panels not equal to number of ' 'phases.') elif not with_pulse_profile_and_spectrum and time_is_space: if len(panel_indices) != len(_c_pidxs): yield ('Warning: number of panels not equal to number of ' 'phases.') if not with_pulse_profile_and_spectrum and not time_is_space and phase_average: with verbose(True, 'Averaging (specific) intensity over rotational phase', 'Averaged (specific) intensity over rotational phase'): if phases[0] != 0.0 or phases[-1] != _2pi: raise ValueError('Minimum and maximum phases at which ' 'star is imaged must be zero and unity ' 'if you are phase averaging.') averaged = _np.zeros((1, images.shape[1], images.shape[2]), dtype = _np.double) intensities = _np.zeros((images.shape[2], images.shape[0]), dtype = _np.double) for i in range(images.shape[1]): # energies for j in range(images.shape[2]): # sky directions intensities[j,:] = images[:,i,j] _averaged = phase_integrator(1.0, # exposure time _np.array([0.0, 1.0]), intensities, phases / _2pi, 0.0) # phase shift for j in range(images.shape[2]): averaged[:,i,j] = _averaged[j,:] images = averaged if not with_pulse_profile_and_spectrum and time_is_space and phase_bounds: with verbose(True, 'Integrating specific intensity over phase intervals', 'Integrated specific intensity over phase intervals'): for bounds in phase_bounds: if bounds[0] > bounds[1]: raise ValueError('Phase bounds in a tuple must be ' 'ordered.') for bound in bounds: if not phases[0] <= bound <= phases[-1]: raise ValueError('Extrapolation would be required.') if len(panel_indices) < len(phase_bounds): yield 'Warning: fewer panels than phase intervals.' integrated = _np.zeros((images.shape[0], len(phase_bounds), images.shape[2]), dtype=_np.double) intensities = _np.zeros((images.shape[2], images.shape[1]), dtype=_np.double) for i in range(images.shape[0]): # energies intensities[...] = images[i,...].T # sky directions for k in range(len(phase_bounds)): bounds = _np.array(phase_bounds[k]) _integrated = phase_integrator(1.0, bounds, intensities, phases / _2pi, 0.0) integrated[i,k,:] = _integrated[:,0] images = integrated if not with_pulse_profile_and_spectrum and time_is_space and bolometric: with verbose(True, 'Integrating bolometric intensity', 'Averaged bolometric intensity'): if phases[0] != 0.0 or phases[-1] != _2pi: raise ValueError('Minimum and maximum phases at which ' 'star is imaged must be zero and unity ' 'if you are phase averaging.') integrated = _np.zeros((images.shape[0], 1, images.shape[2]), dtype = _np.double) intensities = _np.zeros((images.shape[1], images.shape[2]), dtype = _np.double) for i in range(images.shape[0]): # phases for j in range(images.shape[2]): # sky directions intensities[:,j] = images[i,:,j] bounds = _np.log10( _np.array([energies[0], energies[-1]]) ) _integrated = energy_integrator(threads, intensities, _np.log10(energies), bounds) for j in range(images.shape[2]): integrated[i,:,j] = _integrated[:,j] images = integrated if normalise_each_panel: with verbose(True, 'Normalising each sky map panel separately', 'Normalised sky map panels separately'): # normalise intensity for each individual panel levels = [] if not time_is_space: for j in range(images.shape[1]): # at each energy # find extreme intensities over discrete set of image # phases and sky directions MIN = _np.min(images[:,j,:][images[:,j,:] > 0.0]) MAX = _np.max(images[:,j,:]) levels.append(_np.linspace(MIN, MAX, num_levels)) if add_zero_intensity_level: levels[-1] = _np.array([0.0, 0.001*MIN] + list(levels[-1])) else: for j in range(images.shape[0]): # at each energy # find extreme intensities over discrete set of image # phases and sky directions MIN = _np.min(images[j,:,:][images[j,:,:] > 0.0]) MAX = _np.max(images[j,:,:]) levels.append(_np.linspace(MIN, MAX, num_levels)) if add_zero_intensity_level: levels[-1] = _np.array([0.0, 0.001*MIN] + list(levels[-1])) else: with verbose(True, 'Normalising sky map panels globally' 'Normalised sky map panels globally'): MIN = _np.min(images[:,:,:][images[:,:,:] > 0.0]) MAX = _np.max(images[:,:,:]) levels = _np.linspace(MIN, MAX, num_levels) if add_zero_intensity_level: levels = _np.array([0.0, 0.001*MIN] + list(levels)) # because of default tick formatting and a minus sign, # the left and bottom margins need to be different left = 0.09 * (fontsize/14.0) bottom = 0.11 * (fontsize/14.0) right = 0.975 top = bottom + (right - left) ref = self._spacetime fig = Figure(figsize = figsize) canvas = FigureCanvas(fig) if with_pulse_profile_and_spectrum: if not time_is_space and colormap == 'RedGreyBlue': cmap = [cm.Reds_r, cm.Greys_r, cm.Blues_r] _line_colors = [cm.Reds_r(0.25), cm.Greys_r(0.0), cm.Blues_r(0.25)] elif not isinstance(colormap, (list, tuple)): cmap = [cm.Greys_r] * 3 _line_colors = [cm.Greys_r(0.0)] * 3 else: cmap = colormap _line_colors = [cmap[j](0.5) for j in range(len(cmap))] gs = gridspec.GridSpec(panel_layout[0], panel_layout[1], left=left, right=right, bottom=bottom, top=top, wspace=0.2, hspace=_hspace) axes = [fig.add_subplot(gs[j]) for j in range(len(panel_indices))] pp_ax = fig.add_subplot(gs[len(panel_indices):-1]) spec_ax = fig.add_subplot(gs[-1]) line_styles = ['-', '--', '-.'] else: cmap = colormap or (cm.Greys if invert else cm.Greys_r) gs = gridspec.GridSpec(panel_layout[0], panel_layout[1], left=left, right=right, bottom=bottom, top=top, wspace=0.2, hspace=_hspace) axes = [fig.add_subplot(gs[j]) for j in range(len(panel_indices))] if with_pulse_profile_and_spectrum: if not time_is_space: annotate_energies = True else: annotate_phases = True _I = 10 for i in range(images.shape[0]): if phase_average: yield 'Rendering phase-averaged images' elif bolometric: yield 'Rendering bolometric images' elif i == 0 and images.shape[0] < 10: yield 'Rendering images' elif i == 0 and images.shape[0] >= 10: yield 'Rendering image numbers [%i, %i]'%(i+1, i+_I) elif i%_I == 0: yield 'Rendering image numbers (%i, %i]'%(i, i+_I) for j, idx in enumerate(panel_indices): ax = axes[j] if with_pulse_profile_and_spectrum: _cmap = cmap[j] else: _cmap = cmap if (with_pulse_profile_and_spectrum or _np.product(panel_layout) - j - 1 < panel_layout[1]): if ref.R < 1.5 * ref.r_s: ax.set_xlabel(r'$(2x/(3\sqrt{3}r_{\rm s}))$') else: ax.set_xlabel(r'$(x/R_{\rm eq})\sqrt{1-r_{\rm s}/R_{\rm eq}}$') if j % panel_layout[1] == 0: if ref.R < 1.5 * ref.r_s: ax.set_ylabel(r'$(2y/(3\sqrt{3}r_{\rm s}))$') else: ax.set_ylabel(r'$(y/R_{\rm eq})\sqrt{1-r_{\rm s}/R_{\rm eq}}$') _veneer(tick_spacing, tick_spacing, ax, length = tick_length) if with_pulse_profile_and_spectrum: ax.set_facecolor(_cmap(0.0)) else: ax.set_facecolor('white' if invert else 'black') if not time_is_space: lvls = levels if isinstance(levels, _np.ndarray) else levels[idx] else: lvls = levels if isinstance(levels, _np.ndarray) else levels[i] ax.tricontourf(X, Y, images[i,idx,:], cmap = _cmap, levels = lvls) # correct the aspect ratio x_view = ax.xaxis.get_view_interval() diff = x_view[1] - x_view[0] ax.xaxis.set_view_interval(x_view[0] - diff * 0.025, x_view[1] + diff * 0.025) y_view = ax.yaxis.get_view_interval() ax.yaxis.set_view_interval(y_view[1] - diff * 1.025, y_view[1] + diff * 0.025) if not time_is_space and annotate_energies: ax.text(annotate_location[0], annotate_location[1], s=energy_annotation_format % energies[_c_idxs[idx]], fontdict={'color': 'black' if invert else 'white'}, transform=ax.transAxes) if time_is_space and annotate_phases: ax.text(annotate_location[0], annotate_location[1], s=phase_annotation_format % (phases[_c_pidxs[idx]]/_2pi), fontdict={'color': 'black' if invert else 'white'}, transform=ax.transAxes) if with_pulse_profile_and_spectrum: # plot summaries _upper_view_lim = _np.max(flux) * 100.0 _lower_view_lim = _np.min(flux) _view_lim_diff = _np.log10(_upper_view_lim) _view_lim_diff -= _np.log10(_lower_view_lim) for j, idx in enumerate(panel_indices): _idx = _c_idxs[idx] _diff = _np.log10(flux[_idx,i]) _diff -= _np.log10(_lower_view_lim) spec_ax.axvline(energies[_idx], 0.0, _diff/_view_lim_diff, color=_line_colors[j], ls=line_styles[j], lw=1.0) spec_ax.set_xscale('log') spec_ax.set_yscale('log') spec_ax.set_xlim(energies[0], energies[-1]) spec_ax.set_ylim(_lower_view_lim, _upper_view_lim) _veneer(None, None, spec_ax, length = tick_length, log=(True, True)) spec_ax.set_xlabel('Photon energy [keV]') if not time_is_space: spec_ax.plot(energies, flux[:,i], 'k-', lw=1.0) else: for j, idx in enumerate(panel_indices): _idx = _c_pidxs[idx] spec_ax.plot(energies, flux[_idx,:i], color=_line_colors[j], ls=line_styles[j], lw=1.0, label=phase_annotation_format % phases[_idx]) if not time_is_space: for _cycle in range(cycles): # plot pulse-profiles _ext_phases = [] for _i in range(_cycle): _ext_phases += list(phases[1 if _i > 0 else 0:] + _i * _2pi) _ext_phases += list(phases[1 if _cycle > 0 else 0:i+1] + _cycle * _2pi) _ext_phases = _np.array(_ext_phases) for j, idx in enumerate(panel_indices): _idx = _c_idxs[idx] _max = _np.max(flux[_idx,:]) _ext_flux = [] for _i in range(_cycle): _ext_flux += list(flux[_idx, 1 if _i > 0 else 0:]) _ext_flux += list(flux[_idx, 1 if _cycle > 0 else 0:i+1]) _ext_flux = _np.array(_ext_flux) pp_ax.plot(_ext_phases/_2pi, _ext_flux/_max, color=_line_colors[j], ls=line_styles[j], lw=1.0, label=energy_annotation_format % energies[_idx]) pp_ax.set_xlim(0.0, float(cycles)) pp_ax.set_ylim(0.0,1.2) _veneer((0.1,0.5), (0.05,0.2), pp_ax, length = tick_length) pp_ax.legend(loc='upper center', ncol=3, mode='expand', handlelength=4.0, frameon=False, fancybox=False) pp_ax.set_xlabel('Phase [$2\pi$ radians]') pp_ax.set_ylabel('photons/cm$^2$/s/keV') fig.savefig(file_root + '_%i.png' % (len(_ext_phases) - 1), dpi=dpi) pp_ax.clear() spec_ax.clear() else: pp_ax.plot(phases/_2pi, flux[i,:]/_np.max(flux[i,:]), 'k-', lw=1.0) pp_ax.set_xlim(0.0, 1.0) pp_ax.set_ylim(0.0,1.2) _veneer((0.1,0.5), (0.05,0.2), pp_ax, length = tick_length) pp_ax.legend(loc='upper center', ncol=3, mode='expand', handlelength=4.0, frameon=False, fancybox=False) pp_ax.set_xlabel('Phase [$2\pi$ radians]') pp_ax.set_ylabel('photons/cm$^2$/s/keV') fig.savefig(file_root + '_%i.png' % i, dpi=dpi) else: fig.savefig(file_root + '_%i.png' % i, dpi=dpi) for ax in axes: ax.clear() for ax in axes: ax.cla() plt.close(fig) yield figsize, dpi, num_frames
[docs] @staticmethod @make_verbose('Animating intensity sky maps', 'Intensity sky maps animated') def _animate(_file_root, _num_frames, _figsize, _dpi, cycles = 1, fps = None, **kwargs): """ Helper method to animate the intensity sky maps. :param str _file_root: Reserved for internal use. :param int _num_frames: Reserved for internal use. :param int _figsize: Reserved for internal use. :param int _dpi: Reserved for internal use. :param int cycles: Number of explicit cycles to generate frames for. The frames from the principal cycle are reused, so the images are periodic in their phase evolution. There is no delay between cycles in this type of repitition. :param int fps: Frames per second. If ``None``, then one cycle (assuming images have been precomputed for a complete cycle), consisting of so many frames, will exhibit a period of one second. :param bool repeat: Inform *ffmpeg* to enter a loop when video play back commences. Deprecated. :param repeat_delay: Delay between repeats in milliseconds. Deprecated. :param str ffmpeg_path: Absolute path to *ffmpeg* executable. If ``None``, defaults to matplotlib rcParams settings, but no guarantee that the package will be found even if installed on system. Deprecated. """ file_root = _file_root num_frames = _num_frames figsize = _figsize dpi = _dpi fig = plt.figure(figsize = figsize) ax = plt.subplot(111) plt.axis('off') fig.subplots_adjust(left=0, bottom=0, right=1, top=1, wspace=None, hspace=None) # animation code based on: # http://jakevdp.github.io/blog/2012/08/18/matplotlib-animation-tutorial/ filename = file_root + '_0.png' img = mgimg.imread(filename) imgplot = ax.imshow(img, aspect='auto') cycles = int(cycles) if cycles < 1: cycles = 1 elif cycles > 1: num_frames += (cycles - 1) * (_num_frames - 1) class _context: # mutable nonlocal namespace in closure _cycle_idx = 0 # track cycle _j = -1 # track image index _last = -1 def _update(i): # load phase-ordered set of images if _context._last == i: return imgplot, if _context._cycle_idx == 0 and i == _num_frames: _context._j = 1 _context._cycle_idx = 1 elif i == _num_frames + _context._cycle_idx * (_num_frames - 1): _context._j = 1 _context._cycle_idx += 1 else: _context._j += 1 _context._last = i filename = file_root + '_%i.png' % _context._j img = mgimg.imread(filename) imgplot.set_data(img) return imgplot, ani = animation.FuncAnimation(fig, _update, frames=num_frames, blit=True) if fps is None: fps = num_frames # all frames span one second # secret keyword argument; not clear whether should be exposed to user bitrate = kwargs.get('bitrate', -1) # default is let _mpl choose filename = file_root + '_animated.mp4' yield 'Writing to disk: %s' % filename ani.save(filename, writer = 'ffmpeg', dpi = dpi, fps = fps, bitrate = bitrate, extra_args=['-vcodec', 'libx264']) fig.clf() # this or ax.cla() needed to free memory plt.close(fig) yield None
Photosphere._update_doc() def _veneer(x, y, axes, lw=1.0, length=8, log=(False, False)): """ Make the plots a little more aesthetically pleasing. """ if x is not None: if x[1] is not None: axes.xaxis.set_major_locator(MultipleLocator(x[1])) if x[0] is not None: axes.xaxis.set_minor_locator(MultipleLocator(x[0])) elif not log[0]: axes.xaxis.set_major_locator(AutoLocator()) axes.xaxis.set_minor_locator(AutoMinorLocator()) if y is not None: if y[1] is not None: axes.yaxis.set_major_locator(MultipleLocator(y[1])) if y[0] is not None: axes.yaxis.set_minor_locator(MultipleLocator(y[0])) elif not log[1]: axes.yaxis.set_major_locator(AutoLocator()) axes.yaxis.set_minor_locator(AutoMinorLocator()) axes.tick_params(which='major', colors='black', length=length, width=lw) axes.tick_params(which='minor', colors='black', length=int(length/2), width=lw) plt.setp(axes.spines.values(), linewidth=lw, color='black')