.. _surface_brightness_map-label: Computing surface brightness maps ================================= .. currentmodule:: pNbody.Mockimgs If an N-body model ``snapshot.hdf5`` contains all necessary informations (see :ref:`mockimgs_preparation-label`), the ``Mockimgs`` modules provides a script to directly compute surface brightness maps as a ``.fits`` file:: mockimgs_sb_compute_images snapshot.hdf5 --instrument arrakihs_vis_G --distance 35 --rsp_mode const --rsp_val 1.5 -o output.fits See :ref:`mockimgs_scripts-label` for the full documentation. It is possible to save several maps obtained from different line of sight using the ``--nlos`` option:: mockimgs_sb_compute_images snapshot.hdf5 --instrument arrakihs_vis_G --distance 35 --rsp_mode const --rsp_val 1.5 --nlos 4 -o output.fits Sightlines ---------- There are different ways to define one or multiple lines-of-sight. * by default the observer is along the z axis looking at the origin. * one specific line-of-sight is defined with the option:: --los 0 0 1 where the three number gives the position of the observer looking at the origin. * using:: --nlos 9 --random_los will use 9 randomly defined lines-of-sight. A random seed may be defined with the option ``--random_seed`` * if only ``--nlos`` is provided, different lines-of-sight are generated by splitting an hemisphere equally. Note that the effective number of line of sights ``nlos`` will be different than the integer provided, as the procedure guarantee that the number of division in phi is twice the one in theta and keep only one polar ligne of sight. * if a yaml parameter file is provided with the ``-p`` option, additional parameter may be provided. With the following example:: LineOfSights: los : [0,1,0] # line of sight (position of the observer looking at [0,0,0]) grid: n_phi : 3 # number of divisions in azimuthal angle n_theta : 3 # number of divisions in elevation angle d_phi : 10 # azimuthal variation in degree [-d_phi,-d_phi] d_theta : 10 # elevation variation in degree [-d_theta,-d_theta] 9 sightlines will be used, centered on [0,1,0] and slightly rotated by 10 degree (covering a 3x3 grid in angle space). List of instruments ------------------- The list of available implemented instruments is obtained with the ``--instruments_list`` command:: mockimgs_sb_compute_images --instruments_list Multiple outputs ---------------- If multiples lines-of-sight are used, a similar number of fits files will be generated and named lines-of-sight ``output_LOS0i.fits``, with ``i`` from ``0`` to the number of ``-1``. Concatenated outputs -------------------- Instead of saving muliple fits files, it is possible to dump maps in a ``.pkl`` file:: mockimgs_sb_compute_images snapshot.hdf5 --instrument arrakihs_vis_G --distance 35 --rsp_mode const --rsp_val 1.5 --nlos 4 -o output.pkl .. warning:: In ``.pkl`` files stored images are not surface brightness images but fluxes. To convert into a surface brightness image, the following procedure can be used, where ``images`` is a list of images contained in the file:: # split into a matrix and parameter FluxMap,params = images[i] # apply some filters FluxMap = libMockimgs.FilterFluxMap(FluxMap,None) # compute the surface brightness map from the luminosity/flux map SBmap = libMockimgs.FluxToSurfaceBrightness(FluxMap,params["object_distance"],params["ccd_pixel_area"]) The ``.pkl`` file can be displayed with:: mockimgs_sb_show_images output.pkl The list of pre-defined instruments is given here :ref:`mockimgs_instruments-label`. Defining an instrument ---------------------- It is possible to provide ``mockimgs_sb_compute_images`` with an instrument which is not pre-defined. To do so, just define the properties of an instrument in a text file (e.g. ``instrument.py``) with the following content like:: instrument = instrument.Instrument( name = "myInstrument", telescope = telescope.Telescope(name="iSIM-170",focal=1500*u.mm), ccd = ccd.CCD(name="arrakihs_vis",shape=[4096,4096],size=[41*u.mm,41*u.mm]), filter_type = filters.Filter("GAEA_J"), ) as explained in :ref:`mockimgs-label`. To create the suface brightness map, provide the name of your file to the option ``--instrument``:: mockimgs_sb_compute_images snapshot.hdf5 --instrument instrument.py --distance 35 --rsp_mode const --rsp_val 1.5 -o output.fits Note that is it convenient to check the properties of the instrument using the ``--info`` option:: mockimgs_sb_compute_images --instrument instrument.py --info Surface brigthness smoothing ---------------------------- To get more realistic surface brigthness images, the flux of each particle is spread among several pixels at projection time using a 2D gaussian of fwmh of :math:`r_h`. There are several parameters we can use to influence this smoothing. .. list-table:: Smoothing modes :widths: 50 50 50 :header-rows: 1 * - mode option - additional options - description * - ``--rsp_mode=None`` - ``--rsp_fac=1`` - the softening radius of each particles (stored in ``nb.rsp``) is multiplied by the value provided by ``--rsp_fac``. * - ``--rsp_mode=const`` - ``--rsp_val=1.5`` - the softening radius of each particles is fixed to value provided by ``--rsp_fac``. * - ``--rsp_mode=arctan`` - ``--rsp_max=5 --rsp_sca=3`` - the softening radius of each particles is set to :math:`s_{\max}/(\pi/2) \arctan(r_h/s)`, where :math:`r_h=\rm{np.rsp}`, :math:`s_{\max}` is given by ``--rsp_max`` and :math:`s` by ``--rsp_sca``. * - ``--rsp_mode=ln`` - ``--rsp_sca=3`` - the softening radius of each particles is set to :math:`\log(r_h/s+1)`, where :math:`r_h=\rm{np.rsp}` and :math:`s` is given by ``--rsp_sca``. Store the IDs of projected particles ------------------------------------ In certain cases, it can be useful to store the IDs of particles that fall within a given pixel during projection. This can, for example, help trace back the 3D origin of a structure identified on a 2D surface brightness map. The ``mockimgs_sb_compute_images`` script offers this possibility via the option ``--IDsMap``. If this option is added and supplemented with a filename as for example:: mockimgs_sb_compute_images snapshot.hdf5 --instrument arrakihs_vis_G --distance 35 --IDsMap IDsMap.npy the file ``IDsMap.npy`` (numpy format) is created and stors the IDs of particles according to the pixel they fall-in. The following code demonstrates how to retrieve the list of particle IDs that fall within the central pixel of the image, defined by a pixel radius of 10:: import numpy as np mat = np.load('IDsMap.npy',allow_pickle=True) # define a grid nx = mat.shape[0] ny = mat.shape[1] x = np.arange(nx) y = np.arange(ny) x,y = np.meshgrid(x,y) # create a mask rpix = 10 mask = np.sqrt((x-xc)**2 + (y-yc)**2) < rpix # ravel the data x = x.ravel() y = y.ravel() mask = mask.ravel() # keep only pixels in the mask x = np.compress(mask,x) y = np.compress(mask,y) # this will become the list of IDs IDs = np.array([]) for k in range(len(x)): ix = int(x[k]) iy = int(y[k]) if mat[ix,iy] is None: continue IDs = np.append(IDs,mat[ix,iy]) This second example shows how to create a new matrix that contains the number of particles that fall in each pixel:: import numpy as np mat = np.load('IDsMap.npy',allow_pickle=True) matn = np.zeros(mat.shape) for ix in range(mat.shape[0]): for iy in range(mat.shape[1]): if mat[ix,iy] is None: matn[ix,iy] = 0 else: matn[ix,iy] = len(mat[ix,iy]) matn = matn.astype(float) Extracting particles from a snapshot using a 2D mask ---------------------------------------------------- The list of IDs attached to the pixels of a surface brightness image can be used to extract from a 3D model, any feature detected in 2D in this image. Lets start with a snapshot ``snapshot.hdf5`` from which we created a surface brightness images ``snapshot.fits.gz`` using ``mockimgs_sb_compute_images`` as well as an ID file ``snapshot.npy`` using the ``--IDsMap`` option of the ``mockimgs_sb_compute_images``. Imaging now that we created masks on 2D features identified on the ``snapshot.fits.gz`` image. Those mask serves to tag features. They are saved on a new fits file ``masks.fits.gz`` with pixels belonging to a features having a value that corresponds to the ID of the feature. To extract the particles from ``snapshot.hdf5`` that fall on a mask, and are responsible to creating the feature, we can use the following command:: mockimgs_sb_extract_3Dfeatures_from_2Dmask snapshot.npy --masks-file masks.fits.gz -n 1 --nb-file snapshot.hdf5 -o snapshot-LSB1.hdf5 The output file ``snapshot-LSB1.hdf5`` will contains only the extracted particles. Detail of surface brightness computation ---------------------------------------- Analytical formulae ................... We first define the following quantities: * :math:`L` is the luminosity in :math:`\rm{erg/s}` * :math:`d` is the distance in unit of :math:`10\,\rm{pc}` * :math:`A` is the pixel area in :math:`arcsec^2` * :math:`M_0` is the magnitude zero point. The absolute magnitude :math:`M` is related to the luminosity by: .. math:: M = M_0 -2.5 \log_{10}(L) The apparent magnitude :math:`m` is then: .. math:: m = M + 5 \log_{10}(d) The surface brightness magnitude is defined as .. math:: SB = m + 2.5 \log_{10}(A) or equivalently .. math:: SB = M_0 -2.5 \log_{10}(\frac{L}{d^2 A}) = M_0 -2.5 \log_{10}(L) + 5 \log_{10}(d) + 2.5 \log_{10}(A) Implementation .............. The goal is to create a surface brightness map from particles for which we know the magnitude. The first step is to convert those magnitudes into an extensive quantity akin to the flux at :math:`10\,\rm{pc}`. As the flux is proportional to the luminosity, we can use the previous relation, moving the proportionality constant to an additive constant that we can mergre with the zero point. Moreover, we ignore the zero point. This is possible a long as we carry the exact inverse transformation once moving from flux to magnitudes. In ``instrument.ComputeFLuMap`` this is coded the following way:: # convert to flux (ignore the zero point) F = 10**(-M/2.5) Where ``M`` is the magnitude for a given particle. Once this extensive quantity is defined, we can project particles on a grid and sum their contribution and get a map akin to a flux at :math:`10\,\rm{pc}`:: # store in the mass field nb.mass = F.astype(np.float32) # create the map FluxMap = self.ccd.MapUsingGauss(nb.pos,nb.mass,nb.rsp) The surface brightness is then computed in ``SurfaceBrightnessMap`` (``lib.FluxToSurfaceBrightness``). The three following steps are performed following the above derivation: Transform the extensive quantity akin to a flux (the zero point is ignored) into an absolute magnitude .. math:: M = -2.5 \log_{10}(F) Compute the apparent magnitude, i.e., add the contribution to the distance :math:`d` in unit of :math:`10\,\rm{pc}` .. math:: m = M - 5 \log_{10}(d) Compute the surface brightness magnitude .. math:: SB = m + 2.5 \log_{10}(A) where :math:`A` is the pixel area in :math:`arcsec^2` The corresponding code is:: # compute the absolute magnitude in each pixel (as before we ignore the zero point) c = FluxMap==0 FluxMap = np.where(FluxMap==0,1e-40, FluxMap) MagMap = np.where(c,0,-2.5*np.log10(FluxMap)) # compute the apparent magnitude in each pixel # Mpc -> 10pc d = object_distance.to(10*u.pc).value magMap = MagMap + 5*np.log10(d) # compute the surface brightness in each pixel SBMap = np.where(MagMap==0, 100, magMap + 2.5*np.log10(ccd_pixel_area.to(u.arcsec**2).value) ) .. _mockimgs_scripts-label: Scripts ------- ``mockimgs_sb_compute_images`` .............................. .. program-output:: mockimgs_sb_compute_images -h ``mockimgs_sb_extract_3Dfeatures_from_2Dmask`` .............................................. .. program-output:: mockimgs_sb_extract_3Dfeatures_from_2Dmask -h ``mockimgs_sb_show_images`` ........................... .. program-output:: mockimgs_sb_show_images -h ``mockimgs_sb_display_fits`` ............................ .. program-output:: mockimgs_sb_display_fits -h