# encoding: utf-8
#
# @Author: <NAME>, <NAME>
# @Date: Nov 15, 2021
# @Filename: ism.py
# @License: BSD 3-Clause
# @Copyright: <NAME>, <NAME>

import os.path
from astropy import units as u
from astropy import constants as c
import numpy as np
from astropy.io import fits, ascii
from astropy.table import Table
from scipy.special import sph_harm
from astropy.wcs import WCS
from astropy.wcs.utils import proj_plane_pixel_scales
from astropy.coordinates import SkyCoord
from astropy.modeling.models import Sersic2D
from dataclasses import dataclass
import sys
if (sys.version_info[0]+sys.version_info[1]/10.) < 3.8:
    from backports.cached_property import cached_property
else:
    from functools import cached_property
from scipy.ndimage.interpolation import map_coordinates
from scipy.interpolate import interp1d, interp2d
import lvmdatasimulator
from lvmdatasimulator import log
import progressbar
from joblib import Parallel, delayed
from astropy.convolution import convolve_fft, kernels
from lvmdatasimulator.utils import calc_circular_mask, convolve_array, set_default_dict_values, \
    ism_extinction, check_overlap, assign_units
fluxunit = u.erg / (u.cm ** 2 * u.s * u.arcsec ** 2)
velunit = u.km / u.s


def brightness_inhomogeneities_sphere(harm_amplitudes, ll, phi_cur, theta_cur, rho, med, radius, thickness):
    """
    Auxiliary function producing the inhomogeneities on the brightness distribution for the Cloud of Bubble objects
    using the spherical harmonics.
    """
    brt = theta_cur * 0
    for m in np.arange(-ll, ll + 1):
        brt += (harm_amplitudes[m + ll * (ll + 1) - 1] * sph_harm(m, ll, phi_cur, theta_cur).real * med *
                (1 - np.sqrt(abs(rho.value ** 2 / radius.value ** 2 - (1 - thickness / 2) ** 2))))
    return brt


def sphere_brt_in_line(brt_3d, rad_3d, rad_model, flux_model):
    """
    Auxiliary function computing the brightness of the Cloud or Bubble at given radii and in given line
    according to the Cloudy models
    """
    p = interp1d(rad_model, flux_model, fill_value='extrapolate', assume_sorted=True)
    return p(rad_3d) * brt_3d


def interpolate_sphere_to_cartesian(spherical_array, x_grid=None, y_grid=None, z_grid=None,
                                    rad_grid=None, theta_grid=None, phi_grid=None, pxscale=1. * u.pc):
    """
    Auxiliary function to project the brightness or velocities from the spherical to cartesian coordinates
    """
    x, y, z = np.meshgrid(x_grid, y_grid, z_grid, indexing='ij')
    phi_c, theta_c, rad_c = xyz_to_sphere(x, y, z, pxscale=pxscale)
    ir = interp1d(rad_grid, np.arange(len(rad_grid)), bounds_error=False)
    ith = interp1d(theta_grid, np.arange(len(theta_grid)))
    iphi = interp1d(phi_grid, np.arange(len(phi_grid)))
    new_ir = ir(rad_c.ravel())
    new_ith = ith(theta_c.ravel())
    new_iphi = iphi(phi_c.ravel())
    cart_data = map_coordinates(spherical_array, np.vstack([new_ir, new_ith, new_iphi]),
                                order=1, mode='constant', cval=0)
    return cart_data.reshape([len(x_grid), len(y_grid), len(z_grid)]).T


def limit_angle(value, bottom_limit=0, top_limit=np.pi):
    """
    Auxiliary function to limit the angle values to the range of [0, pi]
    """
    value[value < bottom_limit] += (top_limit - bottom_limit)
    value[value > top_limit] -= (top_limit - bottom_limit)
    return value


def xyz_to_sphere(x, y, z, pxscale=1. * u.pc):
    """
    Auxiliary function to map the coordinates from cartesian to spherical system
    """
    phi_c = np.arctan2(y, x)
    rad_c = (np.sqrt(x ** 2 + y ** 2 + z ** 2))
    rad_c[rad_c == 0 * u.pc] = 1e-3 * pxscale
    theta_c = (np.arccos(z / rad_c))
    phi_c = limit_angle(phi_c, 0 * u.radian, 2 * np.pi * u.radian)
    theta_c = limit_angle(theta_c, 0 * u.radian, np.pi * u.radian)
    return phi_c, theta_c, rad_c


def find_model_id(file=lvmdatasimulator.CLOUDY_MODELS,
                  check_id=None, params=lvmdatasimulator.CLOUDY_SPEC_DEFAULTS['id']):
    """
    Checks the input parameters of the pre-computed Cloudy model and return corresponding index in the grid
    """
    with fits.open(file) as hdu:
        if check_id is None:
            if params is None:
                check_id = lvmdatasimulator.CLOUDY_SPEC_DEFAULTS['id']
                log.warning(f'Default Cloudy model will be used (id = {check_id})')
            else:
                summary_table = Table(hdu['Summary'].data)
                indexes = np.arange(len(summary_table)).astype(int)
                rec_table = np.ones(shape=len(summary_table), dtype=bool)

                def closest(rec, prop, val):
                    unique_col = np.unique(summary_table[prop][rec])
                    if isinstance(val, str):
                        res = unique_col[unique_col == val]
                        if len(res) == 0:
                            return ""
                        return res
                    else:
                        return unique_col[np.argsort(np.abs(unique_col - val))[0]]

                for p in params:
                    if p not in summary_table.colnames or params[p] is None or \
                            ((isinstance(params[p], float) or isinstance(params[p], int)) and ~np.isfinite(params[p])):
                        continue
                    rec_table = rec_table & (summary_table[p] == closest(indexes, p, params[p]))
                    indexes = np.flatnonzero(rec_table)
                    if len(indexes) == 0:
                        break
                if len(indexes) == 0 or len(indexes) == len(summary_table):
                    check_id = lvmdatasimulator.CLOUDY_SPEC_DEFAULTS['id']
                    log.warning('Input parameters do not correspond to any pre-computed Cloudy model.'
                                'Default Cloudy model will be used (id = {0})'.format(check_id))
                elif len(indexes) == 1:
                    check_id = summary_table['Model_ID'][indexes[0]]
                    for p in params:
                        if p not in summary_table.colnames or params[p] is None or \
                                ((isinstance(params[p], float) or
                                  isinstance(params[p], int)) and ~np.isfinite(params[p])):
                            continue
                        if params[p] != summary_table[p][indexes[0]]:
                            log.warning(f'Use the closest pre-computed Cloudy model with id = {check_id}')
                            break
                else:
                    check_id = summary_table['Model_ID'][indexes[0]]
                    log.warning(f'Select one of the closest pre-computed Cloudy model with id = {check_id}')
                #
                # for cur_ext in range(len(hdu)):
                #     if cur_ext == 0:
                #         continue
                #     found = False
                #     for p in params:
                #         if p == 'id':
                #             continue
                #         precision = 1
                #         if p == 'Z':
                #             precision = 2
                #         if np.round(params[p], precision) != np.round(hdu[cur_ext].header[p], precision):
                #             break
                #     else:
                #         found = True
                #     if found:
                #         return cur_ext, check_id
                # check_id = lvmdatasimulator.CLOUDY_SPEC_DEFAULTS['id']
                # log.warning('Input parameters do not correspond to any pre-computed Cloudy model.'
                #             'Default Cloudy model will be used (id = {0})'.format(check_id))

        extension_index = None
        while extension_index is None:
            extension_index = [cur_ext for cur_ext in range(len(hdu)) if (
                    check_id == hdu[cur_ext].header.get('MODEL_ID'))]
            if len(extension_index) == 0:
                if check_id == lvmdatasimulator.CLOUDY_SPEC_DEFAULTS['id']:
                    log.warning('Model_ID = {0} is not found in the Cloudy models grid. '
                                'Use the first one in the grid instead'.format(check_id))
                    extension_index = 1
                else:
                    log.warning('Model_ID = {0} is not found in the Cloudy models grid. '
                                'Use default ({1}) instead'.format(check_id,
                                                                   lvmdatasimulator.CLOUDY_SPEC_DEFAULTS['id']))
                    check_id = lvmdatasimulator.CLOUDY_SPEC_DEFAULTS['id']
                    extension_index = None
            else:
                extension_index = extension_index[0]
    return extension_index, check_id


@dataclass
class Nebula:
    """
    Base class defining properties of every nebula type.
    By itself it describes the rectangular nebula (e.g. DIG)
    Constructed nebula has 4 dimensions, where 4th derive its appearance in different lines
        (if spectrum_id is None, or if it is dark nebula => only one line)
    """
    xc: int = None  # Center of the region in the field of view, pix
    yc: int = None  # Center of the region in the field of view, pix
    x0: int = 0  # Coordinates of the bottom-left corner in the field of view, pix
    y0: int = 0  # Coordinates of the bottom-left corner in the field of view, pix
    pix_width: int = None  # full width of cartesian grid, pix (should be odd)
    pix_height: int = None  # full height of cartesian grid, pix (should be odd)
    width: u.pc = 0 * u.pc  # width of the nebula in pc (not used if pix_width is set up)
    height: u.pc = 0 * u.pc  # height of the nebula in pc (not used if pix_height is set up)
    pxscale: u.pc = 0.01 * u.pc  # pixel size in pc
    spectrum_id: int = None  # ID of a template Cloudy emission spectrum for this nebula
    n_brightest_lines: int = None  # limit the number of the lines to the first N brightest
    sys_velocity: velunit = 0 * velunit  # Systemic velocity
    turbulent_sigma: velunit = 10 * velunit  # Velocity dispersion due to turbulence; included in calculations of LSF
    max_brightness: fluxunit = 1e-15 * fluxunit
    max_extinction: u.mag = 0 * u.mag
    perturb_scale: int = 0 * u.pc  # Spatial scale of correlated perturbations
    perturb_amplitude: float = 0.1  # Maximal amplitude of perturbations
    _npix_los: int = 1  # full size along line of sight in pixels
    nchunks: int = -1  # number of chuncks to use for the convolution. If negative, select automatically
    vel_gradient: (velunit / u.pc) = 0  # velocity gradient along the nebula
    vel_pa: u.degree = 0  # Position angle of the kinematical axis (for the velocity gradient or rotation velocity)

    def __post_init__(self):
        self._assign_all_units()
        self._assign_position_params()
        self._ref_line_id = 0
        self.linerat_constant = True  # True if the ratio of line fluxes shouldn't change across the nebula

    def _assign_all_units(self):
        whole_list_properties = ['pxscale', 'sys_velocity', 'turbulent_sigma', 'max_brightness', 'max_extinction',
                                 'perturb_scale', 'radius', 'PA', 'length', 'width', 'vel_gradient', 'r_eff',
                                 'vel_rot', 'expansion_velocity', 'spectral_axis', 'vel_pa']
        whole_list_units = [u.pc, velunit, velunit, fluxunit, u.mag, u.pc, u.pc, u.degree, u.pc, u.pc,
                            (velunit / u.pc), u.kpc, velunit, velunit, velunit, u.degree]
        cur_list_properties = []
        cur_list_units = []
        for prp, unit in zip(whole_list_properties, whole_list_units):
            if hasattr(self, prp):
                cur_list_properties.append(prp)
                cur_list_units.append(unit)
        assign_units(self, cur_list_properties, cur_list_units)

    def _assign_position_params(self, conversion_type='rect'):
        if conversion_type == 'rect':
            for v in ['height', 'width']:
                if self.__getattribute__(f'pix_{v}') is None:
                    val = np.round((self.__getattribute__(v) / self.pxscale).value / 2.).astype(int) * 2 + 1
                else:
                    val = np.round(self.__getattribute__(f'pix_{v}') / 2.).astype(int) * 2 + 1
                setattr(self, f'pix_{v}', val)
        elif conversion_type == 'ellipse':
            self.pix_width = (np.round(np.abs(self.radius / self.pxscale * np.sin(self.PA)) +
                              np.abs(self.radius / self.pxscale *
                                     self.ax_ratio * np.cos(self.PA))).astype(int) * 2 + 1).value
            self.pix_height = (np.round(np.abs(self.radius / self.pxscale * np.cos(self.PA)) +
                               np.abs(self.radius / self.pxscale *
                                      self.ax_ratio * np.sin(self.PA))).astype(int) * 2 + 1).value
        elif conversion_type == 'galaxy':
            self.pix_width = (np.round(np.abs(self.r_max * np.sin(self.PA)) +
                                       np.abs(self.r_max * self.ax_ratio * np.cos(self.PA))).astype(int) * 2 + 1).value
            self.pix_height = (np.round(np.abs(self.r_max * np.cos(self.PA)) +
                                        np.abs(self.r_max * self.ax_ratio * np.sin(self.PA))).astype(int) * 2 + 1).value
        elif conversion_type == 'cylinder':
            self.pix_width = (np.ceil((self.length * np.abs(np.sin(self.PA)) +
                                       self.width * np.abs(np.cos(self.PA))) / self.pxscale / 2.
                                      ).astype(int) * 2 + 1).value
            self.pix_height = (np.ceil((self.length * np.abs(np.cos(self.PA)) +
                                       self.width * np.abs(np.sin(self.PA))) / self.pxscale / 2.
                                       ).astype(int) * 2 + 1).value

        if (self.xc is not None) and (self.yc is not None):
            self.x0 = self.xc - np.round((self.pix_width - 1) / 2).astype(int)
            self.y0 = self.yc - np.round((self.pix_height - 1) / 2).astype(int)
        elif (self.x0 is not None) and (self.y0 is not None):
            self.xc = self.x0 + np.round((self.pix_width - 1) / 2).astype(int)
            self.yc = self.y0 + np.round((self.pix_height - 1) / 2).astype(int)

    @cached_property
    def _cartesian_x_grid(self):
        return np.arange(self.pix_width) * self.pxscale

    @cached_property
    def _cartesian_y_grid(self):
        return np.arange(self.pix_height) * self.pxscale

    @cached_property
    def _cartesian_z_grid(self):
        return np.arange(self._npix_los) * self.pxscale

    @cached_property
    def _max_density(self):
        return self.max_extinction * (1.8e21 / (u.cm ** 2 * u.mag))

    @cached_property
    def _brightness_3d_cartesian(self):
        """
        Method to obtain the brightness (or density) distribution of the nebula in cartesian coordinates
        """
        brt = np.ones(shape=(self.pix_height, self.pix_width, self._npix_los), dtype=float) / self._npix_los
        if (self.perturb_scale > 0) and (self.perturb_amplitude > 0):
            pertscale = (self.perturb_scale / self.pxscale).value
            perturb = np.random.uniform(-1, 1, (self.pix_height, self.pix_width)
                                        ) * self.perturb_amplitude / self._npix_los
            xx, yy = np.meshgrid(np.arange(self.pix_width), np.arange(self.pix_height))
            f = np.exp(-2 * (xx ** 2 + yy ** 2) / pertscale)
            perturb = 4 / np.sqrt(np.pi) / pertscale * np.fft.ifft2(np.fft.fft2(perturb) * np.fft.fft2(f)).real
            brt += (perturb[:, :, None] - np.median(perturb))
        return brt

    @cached_property
    def _brightness_4d_cartesian(self):
        """
        Derive the brightness (or density) distribution of the nebula for each emission line in cartesian coordinates
        """
        if self.spectrum_id is None or self.linerat_constant:
            flux_ratios = np.array([1.])
        else:
            with fits.open(lvmdatasimulator.CLOUDY_MODELS) as hdu:
                flux_ratios = hdu[self.spectrum_id].data[1:, 1]
                index_ha = np.flatnonzero(hdu[self.spectrum_id].data[1:, 0] == 6562.81)
                if self.n_brightest_lines is not None and \
                        (self.n_brightest_lines > 0) and (self.n_brightest_lines < len(flux_ratios)):
                    indexes_sorted = np.argsort(flux_ratios)[::-1]
                    flux_ratios = flux_ratios[indexes_sorted[: self.n_brightest_lines]]
                    index_ha = np.flatnonzero(hdu[self.spectrum_id].data[1:, 0][indexes_sorted] == 6562.81)
                if len(index_ha) == 1:
                    self._ref_line_id = index_ha[0]
        return self._brightness_3d_cartesian[None, :, :, :] * flux_ratios[:, None, None, None]

    @cached_property
    def brightness_skyplane(self):
        """
        Project the 3D nebula onto sky plane (for emission or continuum sources)
        """
        if self.max_brightness > 0:
            norm_max = self.max_brightness
        else:
            norm_max = 1
        map2d = np.nansum(self._brightness_3d_cartesian, 2)
        return map2d / np.nanmax(map2d) * norm_max

    @cached_property
    def brightness_skyplane_lines(self):
        """
        Project the 3D emission nebula line onto sky plane (return images in each emission line)
        """
        if self.max_brightness > 0:
            map2d = np.nansum(self._brightness_4d_cartesian, 3)
            return map2d / np.nanmax(map2d[self._ref_line_id, :, :]) * self.max_brightness
        else:
            return None

    @cached_property
    def extinction_skyplane(self):
        """
        Project the 3D nebula onto sky plane (for dark clouds)
        """
        if self.max_extinction > 0:
            map2d = np.nansum(self._brightness_3d_cartesian, 2)
            return map2d / np.nanmax(map2d) * self._max_density / (1.8e21 / (u.cm ** 2 * u.mag))
        else:
            return None

    @cached_property
    def vel_field(self):
        return self._get_2d_velocity()
        # if vel_field is None:
        #     return np.atleast_1d(self.sys_velocity)
        # else:
        #     return vel_field + self.sys_velocity

    def _get_2d_velocity(self):
        if hasattr(self, 'vel_gradient') and (self.vel_gradient is not None) and (self.vel_gradient != 0):
            xx, yy = np.meshgrid(np.arange(self.pix_width), np.arange(self.pix_height))
            vel_field = (- (xx - (self.pix_width - 1) / 2) * np.sin(self.vel_pa) +
                         (yy - (self.pix_height - 1) / 2) * np.cos(self.vel_pa)) * self.pxscale * self.vel_gradient
            return vel_field
        else:
            return None

    # @cached_property
    # def line_profile(self):
    #     lprf = np.zeros(shape=len(self.los_velocity), dtype=float)
    #     lprf[np.floor(len(lprf) / 2.).astype(int)] = 1.
    #     return lprf


@dataclass
class Rectangle(Nebula):
    """
    Class defining a simple rectangular component.
    This is equal to Nebula, but no perturbations and turbulence by default
    """
    perturb_amplitude: float = 0.0  # Maximal amplitude of perturbations
    turbulent_sigma: velunit = 0 * velunit  # Velocity dispersion due to turbulence; included in calculations of LSF

    def __post_init__(self):
        self._assign_all_units()
        self._assign_position_params()
        self._ref_line_id = 0
        self.linerat_constant = True  # True if the ratio of line fluxes shouldn't change across the nebula


@dataclass
class Ellipse(Nebula):
    """
    Class defining a simple elliptical component.
    No perturbations and turbulence by default
    """
    perturb_amplitude: float = 0.0  # Maximal amplitude of perturbations
    turbulent_sigma: velunit = 0 * velunit  # Velocity dispersion due to turbulence; included in calculations of LSF
    radius: u.pc = 1.0 * u.pc  # Radius along the major axis of the ellipse (or radius of the circle)
    PA: u.degree = 90 * u.degree  # position angle of the major axis
    ax_ratio: float = 1.  # ratio of minor/major axes

    def __post_init__(self):
        self._assign_all_units()
        self._npix_los = 1
        self._assign_position_params(conversion_type='ellipse')
        self._ref_line_id = 0
        self.linerat_constant = True  # True if the ratio of line fluxes shouldn't change across the nebula

    @cached_property
    def _brightness_3d_cartesian(self):
        """
        Method to obtain the brightness (or density) distribution of the nebula in cartesian coordinates
        """
        xx, yy = np.meshgrid(np.arange(self.pix_width), np.arange(self.pix_height))
        brt = np.ones(shape=(self.pix_height, self.pix_width), dtype=np.float32)
        angle = (self.PA + 90 * u.degree).to(u.radian).value
        xct = (xx - (self.pix_width - 1) / 2) * np.cos(angle) + \
              (yy - (self.pix_height - 1) / 2) * np.sin(angle)
        yct = (xx - (self.pix_width - 1) / 2) * np.sin(angle) - \
              (yy - (self.pix_height - 1) / 2) * np.cos(angle)
        rmaj = (self.radius.to(u.pc) / self.pxscale.to(u.pc)).value
        rmin = (self.radius.to(u.pc) / self.pxscale.to(u.pc)).value * self.ax_ratio
        rec = (xct ** 2 / rmaj ** 2) + (yct ** 2 / rmin ** 2) >= 1
        brt[rec] = 0
        brt = brt.reshape((self.pix_height, self.pix_width, 1))
        return brt


@dataclass
class Circle(Ellipse):
    """
    Class defining a simple circular component.
    """
    def __post_init__(self):
        self._assign_all_units()
        self.ax_ratio = 1.
        self._npix_los = 1
        self._assign_position_params(conversion_type='ellipse')
        self._ref_line_id = 0
        self.linerat_constant = True  # True if the ratio of line fluxes shouldn't change across the nebula


@dataclass
class Filament(Nebula):
    """
    Class of an isotropic cylindrical shape filament.
    Defined by its position, length, PA, radius, maximal optical depth.
    If it is emission-type filament, then also maximal brightness is required.
    Velocity gradient also can be set up
    """
    PA: u.degree = 90 * u.degree  # position angle of the filament
    length: u.pc = 10 * u.pc  # full length of the filament
    width: u.pc = 0.1 * u.pc  # full width (diameter) of the filament

    def __post_init__(self):
        self._assign_all_units()
        self._assign_position_params(conversion_type='cylinder')
        self._npix_los = 1
        self._ref_line_id = 0
        self.linerat_constant = True  # True if the ratio of line fluxes shouldn't change across the nebula

    @cached_property
    def _brightness_3d_cartesian(self):
        """
        Method to obtain the brightness (or density) distribution of the nebula in cartesian coordinates
        """
        xx, yy = np.meshgrid(np.arange(self.pix_width), np.arange(self.pix_height))
        brt = np.zeros_like(xx, dtype=np.float32)

        xct = (xx - (self.pix_width - 1) / 2) * np.cos(self.PA + 90 * u.degree) + \
              (yy - (self.pix_height - 1) / 2) * np.sin(self.PA + 90 * u.degree)
        yct = (xx - (self.pix_width - 1) / 2) * np.sin(self.PA + 90 * u.degree) - \
              (yy - (self.pix_height - 1) / 2) * np.cos(self.PA + 90 * u.degree)

        rad = ((self.width / self.pxscale).value / 2.)
        len_px = ((self.length / self.pxscale).value / 2.)
        rec = (np.abs(yct) <= rad) & (np.abs(xct) <= len_px)
        brt[rec] = np.sqrt(1. - (yct[rec] / rad) ** 2)
        brt = brt.reshape((self.pix_height, self.pix_width, 1))
        return brt


@dataclass
class _ObsoleteFilament(Nebula):
    """
    Class of an isotropic cylindrical shape filament.
    Defined by its position, length, PA, radius, maximal optical depth
    if it is emission-type filament, then maximal brightness
    NB: this class is obsolete, but might be considered later in case of implementation of varying line ratios
    """
    PA: u.degree = 90 * u.degree  # position angle of the filament
    length: u.pc = 10 * u.pc  # full length of the filament
    width: u.pc = 0.1 * u.pc  # full width (diameter) of the filament
    vel_gradient: (velunit / u.pc) = 0  # velocity gradient along the filament (to be added)
    _theta_bins: int = 50
    _rad_bins: int = 0
    _h_bins: int = 2
    _npix_los: int = 101

    def __post_init__(self):
        self._assign_all_units()
        if self._rad_bins == 0:
            self._rad_bins = np.ceil(self.width.to(u.pc).value / self.pxscale.to(u.pc).value * 5).astype(int)
        if (self.xc is not None) and (self.yc is not None):
            self.x0 = self.xc - np.round((len(self._cartesian_y_grid) - 1) / 2).astype(int)
            self.y0 = self.yc - np.round((len(self._cartesian_z_grid) - 1) / 2).astype(int)
        elif (self.x0 is not None) and (self.y0 is not None):
            self.xc = self.x0 + np.round((len(self._cartesian_y_grid) - 1) / 2).astype(int)
            self.yc = self.y0 + np.round((len(self._cartesian_z_grid) - 1) / 2).astype(int)
        self._ref_line_id = 0
        self.linerat_constant = True  # True if the ratio of line fluxes shouldn't change across the nebula

    @cached_property
    def _theta_grid(self):
        return np.linspace(0, 2 * np.pi, self._theta_bins)

    @cached_property
    def _h_grid(self):
        return np.linspace(0, self.length, self._h_bins)

    @cached_property
    def _rad_grid(self):
        return np.linspace(0, self.width / 2, self._rad_bins)

    @cached_property
    def _cartesian_y_grid(self):
        npix = np.ceil(1.01 * (self.length * np.abs(np.sin(self.PA)) +
                               self.width * np.abs(np.cos(self.PA))) / self.pxscale).astype(int)
        npix_l = npix / 2 - np.ceil(self.length / 2 * np.sin(-self.PA) / self.pxscale).astype(int)
        return (np.linspace(0, npix, npix + 1) - npix_l) * self.pxscale

    @cached_property
    def _cartesian_z_grid(self):
        npix = np.ceil(1.01 * (self.length * np.abs(np.cos(self.PA)) +
                               self.width * np.abs(np.sin(self.PA))) / self.pxscale).astype(int)
        npix_l = npix / 2 - np.ceil(self.length / 2 * np.cos(-self.PA) / self.pxscale).astype(int)
        return (np.linspace(0, npix, npix + 1) - npix_l) * self.pxscale

    @cached_property
    def _cartesian_x_grid(self):
        return np.linspace(-1.01, 1.01, self._npix_los) * self.width / 2

    @cached_property
    def _brightness_3d_cylindrical(self):
        """
        Method to calculate brightness (or opacity) of the cloud at given theta, phi and radii

        theta: float -- azimuthal angle [0, 2 * np.pi]
        rad: float -- radius [0, self.width / 2]
        h: float -- height [0, self.length]
        Returns:
            3D cube of normalized brightness in theta-rad-h grid; total brightness = 1
        """
        rho, theta, h = np.meshgrid(self._rad_grid, self._theta_grid, self._h_grid, indexing='ij')
        brt = np.ones_like(theta)
        brt[rho > (self.width / 2)] = 0
        norm = np.sum(brt)
        if norm > 0:
            brt = brt / np.sum(brt)
        return brt

    @cached_property
    def _brightness_3d_cartesian(self):
        x, y, z = np.meshgrid(self._cartesian_x_grid, self._cartesian_y_grid,
                              self._cartesian_z_grid, indexing='ij')

        h_c = -y * np.sin(self.PA) + z * np.cos(self.PA)
        theta_c = np.arctan2(y * np.cos(self.PA) + z * np.sin(self.PA), x)
        rad_c = np.sqrt(x ** 2 + (y * np.cos(self.PA) + z * np.sin(self.PA)) ** 2)
        rad_c[rad_c == 0 * u.pc] = 1e-3 * self.pxscale

        theta_c = limit_angle(theta_c, 0 * u.radian, 2 * np.pi * u.radian)

        ir = interp1d(self._rad_grid, np.arange(self._rad_bins), bounds_error=False)
        ith = interp1d(self._theta_grid, np.arange(self._theta_bins))
        ih = interp1d(self._h_grid, np.arange(self._h_bins), bounds_error=False)
        new_ir = ir(rad_c.ravel())
        new_ith = ith(theta_c.ravel())
        new_ih = ih(h_c.ravel())

        cart_data = map_coordinates(self._brightness_3d_cylindrical,
                                    np.vstack([new_ir, new_ith, new_ih]),
                                    order=1, mode='constant', cval=0)

        return cart_data.reshape([len(self._cartesian_x_grid),
                                  len(self._cartesian_y_grid),
                                  len(self._cartesian_z_grid)]).T


@dataclass
class Galaxy(Nebula):
    """
    Class defining the galaxy object (set up it as Sersic2D profile assuming it has continuum and emission components)
    """
    PA: u.degree = 90 * u.degree  # position angle of the major axis
    ax_ratio: float = 0.7  # ratio of minor/major axes
    r_eff: u.kpc = 1 * u.kpc  # Effective radius in kpc
    rad_lim: float = 3.  # Maximum radius for calculations (in R_eff)
    n: float = 1.  # Sersic index
    vel_rot: velunit = 0 * velunit  # Rotational velocity (not implemented yet)

    def __post_init__(self):
        self._assign_all_units()
        self._npix_los = 1
        self.r_max = self.r_eff.to(u.pc).value / self.pxscale.to(u.pc).value * self.rad_lim
        self._assign_position_params(conversion_type='galaxy')
        self._ref_line_id = 0
        self.linerat_constant = True  # True if the ratio of line fluxes shouldn't change across the nebula

    @cached_property
    def _brightness_3d_cartesian(self):
        """
        Method to obtain the brightness (or density) distribution of the nebula in cartesian coordinates
        """
        xx, yy = np.meshgrid(np.arange(self.pix_width), np.arange(self.pix_height))
        angle = (self.PA + 90 * u.degree).to(u.radian).value
        mod = Sersic2D(amplitude=1, r_eff=(self.r_eff.to(u.pc) / self.pxscale.to(u.pc)).value,
                       n=self.n, x_0=(self.pix_width - 1) / 2, y_0=(self.pix_height - 1) / 2,
                       ellip=1 - self.ax_ratio, theta=angle)
        brt = mod(xx, yy)

        xct = (xx - (self.pix_width - 1) / 2) * np.cos(angle) + \
              (yy - (self.pix_height - 1) / 2) * np.sin(angle)
        yct = (xx - (self.pix_width - 1) / 2) * np.sin(angle) - \
              (yy - (self.pix_height - 1) / 2) * np.cos(angle)
        rmaj = self.rad_lim * (self.r_eff.to(u.pc) / self.pxscale.to(u.pc)).value
        rmin = self.rad_lim * (self.r_eff.to(u.pc) / self.pxscale.to(u.pc)).value * self.ax_ratio
        mask = np.ones_like(brt, dtype=np.float32)
        rec = (xct ** 2 / rmaj ** 2) + (yct ** 2 / rmin ** 2) >= 1
        mask[rec] = 0
        mask = convolve_fft(mask, kernels.Gaussian2DKernel(3.), fill_value=0, allow_huge=True)
        brt = brt * mask
        brt = brt.reshape(self.pix_height, self.pix_width, 1)
        return brt

    def _get_2d_velocity(self):
        if hasattr(self, 'vel_rot') and (self.vel_rot is not None) and (self.vel_rot != 0):
            xx, yy = np.meshgrid(np.arange(self.pix_width), np.arange(self.pix_height))
            angle = (self.PA + 90 * u.degree).to(u.radian).value
            xct = (xx - (self.pix_width - 1) / 2) * np.cos(angle) + \
                  (yy - (self.pix_height - 1) / 2) * np.sin(angle)
            yct = (xx - (self.pix_width - 1) / 2) * np.sin(angle) - \
                  (yy - (self.pix_height - 1) / 2) * np.cos(angle)
            rad = np.sqrt(xct ** 2 + yct ** 2)
            vel_field = np.zeros_like(xx, dtype=np.float32) * velunit
            rec = rad > 0
            vel_field[rec] = self.vel_rot * np.sqrt(1 - self.ax_ratio ** 2) * xct[rec] / rad[rec]
            return vel_field
        else:
            return None



@dataclass
class DIG(Nebula):
    """
    Class defining the DIG component. For now it is defined just by its brightness (constant)
    """
    max_brightness: fluxunit = 1e-17 * fluxunit
    vel_gradient: (velunit / u.pc) = 0


@dataclass
class Cloud(Nebula):
    """Class of an isotropic spherical gas cloud without any ionization source.
    Defined by its position, radius, density, maximal optical depth"""
    radius: u.pc = 1.0 * u.pc
    max_brightness: fluxunit = 0 * fluxunit
    max_extinction: u.mag = 2.0 * u.mag
    thickness: float = 1.0
    perturb_degree: int = 0  # Degree of perturbations (max. degree of spherical harmonics for cloud)
    linerat_constant: bool = False  # True if the ratio of line fluxes shouldn't change across the nebula
    _phi_bins: int = 90
    _theta_bins: int = 90
    _rad_bins: int = 0
    _npix_los: int = 100

    def __post_init__(self):
        self._assign_all_units()
        if self._rad_bins == 0:
            self._rad_bins = np.ceil(self.radius.to(u.pc).value / self.pxscale.to(u.pc).value * 3).astype(int)
        delta = np.round((len(self._cartesian_y_grid) - 1) / 2).astype(int)
        if (self.xc is not None) and (self.yc is not None):
            self.x0 = self.xc - delta
            self.y0 = self.yc - delta
        elif (self.x0 is not None) and (self.y0 is not None):
            self.xc = self.x0 + delta
            self.yc = self.y0 + delta
        self._ref_line_id = 0

    @cached_property
    def _theta_grid(self):
        return np.linspace(0, np.pi, self._theta_bins)

    @cached_property
    def _phi_grid(self):
        return np.linspace(0, 2 * np.pi, self._phi_bins)

    @cached_property
    def _rad_grid(self):
        return np.linspace(0, self.radius, self._rad_bins)

    @cached_property
    def _cartesian_z_grid(self):
        npix = np.ceil(1.02 * self.radius / self.pxscale).astype(int)
        return np.linspace(-npix, npix, 2 * npix + 1) * self.pxscale

    @cached_property
    def _cartesian_y_grid(self):
        return self._cartesian_z_grid.copy()

    @cached_property
    def _cartesian_x_grid(self):
        return np.linspace(-1.02, 1.02, self._npix_los) * self.radius

    @cached_property
    def _brightness_3d_spherical(self):
        """
        Method to calculate brightness (or opacity) of the cloud at given theta, phi and radii

        theta: float -- polar angle [0, np.pi]
        phi: float -- azimuthal angle [0, 2 * np.pi]
        rad: float -- radius [0, self.radius]
        Returns:
            3D cube of normalized brightness in theta-phi-rad grid; total brightness = 1
        """
        rho, theta, phi = np.meshgrid(self._rad_grid, self._theta_grid, self._phi_grid, indexing='ij')
        brt = np.ones_like(theta)
        brt[rho < (self.radius * (1 - self.thickness))] = 0
        brt[rho > self.radius] = 0
        med = np.median(brt[brt > 0])
        if self.perturb_degree > 0:
            phi_cur = limit_angle(phi + np.random.uniform(0, 2 * np.pi, 1), 0, 2 * np.pi)
            theta_cur = limit_angle(theta + np.random.uniform(0, np.pi, 1), 0, np.pi)
            harm_amplitudes = self.perturb_amplitude * np.random.randn(self.perturb_degree * (self.perturb_degree + 2))

            brt += np.nansum(Parallel(n_jobs=lvmdatasimulator.n_process)(delayed(brightness_inhomogeneities_sphere)
                                                                         (harm_amplitudes, ll, phi_cur, theta_cur,
                                                                         rho, med, self.radius, self.thickness)
                                                                         for ll in np.arange(1,
                                                                                             self.perturb_degree + 1)),
                             axis=0)
            brt[brt < 0] = 0
        if med > 0:
            brt = brt / np.nansum(brt)
        return brt

    @cached_property
    def _brightness_4d_spherical(self):
        """
        Method to calculate brightness of the cloud at given theta, phi and radii for each line

        theta: float -- polar angle [0, np.pi]
        phi: float -- azimuthal angle [0, 2 * np.pi]
        rad: float -- radius [0, self.radius]
        Returns:
            4D cube of brightness in line-theta-phi-rad grid; normalized to the total brightness in Halpha
        """
        s = self._brightness_3d_spherical.shape
        if self.spectrum_id is None or self.linerat_constant:
            return self._brightness_3d_spherical.reshape((1, s[0], s[1], s[2]))
        rho, _, _ = np.meshgrid(self._rad_grid, self._theta_grid, self._phi_grid, indexing='ij')
        with fits.open(lvmdatasimulator.CLOUDY_MODELS) as hdu:
            radius = hdu[self.spectrum_id].data[0, 2:] * (self.thickness * self.radius) + \
                     self.radius * (1 - self.thickness)
            fluxes = hdu[self.spectrum_id].data[1:, 2:]
            radius = np.insert(radius, 0, self.radius * (1 - self.thickness))
            fluxes = np.insert(fluxes, 0, fluxes[:, 0], axis=1)
            index_ha = np.flatnonzero(hdu[self.spectrum_id].data[1:, 0] == 6562.81)
            if self.n_brightest_lines is not None and \
                    (self.n_brightest_lines > 0) and (self.n_brightest_lines < len(fluxes)):
                indexes_sorted = np.argsort(hdu[self.spectrum_id].data[1:, 1])[::-1]
                fluxes = fluxes[indexes_sorted[:self.n_brightest_lines], :]
                index_ha = np.flatnonzero(hdu[self.spectrum_id].data[1:, 0][indexes_sorted] == 6562.81)
            if len(index_ha) == 1:
                self._ref_line_id = index_ha[0]

            brt = np.array(Parallel(n_jobs=lvmdatasimulator.n_process)(delayed(sphere_brt_in_line)
                                                                       (self._brightness_3d_spherical, rho,
                                                                        radius, flux)
                                                                       for flux in fluxes)).reshape((fluxes.shape[0],
                                                                                                     s[0], s[1], s[2]))
            return brt / np.nansum(brt[self._ref_line_id])

    @cached_property
    def _brightness_3d_cartesian(self):
        return interpolate_sphere_to_cartesian(self._brightness_3d_spherical, x_grid=self._cartesian_x_grid,
                                               y_grid=self._cartesian_y_grid, z_grid=self._cartesian_z_grid,
                                               rad_grid=self._rad_grid, theta_grid=self._theta_grid,
                                               phi_grid=self._phi_grid, pxscale=self.pxscale)

    @cached_property
    def _brightness_4d_cartesian(self):
        s = self._brightness_4d_spherical.shape
        return np.array(Parallel(n_jobs=lvmdatasimulator.n_process)(delayed(interpolate_sphere_to_cartesian)
                                                                    (cur_line_array,
                                                                     self._cartesian_x_grid, self._cartesian_y_grid,
                                                                     self._cartesian_z_grid, self._rad_grid,
                                                                     self._theta_grid, self._phi_grid, self.pxscale)
                                                                    for cur_line_array in self._brightness_4d_spherical)
                        ).reshape((s[0], len(self._cartesian_z_grid), len(self._cartesian_y_grid),
                                   len(self._cartesian_x_grid)))


@dataclass
class Bubble(Cloud):
    """Class of an isotropic thin expanding bubble."""
    spectral_axis: velunit = np.arange(-20, 20, 10) * velunit
    expansion_velocity: velunit = 20 * velunit
    max_brightness: fluxunit = 1e-15 * fluxunit
    max_extinction: u.mag = 0 * u.mag
    thickness: float = 0.2

    @cached_property
    def _velocity_3d_spherical(self) -> velunit:
        """
        Calculate line of sight velocity at given radius, phi, theta

        V ~ 1/brightness (given that v~1/n_e^2 and brightness~ne^2)
        """
        rho, theta, phi = np.meshgrid(self._rad_grid, self._theta_grid, self._phi_grid, indexing='ij')
        vel_cube = np.zeros_like(self._brightness_3d_spherical)
        rec = (rho <= self.radius) & (rho >= (self.radius * (1 - self.thickness)))
        vel_cube[rec] = \
            np.sin(theta[rec]) * \
            np.cos(phi[rec]) * \
            self.expansion_velocity / self._brightness_3d_spherical[rec] * \
            np.median(self._brightness_3d_spherical[self._brightness_3d_spherical > 0])
        return vel_cube

    @cached_property
    def _velocity_3d_cartesian(self) -> velunit:
        return interpolate_sphere_to_cartesian(self._velocity_3d_spherical, x_grid=self._cartesian_x_grid,
                                               y_grid=self._cartesian_y_grid, z_grid=self._cartesian_z_grid,
                                               rad_grid=self._rad_grid, theta_grid=self._theta_grid,
                                               phi_grid=self._phi_grid, pxscale=self.pxscale)

    def _turbulent_lsf(self, velocity):
        """Line spread function as a function of coorinates, including the velocity center shift"""
        # mu = self.velocity(theta, phi)
        mu = self._velocity_3d_cartesian[:, :, :, None] * velunit + self.sys_velocity
        sig = self.turbulent_sigma
        return 1. / (np.sqrt(2. * np.pi) * sig) * np.exp(-np.power((velocity - mu) / sig, 2.) / 2)

    def _d_spectrum_cartesian(self, velocity: velunit):
        """Returns local spectrum, per pc**3 of area"""
        return (self._brightness_3d_cartesian[:, :, :, None] * (
                fluxunit / u.pc ** 3) * self._turbulent_lsf(velocity)).to(fluxunit / velunit / u.pc ** 3)

    @cached_property
    def line_profile(self) -> (fluxunit / velunit):
        """
        Produces the distribution of the observed line profiles in each pixels of the sky plane
        """
        vel_axis = self.spectral_axis.to(velunit, equivalencies=u.spectral())
        _, _, _, vels = np.meshgrid(self._cartesian_z_grid,
                                    self._cartesian_y_grid,
                                    self._cartesian_x_grid,
                                    vel_axis, indexing='ij')
        spectrum = (
                np.sum(self._d_spectrum_cartesian(vels), axis=2
                       ).T * (self._cartesian_x_grid[1] - self._cartesian_x_grid[0]
                              ) * (self._cartesian_y_grid[1] - self._cartesian_y_grid[0]
                                   ) * (self._cartesian_z_grid[1] - self._cartesian_z_grid[0])
        )
        return spectrum / np.sum(spectrum, axis=0)


@dataclass
class CustomNebula(Nebula):
    """
        Class defining the custom nebulae with the user-defined distribution of the brighntess, continuum
        and line shapes in different lines.
    """
    pass


@dataclass
class ISM:
    """
    Class defining the ISM contribution to the field of view
    """
    wcs: WCS
    width: int = 401  # Width of field of view in pixels
    height: int = 401  # Width of field of view in pixels
    spec_resolution: u.Angstrom = 0.06 * u.Angstrom  # Spectral resolution of the simulation
    npix_line: int = 1  # Minimal number of pixels for a resolution element at wl = 10000A for construction of vel.grid
    distance: u.kpc = 50 * u.kpc  # Distance to the object for further conversion between arcsec and pc
    sys_velocity: velunit = 0 * velunit  # Systemic velocity to center the vel.grid on
    vel_amplitude: velunit = 150 * velunit  # Maximal deviation from the systemic velocity to setup vel.grid
    turbulent_sigma: velunit = 20. * velunit  # turbulence vel. disp. to be used for every nebula unless other specified
    R_V: float = 3.1  # R_V value defining the reddening curve (to be used unless other value is provided for a nebula)
    ext_law: str = 'F99'  # Reddening law (to be used unless other value is provided for a nebula)

    def __post_init__(self):
        assign_units(self, ['vel_amplitude', 'turbulent_sigma', 'sys_velocity', 'distance', 'spec_resolution'],
                     [velunit, velunit, velunit, u.kpc, u.Angstrom])
        self.content = fits.HDUList()
        self.content.append(fits.PrimaryHDU(header=self.wcs.to_header(), data=np.zeros(shape=(2, 2), dtype=int)))
        self.vel_grid = np.linspace(-self.vel_amplitude + self.sys_velocity,
                                    self.vel_amplitude + self.sys_velocity,
                                    np.ceil(self.vel_amplitude / self.vel_resolution).astype(int) * 2 + 1)
        self.pxscale = proj_plane_pixel_scales(self.wcs)[0] * 3600 * self.distance.to(u.pc) / 206265.
        self.content[0].header['width'] = (self.width, "Width of field of view, px")
        self.content[0].header['height'] = (self.height, "Height of field of view, px")
        self.content[0].header['PhysRes'] = (self.pxscale.value, "Physical resolution, pc/px")
        self.content[0].header['Dist'] = (self.distance.value, "Distance, kpc")
        self.content[0].header['Vsys'] = (self.distance.value, "Systemic Velocity, km/s")
        self.content[0].header['VelRes'] = (self.vel_resolution.value, "Velocity resolution, km/s/px")
        self.content[0].header['TurbSig'] = (self.turbulent_sigma.value, "Default turbulent velocity dispersion, km/s")
        self.content[0].header['Nobj'] = (0, "Number of generated nebulae")

    @cached_property
    def vel_resolution(self):
        return (self.spec_resolution / self.npix_line / (10000 * u.Angstrom) * c.c).to(velunit)

    def _get_continuum(self, my_comp, wl_grid):
        """
        Properly extracts continuum for current nebula taking into account its shape and surface brightness
        :param my_comp:
        :param wl_grid:
        :return: continuum
        """
        cont_type = self.content[my_comp + "_CONTINUUM"].header.get("CONTTYPE")
        continuum = self.content[my_comp + "_CONTINUUM"].data
        cont_norm = self.content[my_comp + "_CONTINUUM"].header.get("CONTFLUX")
        cont_norm_wl = self.content[my_comp + "_CONTINUUM"].header.get("CONTWL")
        if cont_type.lower() == 'model':
            cont_wl_fullrange = continuum[0, :]
            cont_fullrange = continuum[1, :]
            p = interp1d(continuum[0, :], continuum[1, :], assume_sorted=True,
                         bounds_error=False, fill_value='extrapolate')
            continuum = p(wl_grid)
        elif cont_type.lower() == 'poly':
            p = np.poly1d(continuum)
            cont_wl_fullrange = np.linspace(3500, 10000, 501)
            cont_fullrange = p(cont_wl_fullrange)
            continuum = p(wl_grid)
        elif cont_type.lower() == 'bb':
            cont_wl_fullrange = np.linspace(3500, 10000, 501)
            cont_fullrange = 1 / cont_wl_fullrange ** 5 / (
                    np.exp(6.63e-27 * 3e10 / cont_wl_fullrange / 1e-8 / continuum / 1.38e-16) - 1)
            continuum = 1 / wl_grid ** 5 / (np.exp(6.63e-27 * 3e10 / wl_grid / 1e-8 / continuum / 1.38e-16) - 1)

        t_filter = None
        if type(cont_norm_wl) == str:
            file_filter = os.path.join(lvmdatasimulator.ROOT_DIR, 'data', 'instrument', 'filters', cont_norm_wl+".dat")
            if not os.path.isfile(file_filter):
                log.warning("Cannot find filter {0}. "
                            "Default Wavelength = 5500A will be used for continuum normalization".format(cont_norm_wl))
                cont_norm_wl = 5500.
                t_filter = None
            else:
                t_filter = ascii.read(file_filter, names=['lambda', 'transmission'])
                cont_norm_wl = np.sum(t_filter['lambda'] * t_filter['transmission']) / np.sum(t_filter['transmission'])
        if t_filter is None:
            cont_model_max = cont_fullrange[np.argmin(abs(cont_wl_fullrange - cont_norm_wl))]
        else:
            dl = np.roll(t_filter['lambda'], -1) - t_filter['lambda']
            dl[-1] = dl[-2]
            w_filter = np.sum(dl * t_filter['transmission'])/np.max(t_filter['transmission'])
            p = interp1d(t_filter['lambda'], t_filter['transmission'], assume_sorted=True,
                         fill_value=0, bounds_error=False)
            cont_model_max = np.sum(cont_fullrange * p(cont_wl_fullrange)) / w_filter
        if ~np.isfinite(cont_norm) or cont_norm <= 0:
            cont_norm = self.content[my_comp + "_CONTINUUM"].header.get("CONTMAG") * u.ABmag
            cont_norm = cont_norm.to(u.STmag, u.spectral_density(cont_norm_wl * u.AA)).to(u.erg/u.s/u.cm**2/u.AA).value
        return continuum / cont_model_max * cont_norm * (wl_grid[1] - wl_grid[0])

    def _add_fits_extension(self, name, value, obj_to_add, zorder=0, cur_wavelength=0, add_fits_kw=None,
                            add_counter=False):
        self.content.append(fits.ImageHDU(np.atleast_1d(value), name=name))
        self.content[-1].header['Nebtype'] = (type(obj_to_add).__name__, "Type of the nebula")
        is_dark = ((obj_to_add.max_brightness <= 0) and (obj_to_add.max_extinction > 0))
        self.content[-1].header['Dark'] = (is_dark, " Emitting or absorbing nebula?")
        self.content[-1].header['X0'] = (obj_to_add.x0, "Position in the field of view")
        self.content[-1].header['Y0'] = (obj_to_add.y0, "Position in the field of view")
        self.content[-1].header['Zorder'] = (zorder, "Z-order in the field of view")
        self.content[-1].header['NCHUNKS'] = (obj_to_add.nchunks, "N of chunks for convolution")
        if type(obj_to_add) in [Bubble, Cloud]:
            self.content[-1].header['Radius'] = (obj_to_add.radius.to_value(u.pc), "Radius of the nebula, pc")
            self.content[-1].header['PertOrd'] = (obj_to_add.perturb_degree, "Degree to produce random perturbations")
            self.content[-1].header['PertAmp'] = (obj_to_add.perturb_amplitude, "Max amplitude of random perturb.")
        if type(obj_to_add) in [DIG, Nebula, Rectangle, Ellipse, Circle]:
            self.content[-1].header['PertScl'] = (obj_to_add.perturb_scale.to_value(u.pc),
                                                  "Scale of the random perturbations, pc")
            self.content[-1].header['PertAmp'] = (obj_to_add.perturb_amplitude, "Max amplitude of random perturb.")
        if type(obj_to_add) in [Filament]:
            self.content[-1].header['Width'] = (obj_to_add.width.to_value(u.pc), 'Width of the filament, pc')
            self.content[-1].header['Length'] = (obj_to_add.length.to_value(u.pc), 'Length of the filament, pc')

        if type(obj_to_add) not in [Bubble, Cloud, Galaxy]:
            if obj_to_add.vel_gradient is not None:
                self.content[-1].header['Vgrad'] = (obj_to_add.vel_gradient.to_value(velunit / u.pc),
                                                    'Velocity gradient, km/s per pc')
                if obj_to_add.vel_pa is not None:
                    self.content[-1].header['PAkin'] = (obj_to_add.vel_pa.to_value(u.degree),
                                                        'Kinematical PA, degree')
        if type(obj_to_add) in [Nebula, Rectangle]:
            self.content[-1].header['Width'] = (obj_to_add.width.to_value(u.pc), 'Width of the nebula, pc')
            self.content[-1].header['Height'] = (obj_to_add.height.to_value(u.pc), 'Height of the nebula, pc')
        if type(obj_to_add) in [Ellipse, Filament, Galaxy]:
            self.content[-1].header['PA'] = (obj_to_add.PA.to_value(u.degree), 'Position angle, degree')
        if type(obj_to_add) in [Galaxy]:
            self.content[-1].header['Reff'] = (obj_to_add.r_eff.to_value(u.kpc), 'Effective radius of the galaxy, kpc')
            self.content[-1].header['NSersic'] = (obj_to_add.n, 'Sersic index')
            self.content[-1].header['Rlim'] = (obj_to_add.rad_lim, 'Limiting distance in Reff')
            self.content[-1].header['Vrot'] = (obj_to_add.vel_rot.to_value(velunit),
                                               'Rotational velocity, km/s')
            self.content[-1].header['PAkin'] = (obj_to_add.vel_pa.to_value(u.degree),
                                                'Kinematical PA, degree')
        if type(obj_to_add) in [Ellipse, Circle]:
            self.content[-1].header['Radius'] = (obj_to_add.radius.to_value(u.kpc), 'Radius (major axis), pc')
        if type(obj_to_add) in [Ellipse, Galaxy]:
            self.content[-1].header['AxRat'] = (obj_to_add.ax_ratio, "Axis ratio")
        if (obj_to_add.max_brightness <= 0) and (obj_to_add.max_extinction > 0):
            self.content[-1].header['MaxExt'] = (obj_to_add.max_extinction.value, "Max extinction, mag/pix")
        elif obj_to_add.max_brightness > 0:
            self.content[-1].header['MaxBrt'] = (obj_to_add.max_brightness.value, "Max brightness, erg/s/cm^2/arcsec^2")
        if type(obj_to_add) == Bubble:
            self.content[-1].header['Vexp'] = (obj_to_add.expansion_velocity.to_value(velunit),
                                               'Expansion velocity, km/s')
        self.content[-1].header['SysVel'] = (obj_to_add.sys_velocity.to_value(velunit), "Systemic velocity, km/s")
        self.content[-1].header['TurbVel'] = (obj_to_add.turbulent_sigma.to_value(velunit),
                                              "ISM Velocity dispersion, km/s")
        self.content[-1].header['SpecID'] = (obj_to_add.spectrum_id, "Ref. spectrum ID in model grid")
        self.content[-1].header['NLines'] = (obj_to_add.n_brightest_lines, "Maximal number of lines to use")
        if cur_wavelength:
            self.content[-1].header['Lambda'] = (cur_wavelength, "Current line wavelength")
        if add_fits_kw is not None:
            for kw in add_fits_kw:
                self.content[-1].header[kw] = add_fits_kw[kw]
        if add_counter:
            self.content[0].header['Nobj'] = (self.content[0].header['Nobj'] + 1, "Total number of nebulae")

    def add_nebula(self, obj_to_add, obj_id=0, zorder=0, add_fits_kw=None, continuum=None):
        """
        Method to add the particular nebula to the ISM object and to the output multi-extensions fits file
        """
        if type(obj_to_add) not in [Nebula, Bubble, Filament, DIG, Cloud, Galaxy, Ellipse, Circle,
                                    Rectangle, CustomNebula]:
            log.warning('Skip nebula of wrong type ({0})'.format(type(obj_to_add)))
            return

        if (obj_to_add.max_brightness <= 0) and (obj_to_add.max_extinction <= 0) and (continuum is None):
            log.warning('Skip nebula with zero extinction and brightness')
            return
        if obj_to_add.max_brightness > 0:
            brt = obj_to_add.brightness_skyplane.value
            if obj_to_add.spectrum_id is not None and not obj_to_add.linerat_constant:
                brt_4d = obj_to_add.brightness_skyplane_lines.value
            else:
                brt_4d = None
        elif obj_to_add.max_extinction > 0:
            brt = obj_to_add.extinction_skyplane.value
            brt_4d = None
        elif continuum is not None:
            brt = obj_to_add.brightness_skyplane
            brt_4d = None

        self._add_fits_extension(name="Comp_{0}_Brightness".format(obj_id), value=brt,
                                 obj_to_add=obj_to_add, zorder=zorder, add_fits_kw=add_fits_kw, add_counter=True)
        if obj_to_add.max_brightness > 0:
            if brt_4d is not None:
                with fits.open(lvmdatasimulator.CLOUDY_MODELS) as hdu:
                    wl_list = hdu[obj_to_add.spectrum_id].data[1:, 0]
                    if obj_to_add.n_brightest_lines is not None and \
                            (obj_to_add.n_brightest_lines > 0) and (obj_to_add.n_brightest_lines < len(wl_list)):
                        wl_list = wl_list[np.argsort(hdu[obj_to_add.spectrum_id].data[1:, 1]
                                                     )[::-1][: obj_to_add.n_brightest_lines]]
                for line_ind in range(brt_4d.shape[0]):
                    self._add_fits_extension(name="Comp_{0}_Flux_{1}".format(obj_id, wl_list[line_ind]),
                                             value=brt_4d[line_ind],
                                             obj_to_add=obj_to_add, zorder=zorder, add_fits_kw=add_fits_kw,
                                             cur_wavelength=wl_list[line_ind])
            elif obj_to_add.spectrum_id is not None and obj_to_add.linerat_constant:
                with fits.open(lvmdatasimulator.CLOUDY_MODELS) as hdu:
                    data_save = hdu[obj_to_add.spectrum_id].data[1:, :2]
                    if obj_to_add.n_brightest_lines is not None and \
                            (obj_to_add.n_brightest_lines > 0) and \
                            (obj_to_add.n_brightest_lines < len(hdu[obj_to_add.spectrum_id].data[1:, 0])):
                            data_save = data_save[np.argsort(data_save[:, 1])[::-1][: obj_to_add.n_brightest_lines]]
                    self._add_fits_extension(name="Comp_{0}_FluxRatios".format(obj_id),
                                             value=data_save.T,
                                             obj_to_add=obj_to_add, zorder=zorder, add_fits_kw=add_fits_kw)

        if type(obj_to_add) == Bubble:
            self._add_fits_extension(name="Comp_{0}_LineProfile".format(obj_id), value=obj_to_add.line_profile.value,
                                     obj_to_add=obj_to_add, zorder=zorder, add_fits_kw=add_fits_kw)

        if obj_to_add.vel_field is not None:
            self._add_fits_extension(name="Comp_{0}_Vel".format(obj_id), value=obj_to_add.vel_field.value,
                                     obj_to_add=obj_to_add, zorder=zorder, add_fits_kw=add_fits_kw)

        if continuum is not None:
            self._add_fits_extension(name="Comp_{0}_Continuum".format(obj_id), value=continuum,
                                     obj_to_add=obj_to_add, zorder=zorder, add_fits_kw=add_fits_kw)
        return self.content

    def save_ism(self, filename):
        self.content.writeto(filename, overwrite=True)
        log.info("Generated ISM saved to {0}".format(filename))

    def generate(self, all_objects):
        """
        Generate all the Nebulae from the input list

        Args:
            all_objects: list -- contains a dictionary describing the nebula to add:
                example:
                all_objects = [{type: "Bubble, Filament, DIG, ....",
                                sys_velocity: 0 * u.km/u.s,
                                expansion_velocity: 30 * u.km/u.s,
                                turbulent_sigma: 10 * u.km/u.s,
                                radius: 5 * u.pc,
                                max_brightness: 1e-16 * u.erg / u.cm**2 / u.s / u.arcsec ** 2,
                                RA: "08h12m13s",
                                DEC: "-20d14m13s",
                                'perturb_degree': 8, # max. order of spherical harmonics to generate inhomogeneities
                                'perturb_amplitude': 0.1, # relative max. amplitude of inhomogeneities,
                                'perturb_scale': 200 * u.pc,  # spatial scale to generate inhomogeneities (DIG only)
                                'distance': 50 * u.kpc,  # distance to the nebula (default is from ISM)
                                'cloudy_id': None,  # id of pre-generated Cloudy model
                                'cloudy_params': {'Z': 0.5, 'qH': 49., 'nH': 10, 'Teff': 30000, 'Geometry': 'Sphere'}, #
                                    parameters defining the pre-generated Cloudy model (used if spectrum_id is None)
                                'n_brightest_lines': 10,  # only this number of the brightest lines will be evaluated
                                'linerat_constant': False #  if True -> lines ratios don't vary across Cloud/Bubble
                                'continuum_type': 'BB' or 'Model' or 'Poly' # type of the continuum model
                                'continuum_data': model_id or [poly_coefficients] or Teff or dict with "wl" and "flux"
                                                # value defining cont. shape
                                'continuum_flux': 1e-16 * u.erg / u.cm ** 2 / u.s / u.arcsec **2 / u.AA,
                                'continuum_mag': 22 * u.mag,
                                'continuum_wl': 5500, # could be also R, V, B,
                                'ext_law': 'F99',  # Extinction law, one of those used by pyneb (used for dark nebulae)
                                'ext_rv': 3.1,  # Value of R_V for extinction curve calculation (used for dark nebulae)
                                'vel_gradient: 12. * u.km / u.s / u.pc  # Line-of-sight velocity gradient
                                'vel_pa: 30. * u.degree  # PA of kinematical axis (for vel_gradient or vel_rot)
                                }]
        """
        if type(all_objects) is dict:
            all_objects = [all_objects]
        if type(all_objects) not in [list, tuple]:
            log.warning('Cannot generate nebulae as the input is not a list or tuple')
            return None
        all_objects = [cobj for cobj in all_objects if cobj.get('type') in ['Nebula', 'Bubble', 'Galaxy',
                                                                            'Filament', 'DIG', 'Cloud',
                                                                            'Rectangle', 'Circle', 'Ellipse',
                                                                            'CustomNebula']]
        n_objects = len(all_objects)
        log.info("Start generating {} nebulae".format(n_objects))
        bar = progressbar.ProgressBar(max_value=n_objects).start()
        obj_id = self.content[0].header['Nobj']
        obj_id_ini = self.content[0].header['Nobj']
        for ind_obj, cur_obj in enumerate(all_objects):
            bar.update(ind_obj)
            # Setup default parameters for missing keywords
            kin_pa_default = 0
            if 'PA' in cur_obj:
                kin_pa_default = cur_obj['PA']
            for k, v, unit in zip(['max_brightness', 'max_extinction', 'thickness',
                                   'expansion_velocity', 'sys_velocity',
                                   'turbulent_sigma', 'perturb_degree',
                                   'perturb_amplitude', 'perturb_scale', 'radius', 'distance',
                                   'continuum_type', 'continuum_data', 'continuum_flux', 'continuum_mag',
                                   'continuum_wl', 'ext_law', 'ext_rv', 'vel_gradient', 'vel_rot', 'vel_pa',
                                   'n_brightest_lines', 'offset_RA', 'offset_DEC', 'RA', 'DEC'],
                                  [0, 0, 1., 0, self.sys_velocity, self.turbulent_sigma, 0, 0.1, 0, 0, self.distance,
                                   None, None, 0, None, 5500., self.ext_law, self.R_V, 0, 0, kin_pa_default, None,
                                   None, None, None, None],
                                  [fluxunit, u.mag, None, velunit, velunit, velunit, None, None,
                                   u.pc, u.pc, u.kpc, None, None, fluxunit/u.AA, None, u.Angstrom,
                                   None, None, velunit / u.pc, velunit, u.degree, None, u.arcsec, u.arcsec,
                                   u.degree, u.degree]):
                set_default_dict_values(cur_obj, k, v, unit=unit)

            for k in ['max_brightness', 'max_extinction', 'radius', 'continuum_flux']:
                if cur_obj[k] < 0:
                    cur_obj[k] = 0

            if (cur_obj['max_brightness'] == 0) and (cur_obj['max_extinction'] == 0) and \
                    (((cur_obj['continuum_mag'] is None) and (cur_obj['continuum_flux'] == 0)) or
                     (cur_obj['continuum_data'] is None) or (cur_obj['continuum_type'] is None)):
                log.warning("Wrong set of parameters for the nebula #{0}: skip this one".format(ind_obj))
                continue

            if lvmdatasimulator.CLOUDY_MODELS is None or (cur_obj['max_brightness'] <= 0):
                cloudy_model_index = None
                cloudy_model_id = None
            else:
                if cur_obj.get('cloudy_id') is None:
                    if cur_obj.get('cloudy_params') is None or (type(cur_obj.get('cloudy_params')) is not dict):
                        log.warning("Neither of 'cloudy_id' or 'cloudy_params' is set for the nebula #{0}: "
                                    "use default 'cloudy_id={1}'".format(ind_obj,
                                                                         lvmdatasimulator.CLOUDY_SPEC_DEFAULTS['id']))
                        cur_obj['cloudy_id'] = lvmdatasimulator.CLOUDY_SPEC_DEFAULTS['id']
                    else:
                        for p in lvmdatasimulator.CLOUDY_SPEC_DEFAULTS:
                            if p == 'id':
                                continue
                            if cur_obj['cloudy_params'].get(p) is None:
                                cur_obj['cloudy_params'][p] = lvmdatasimulator.CLOUDY_SPEC_DEFAULTS[p]

                cloudy_model_index, cloudy_model_id = find_model_id(file=lvmdatasimulator.CLOUDY_MODELS,
                                                                    check_id=cur_obj.get('cloudy_id'),
                                                                    params=cur_obj.get('cloudy_params'))

            if cur_obj.get('linerat_constant') is None:
                if cur_obj['type'] in ['Bubble', 'Cloud']:
                    cur_obj['linerat_constant'] = False
                else:
                    cur_obj['linerat_constant'] = True

            if cur_obj['type'] == 'DIG':
                if cur_obj.get('max_brightness') is None or cur_obj.get('max_brightness') <= 0:
                    log.warning("Wrong set of parameters for the nebula #{0}: skip this one".format(ind_obj))
                    continue
                if not cur_obj.get('zorder'):
                    cur_obj['zorder'] = -1
                if cur_obj['perturb_scale'] < 0:
                    cur_obj['perturb_scale'] = 0
                generated_object = DIG(max_brightness=cur_obj.get('max_brightness'),
                                       turbulent_sigma=cur_obj['turbulent_sigma'],
                                       sys_velocity=cur_obj['sys_velocity'],
                                       vel_gradient=cur_obj['vel_gradient'],
                                       vel_pa=cur_obj['vel_pa'],
                                       spectrum_id=cloudy_model_index,
                                       n_brightest_lines=cur_obj['n_brightest_lines'],
                                       pxscale=self.pxscale * (cur_obj['distance'].to(u.pc) / self.distance.to(u.pc)),
                                       perturb_scale=cur_obj['perturb_scale'],
                                       perturb_amplitude=cur_obj['perturb_amplitude'],
                                       pix_width=self.width, pix_height=self.height,
                                       )
            else:
                # ==== Check input parameters and do necessary conversions
                if not cur_obj.get('zorder'):
                    cur_obj['zorder'] = 0
                if not ((cur_obj.get('RA') is not None and cur_obj.get('DEC') is not None) or
                        (cur_obj.get('X') is not None and cur_obj.get('Y') is not None) or
                        (cur_obj.get('offset_X') is not None and cur_obj.get('offset_Y') is not None) or
                        (cur_obj.get('offset_RA') is not None and cur_obj.get('offset_DEC') is not None)):
                    log.warning("Wrong set of parameters for the nebula #{0}: skip this one".format(ind_obj))
                    continue
                if cur_obj['type'] in ['Rectangle', 'Nebula'] and not (('width' in cur_obj) and ('height' in cur_obj)):
                    log.warning("Wrong set of parameters for the nebula #{0}: skip this one".format(ind_obj))
                    continue
                if (cur_obj['type'] in ["Bubble", "Cloud", "Ellipse", 'Circle']) and (cur_obj['radius'] == 0):
                    log.warning("Wrong set of parameters for the nebula #{0}: skip this one".format(ind_obj))
                    continue
                if cur_obj['type'] == 'Filament' and not (('length' in cur_obj) and ('PA' in cur_obj)):
                    log.warning("Wrong set of parameters for the nebula #{0}: skip this one".format(ind_obj))
                    continue
                if cur_obj['type'] == 'Galaxy' and not (('r_eff' in cur_obj) and ('PA' in cur_obj) and
                                                        ('ax_ratio' in cur_obj)):
                    log.warning("Wrong set of parameters for the nebula #{0}: skip this one".format(ind_obj))
                    continue
                if cur_obj['type'] == 'Galaxy':
                    if 'n' not in cur_obj:
                        log.info("Set default Sersic index n=1 for the nebula #{0}".format(ind_obj))
                        cur_obj['n'] = 1
                    if 'rad_lim' not in cur_obj:
                        cur_obj['rad_lim'] = 3.

                if not (cur_obj.get('X') is not None and cur_obj.get('Y') is not None):
                    if cur_obj.get('offset_X') is not None and cur_obj.get('offset_Y') is not None:
                        x = (self.width - 1) / 2. + cur_obj.get('offset_X')
                        y = (self.height - 1) / 2. + cur_obj.get('offset_Y')
                    elif cur_obj.get('RA') is not None and cur_obj.get('DEC') is not None:
                        radec = SkyCoord(ra=cur_obj.get('RA'), dec=cur_obj.get('DEC'))
                        x, y = self.wcs.world_to_pixel(radec)
                    elif cur_obj.get('offset_RA') is not None and cur_obj.get('offset_DEC') is not None:
                        x = (self.width - 1) / 2. - (cur_obj.get('offset_RA').to_value(u.degree) /
                                                     proj_plane_pixel_scales(self.wcs)[0])
                        y = (self.height - 1) / 2. + (cur_obj.get('offset_DEC').to_value(u.degree) /
                                                     proj_plane_pixel_scales(self.wcs)[0])
                    x = np.round(x).astype(int)
                    y = np.round(y).astype(int)
                else:
                    x, y = [cur_obj.get('X'), cur_obj.get('Y')]

                if (cur_obj['thickness'] <= 0) or (cur_obj['thickness'] > 1):
                    log.warning("Wrong value of thickness of the nebula #{0}: set it to 1.".format(ind_obj))
                    cur_obj['thickness'] = 1.

                if cur_obj['type'] == "Bubble" and cur_obj.get('expansion_velocity') <= 0:
                    log.warning("Contracting bubbles are not supported (nebula #{0})."
                                " Use non-expanding cloud instead".format(ind_obj))
                    cur_obj['type'] = "Cloud"

                if cur_obj['type'] in ["Bubble", "Cloud"]:
                    if cur_obj['perturb_degree'] < 0:
                        cur_obj['perturb_degree'] = 0

                elif cur_obj['type'] == 'Filament' and ('width' not in cur_obj):
                    log.info("Set default width of the filament 0.1 pc for the nebula #{0}".format(ind_obj))
                    cur_obj['width'] = 0.1 * u.pc

                # ==== Start calculations of different nebulae
                if cur_obj['type'] == "Bubble":
                    generated_object = Bubble(xc=x, yc=y,
                                              max_brightness=cur_obj.get('max_brightness'),
                                              max_extinction=cur_obj.get('max_extinction'),
                                              spectral_axis=self.vel_grid,
                                              expansion_velocity=cur_obj.get('expansion_velocity'),
                                              thickness=cur_obj['thickness'],
                                              radius=cur_obj['radius'],
                                              pxscale=self.pxscale * (cur_obj['distance'].to(u.pc) /
                                                                      self.distance.to(u.pc)),
                                              perturb_degree=cur_obj['perturb_degree'],
                                              perturb_amplitude=cur_obj['perturb_amplitude'],
                                              turbulent_sigma=cur_obj['turbulent_sigma'],
                                              sys_velocity=cur_obj['sys_velocity'],
                                              spectrum_id=cloudy_model_index,
                                              n_brightest_lines=cur_obj['n_brightest_lines'],
                                              linerat_constant=cur_obj['linerat_constant'],
                                              )
                elif cur_obj['type'] == "Cloud":
                    generated_object = Cloud(xc=x, yc=y,
                                             max_brightness=cur_obj.get('max_brightness'),
                                             max_extinction=cur_obj.get('max_extinction'),
                                             thickness=cur_obj['thickness'],
                                             radius=cur_obj['radius'],
                                             pxscale=self.pxscale * (cur_obj['distance'].to(u.pc) /
                                                                     self.distance.to(u.pc)),
                                             perturb_degree=cur_obj['perturb_degree'],
                                             perturb_amplitude=cur_obj['perturb_amplitude'],
                                             spectrum_id=cloudy_model_index,
                                             n_brightest_lines=cur_obj['n_brightest_lines'],
                                             turbulent_sigma=cur_obj['turbulent_sigma'],
                                             sys_velocity=cur_obj['sys_velocity'],
                                             linerat_constant=cur_obj['linerat_constant'],
                                             vel_gradient=cur_obj['vel_gradient'],
                                             vel_pa=cur_obj['vel_pa'],
                                             )

                elif cur_obj['type'] == "Filament":
                    generated_object = Filament(xc=x, yc=y,
                                                max_brightness=cur_obj.get('max_brightness'),
                                                max_extinction=cur_obj.get('max_extinction'),
                                                width=cur_obj['width'],
                                                length=cur_obj['length'],
                                                PA=cur_obj['PA'],
                                                vel_gradient=cur_obj['vel_gradient'],
                                                vel_pa=cur_obj['vel_pa'],
                                                spectrum_id=cloudy_model_index,
                                                n_brightest_lines=cur_obj['n_brightest_lines'],
                                                turbulent_sigma=cur_obj['turbulent_sigma'],
                                                sys_velocity=cur_obj['sys_velocity'],
                                                pxscale=self.pxscale * (cur_obj['distance'].to(u.pc) /
                                                                        self.distance.to(u.pc)),
                                                )

                elif cur_obj['type'] == "Galaxy":
                    generated_object = Galaxy(xc=x, yc=y,
                                              max_brightness=cur_obj.get('max_brightness'),
                                              max_extinction=cur_obj.get('max_extinction'),
                                              r_eff=cur_obj['r_eff'],
                                              rad_lim=cur_obj['rad_lim'],
                                              ax_ratio=cur_obj['ax_ratio'],
                                              PA=cur_obj['PA'],
                                              n=cur_obj['n'],
                                              vel_rot=cur_obj.get('vel_rot'),
                                              vel_pa=cur_obj['vel_pa'],
                                              spectrum_id=cloudy_model_index,
                                              n_brightest_lines=cur_obj['n_brightest_lines'],
                                              turbulent_sigma=cur_obj['turbulent_sigma'],
                                              sys_velocity=cur_obj['sys_velocity'],
                                              pxscale=self.pxscale * (cur_obj['distance'].to(u.pc) /
                                                                      self.distance.to(u.pc)),
                                              )
                elif cur_obj['type'] == "Ellipse":
                    generated_object = Ellipse(xc=x, yc=y,
                                               max_brightness=cur_obj.get('max_brightness'),
                                               max_extinction=cur_obj.get('max_extinction'),
                                               radius=cur_obj['radius'],
                                               ax_ratio=cur_obj['ax_ratio'],
                                               PA=cur_obj['PA'],
                                               spectrum_id=cloudy_model_index,
                                               n_brightest_lines=cur_obj['n_brightest_lines'],
                                               turbulent_sigma=cur_obj['turbulent_sigma'],
                                               sys_velocity=cur_obj['sys_velocity'],
                                               vel_gradient=cur_obj['vel_gradient'],
                                               vel_pa=cur_obj['vel_pa'],
                                               pxscale=self.pxscale * (cur_obj['distance'].to(u.pc) /
                                                                       self.distance.to(u.pc)),
                                               perturb_scale=cur_obj['perturb_scale'],
                                               perturb_amplitude=cur_obj['perturb_amplitude'],
                                               )
                elif cur_obj['type'] == "Circle":
                    generated_object = Circle(xc=x, yc=y,
                                              max_brightness=cur_obj.get('max_brightness'),
                                              max_extinction=cur_obj.get('max_extinction'),
                                              radius=cur_obj['radius'],
                                              spectrum_id=cloudy_model_index,
                                              n_brightest_lines=cur_obj['n_brightest_lines'],
                                              turbulent_sigma=cur_obj['turbulent_sigma'],
                                              sys_velocity=cur_obj['sys_velocity'],
                                              vel_gradient=cur_obj['vel_gradient'],
                                              vel_pa=cur_obj['vel_pa'],
                                              pxscale=self.pxscale * (cur_obj['distance'].to(u.pc) /
                                                                      self.distance.to(u.pc)),
                                              perturb_scale=cur_obj['perturb_scale'],
                                              perturb_amplitude=cur_obj['perturb_amplitude'],
                                              )
                elif cur_obj['type'] == "Rectangle" or (cur_obj['type'] == "Nebula"):
                    generated_object = Rectangle(xc=x, yc=y,
                                                 width=cur_obj.get('width'), height=cur_obj.get('height'),
                                                 max_brightness=cur_obj.get('max_brightness'),
                                                 max_extinction=cur_obj.get('max_extinction'),
                                                 spectrum_id=cloudy_model_index,
                                                 n_brightest_lines=cur_obj['n_brightest_lines'],
                                                 turbulent_sigma=cur_obj['turbulent_sigma'],
                                                 sys_velocity=cur_obj['sys_velocity'],
                                                 vel_gradient=cur_obj['vel_gradient'],
                                                 vel_pa=cur_obj['vel_pa'],
                                                 pxscale=self.pxscale * (cur_obj['distance'].to(u.pc) /
                                                                         self.distance.to(u.pc)),
                                                 perturb_scale=cur_obj['perturb_scale'],
                                                 perturb_amplitude=cur_obj['perturb_amplitude'],
                                                 )
                elif cur_obj['type'] == "CustomNebulae":
                    # generated_object = CustomNebula(xc=x, yc=y,)
                    log.warning("Custom Nebulae will be added soon")
                    continue
                else:
                    log.warning("Unrecognized type of the nebula #{0}: skip this one".format(ind_obj))
                    continue
            if cloudy_model_index is not None:
                if cur_obj['linerat_constant']:
                    lr = "Constant"
                else:
                    lr = "Variable"
                add_fits_kw = {"Model_ID": cloudy_model_id, "LineRat": lr}
            else:
                add_fits_kw = {}

            add_fits_kw["Distance"] = (cur_obj['distance'].to_value(u.kpc), "Distance to the nebula, kpc")

            continuum = None
            if cur_obj['continuum_type'] is not None and cur_obj['continuum_data'] is not None \
                    and cur_obj['continuum_type'].lower() in ['bb', 'poly', 'model']:

                if cur_obj['continuum_type'].lower() == 'model':
                    if isinstance(cur_obj['continuum_data'], dict) and ('wl' in cur_obj['continuum_data']) and \
                            ('flux' in cur_obj['continuum_data']):
                        if len(cur_obj['continuum_data']['wl']) != len(cur_obj['continuum_data']['flux']):
                            log.error("Number of wavelength and flux points is inconsistent for continuum")
                        else:
                            wlscale = cur_obj['continuum_data']['wl']
                            continuum = np.vstack([wlscale, cur_obj['continuum_data']['flux']])
                    elif ~isinstance(cur_obj['continuum_data'], dict) and lvmdatasimulator.CONTINUUM_MODELS is not None:
                        with fits.open(lvmdatasimulator.CONTINUUM_MODELS) as hdu:
                            if cur_obj['continuum_data'] >= hdu[0].data.shape[0]:
                                log.warning("Wrong continuum model ID for nebula #{0}".format(obj_id))
                            else:
                                wlscale = (np.arange(hdu[0].data.shape[1]) - hdu[0].header['CRPIX1'] + 1
                                           ) * hdu[0].header['CDELT1'] + hdu[0].header['CRVAL1']
                                continuum = np.vstack([wlscale, hdu[0].data[cur_obj['continuum_data']]])
                elif cur_obj['continuum_type'].lower() in ['poly', 'bb']:
                    continuum = cur_obj['continuum_data']
            if continuum is not None:
                if add_fits_kw is None:
                    add_fits_kw = {}
                add_fits_kw['CONTTYPE'] = (cur_obj['continuum_type'], "Type of the continuum")
                if cur_obj['continuum_flux'] > 0:
                    contflux = cur_obj['continuum_flux'].to_value(u.erg / u.cm ** 2 / u.s / u.arcsec ** 2 / u.AA)
                else:
                    contflux = 0
                add_fits_kw['CONTFLUX'] = (contflux,
                                           "Continuum brightness (in erg/s/cm^2/asec^2/AA) at ref. wl/Filter")
                if cur_obj['continuum_mag'] is not None:
                    contmag = cur_obj['continuum_mag'].to_value(u.mag)
                else:
                    contmag = None
                add_fits_kw['CONTMAG'] = (contmag,
                                          "Continuum brightness (in mag/asec^2) at ref. wl/Filter")
                if isinstance(cur_obj['continuum_wl'], str):
                    cont_wl = cur_obj['continuum_wl']
                else:
                    cont_wl = cur_obj['continuum_wl'].to_value(u.AA)
                add_fits_kw['CONTWL'] = (cont_wl, 'Reference wavelength/filter for cont. flux/mag')

            if cur_obj.get('max_extinction') > 0:
                if add_fits_kw is None:
                    add_fits_kw = {}
                add_fits_kw['EXT_LAW'] = (cur_obj['ext_law'], "Extinction law according to pyneb list")
                add_fits_kw['EXT_RV'] = (cur_obj['ext_rv'], "R_V value for extinction calculations")

            self.add_nebula(generated_object, obj_id=obj_id, zorder=cur_obj.get('zorder'), add_fits_kw=add_fits_kw,
                            continuum=continuum)
            obj_id += 1
        bar.finish()
        if (obj_id - obj_id_ini) == 0:
            return None
        else:
            return True

    def load_nebulae(self, file):
        """
        Load previously saved fits-file containing the information about all nebulae.
        Note: Grid should be equal to that used when fits-file was generated!
        """
        if not os.path.isfile(file):
            log.warning("ISM doesn't contain any nebula")
            return None
        else:
            with fits.open(file) as hdu:
                wcs = WCS(hdu[0].header)
                cdelt_file = [cdelt.to(u.degree).value for cdelt in wcs.proj_plane_pixel_scales()]
                cdelt_ism = [cdelt.to(u.degree).value for cdelt in self.wcs.proj_plane_pixel_scales()]
                check = ~np.isclose(cdelt_file, cdelt_ism)
                check = np.append(check, ~np.isclose(wcs.wcs.crval, self.wcs.wcs.crval))
                if any(np.array([cur_hdu.header.get('EXTNAME') for cur_hdu in hdu
                                 if (cur_hdu.header.get('EXTNAME') is not None) and
                                    ("_LINEPROFILE" in cur_hdu.header.get('EXTNAME'))])):
                    check = np.append(check, ~np.isclose(hdu[0].header.get("VELRES"), self.vel_resolution.value))
                if any(check):
                    log.warning("Grid of fits file is inconsistent with that defined in ISM")
                    return None
                self.content = fits.HDUList()
                for cur_ind, hh in enumerate(hdu):
                    if cur_ind == 0:
                        self.content.append(fits.PrimaryHDU(header=hh.header))
                    else:
                        self.content.append(fits.ImageHDU(header=hh.header, data=hh.data,
                                                          name=hh.header.get('EXTNAME')))
                return True

    def calc_extinction(self, wavelength=6562.81, x0=0, y0=0, xs=None, ys=None, extension_name=None):
        """
        Calculate coefficient to reduce flux due to extinction at given wavelength(s)

        Args:
            x0: start x-coordinate in the field of view for calculations
            y0: start x-coordinate in the field of view for calculations
            xs: x-size (in pixels) of the area for calculations (if None => then just pixel x0,y0 is considered; xs=1)
            ys: y-size (in pixels) of the area for calculations (if None => then just pixel x0,y0 is considered; ys=1)
            wavelength: in angstrom, particular wavelength (or wavelengths)
                at which the calculations should be performed
            extension_name (str): name of the extension for current dark nebula

        Returns:
            None (if no any dark nebula at particular location) or np.array of (nlines, ys, xs) shape
        """
        if self.content[0].header['Nobj'] == 0 or (x0 > self.width) or (y0 > self.height):
            return None
        if xs is None:
            x1 = x0
        else:
            x1 = np.clip(x0 + xs - 1, 0, self.width - 1)
        if ys is None:
            y1 = y0
        else:
            y1 = np.clip(y0 + ys - 1, 0, self.height - 1)

        def check_in_region(reg_x0, reg_y0, reg_xs, reg_ys):
            if ((reg_x0 + reg_xs) < x0) or (reg_x0 > x1) or ((reg_y0 + reg_ys) < y0) or (reg_y0 > y1):
                return False
            else:
                return True

        if extension_name is None:
            all_dark_nebulae = [hdu.header.get('EXTNAME') for hdu in self.content if
                                hdu.header.get('EXTNAME') is not None and
                                ("BRIGHTNESS" in hdu.header.get('EXTNAME')) and
                                hdu.header.get('DARK') and
                                check_in_region(hdu.header.get('X0'), hdu.header.get('Y0'),
                                                hdu.header.get('NAXIS1'), hdu.header.get('NAXIS2'))]
        else:
            if not check_in_region(self.content[extension_name].header.get('X0'),
                                   self.content[extension_name].header.get('Y0'),
                                   self.content[extension_name].header.get('NAXIS1'),
                                   self.content[extension_name].header.get('NAXIS2')):
                return None
            all_dark_nebulae = [extension_name]
        if len(all_dark_nebulae) == 0:
            return None
        if type(wavelength) in [float, int, np.float64, np.float32]:
            wavelength = np.atleast_1d(wavelength)
        ext_map = np.ones(shape=(len(wavelength), y1 - y0 + 1, x1 - x0 + 1), dtype=np.float32)

        for dark_nebula in all_dark_nebulae:
            cur_neb_av = np.zeros(shape=(y1 - y0 + 1, x1 - x0 + 1), dtype=np.float32)
            cur_neb_x0 = np.clip(x0 - self.content[dark_nebula].header.get('X0'), 0, None)
            cur_neb_y0 = np.clip(y0 - self.content[dark_nebula].header.get('Y0'), 0, None)
            cur_neb_x1 = self.content[dark_nebula].header.get('NAXIS1') - 1 - np.clip(
                self.content[dark_nebula].header.get('X0') + self.content[dark_nebula].header.get('NAXIS1') - 1 - x1,
                0, None)
            cur_neb_y1 = self.content[dark_nebula].header.get('NAXIS2') - 1 - np.clip(
                self.content[dark_nebula].header.get('Y0') + self.content[dark_nebula].header.get('NAXIS2') - 1 - y1,
                0, None)
            cur_neb_av[cur_neb_y0 + self.content[dark_nebula].header.get('Y0') - y0:
                       cur_neb_y1 + self.content[dark_nebula].header.get('Y0') - y0 + 1,
                       cur_neb_x0 + self.content[dark_nebula].header.get('X0') - x0:
                       cur_neb_x1 + self.content[dark_nebula].header.get('X0') - x0 + 1
                       ] = self.content[dark_nebula].data[cur_neb_y0: cur_neb_y1 + 1, cur_neb_x0: cur_neb_x1 + 1]
            cur_extinction_law = self.ext_law
            cur_r_v = self.R_V
            if self.content[dark_nebula].header.get('EXT_LAW'):
                cur_extinction_law = self.content[dark_nebula].header.get('EXT_LAW')
            if self.content[dark_nebula].header.get('EXT_RV'):
                cur_r_v = self.content[dark_nebula].header.get('EXT_RV')
            ext_map = ext_map * ism_extinction(av=cur_neb_av, r_v=cur_r_v,
                                               ext_law=cur_extinction_law, wavelength=wavelength)
        return ext_map

    def get_map(self, wavelength=6562.81, get_continuum=False):
        """
        Method to produce 2D map of all ISM nebulae in the selected line
        Args:
            wavelength (float or iterative): wavelength (in Angstrom) according to the lines list, or wavelength range
            get_continuum (bool): if True, then also counts the flux from the continuum
        """
        wavelength = np.atleast_1d(wavelength)
        if len(wavelength) == 1:
            wavelength = np.array([wavelength[0]-0.01, wavelength[0]+0.01])
        if self.content[0].header['Nobj'] == 0:
            log.warning("ISM doesn't contain any nebula")
            return None
        all_extensions = [hdu.header.get('EXTNAME') for hdu in self.content]
        all_extensions_brt = np.array([extname for extname in all_extensions
                                       if extname is not None and ("BRIGHTNESS" in extname) and
                                       check_overlap(self.content[extname], (self.height, self.width))])

        if all([self.content[cur_ext].header.get("DARK") for cur_ext in all_extensions_brt]):
            # !!!! ADD later accounting of the continuum and extinction from those nebulae
            log.warning("ISM doesn't contain any emission nebula")
            return None

        all_extensions_brt = all_extensions_brt[
            np.argsort([self.content[cur_ext].header.get('ZORDER') for cur_ext in all_extensions_brt])]

        map_2d = np.zeros(shape=(self.height, self.width), dtype=np.float32)
        map_is_empty = True
        for cur_ext in all_extensions_brt:
            my_comp = "_".join(cur_ext.split("_")[:2])

            y0_in_field = np.clip(self.content[cur_ext].header['Y0'], 0, None)
            y1_in_field = np.clip(self.content[cur_ext].header['Y0'] + self.content[cur_ext].header['NAXIS2'] - 1, None,
                                  self.height-1)
            x0_in_field = np.clip(self.content[cur_ext].header['X0'], 0, None)
            x1_in_field = np.clip(self.content[cur_ext].header['X0'] + self.content[cur_ext].header['NAXIS1'] - 1, None,
                                  self.width - 1)

            if get_continuum and (my_comp + "_CONTINUUM" in all_extensions):
                n_wl_bins = int(np.clip(np.ceil((wavelength[1]-wavelength[0])/20.), 10, 200))
                wl_grid = np.linspace(wavelength[0], wavelength[1], n_wl_bins)
                continuum = np.sum(self._get_continuum(my_comp, wl_grid))
                add_continuum = self.content[cur_ext].data / np.max(self.content[cur_ext].data) * continuum
                map_2d[y0_in_field: y1_in_field + 1, x0_in_field: x1_in_field + 1] += \
                    add_continuum[y0_in_field - self.content[cur_ext].header['Y0']:
                                  y1_in_field - self.content[cur_ext].header['Y0'] + 1,
                                  x0_in_field - self.content[cur_ext].header['X0']: x1_in_field -
                                  self.content[cur_ext].header['X0'] + 1
                                  ]
                map_is_empty = False
            if self.content[cur_ext].header.get("DARK"):
                if map_is_empty:
                    continue
                ext_map = self.calc_extinction(wavelength=(wavelength[-1] + wavelength[0])/2., xs=self.width,
                                               ys=self.height,
                                               extension_name=cur_ext)
                if ext_map is not None:
                    map_2d = map_2d * ext_map[0]
                continue

            all_flux_wl = [extname[-7:].strip('_') for extname in all_extensions
                           if extname is not None and (my_comp in extname and "FLUX_" in extname)]
            all_flux_wl = np.array(all_flux_wl)

            # flux_ext = [extname for extname in all_extensions
            #             if extname is not None and (my_comp in extname and
            #                                         "FLUX_{0}".format(np.round(wavelength, 2)) in extname)]
            add_emission = np.zeros(shape=(self.content[cur_ext].header['NAXIS2'],
                                           self.content[cur_ext].header['NAXIS1']), dtype=np.float32)
            if len(all_flux_wl) == 0:
                fluxrat_ext = [extname for extname in all_extensions
                               if extname is not None and (my_comp in extname and "FLUXRATIOS" in extname)]
                if len(fluxrat_ext) == 0:
                    continue
                fluxrat_ext = fluxrat_ext[0]

                wl_indexes = np.flatnonzero((self.content[fluxrat_ext].data[0, :] > (wavelength[0])) &
                                            (self.content[fluxrat_ext].data[0, :] < (wavelength[1])))

                # wl_index = np.flatnonzero(np.isclose(self.content[fluxrat_ext].data[0, :], wavelength))
                if len(wl_indexes) == 0:
                    continue
                for wl_index in wl_indexes:
                    add_emission += (self.content[cur_ext].data * self.content[fluxrat_ext].data[1, wl_index])
            else:
                all_flux_wl_float = np.array(all_flux_wl).astype(np.float32)
                wl_indexes = np.flatnonzero((all_flux_wl_float > (wavelength[0] - 0.01)) &
                                            (all_flux_wl_float < (wavelength[1] + 0.01)))
                flux_ext_wl = all_flux_wl[wl_indexes]
                if len(wl_indexes) == 0:
                    continue

                for cur_wl in flux_ext_wl:
                    add_emission += self.content[my_comp + "_FLUX_" + cur_wl].data

            map_2d[y0_in_field: y1_in_field + 1, x0_in_field: x1_in_field + 1] += \
                add_emission[y0_in_field - self.content[cur_ext].header['Y0']:
                             y1_in_field - self.content[cur_ext].header['Y0'] + 1,
                             x0_in_field - self.content[cur_ext].header['X0']: x1_in_field -
                             self.content[cur_ext].header['X0'] + 1
                             ]

            map_is_empty = False
        return map_2d * (proj_plane_pixel_scales(self.wcs)[0] * 3600) ** 2

    def get_spectrum(self, wl_grid=None, aperture_mask=None, fibers_coords=None, pix_size=1):
        if aperture_mask is None or (np.sum(aperture_mask) == 0) or (self.content[0].header['Nobj'] == 0):
            log.warning("No overlapping detected between the ISM component and the fibers => no spectra extraction")
            return None
        all_extensions = [hdu.header.get('EXTNAME') for hdu in self.content]
        all_extensions_brt = np.array([extname for extname in all_extensions
                                       if extname is not None and ("BRIGHTNESS" in extname) and
                                       check_overlap(self.content[extname], (self.height, self.width))])
        if all([self.content[cur_ext].header.get("DARK") for cur_ext in all_extensions_brt]):
            return None
        all_extensions_brt = all_extensions_brt[np.argsort([self.content[cur_ext].header.get('ZORDER')
                                                            for cur_ext in all_extensions_brt])]

        wl_logscale = np.log(wl_grid.value)
        wl_logscale_highres = np.arange((np.round(wl_logscale[-1] - wl_logscale[0]) * 1e6
                                         ).astype(int)) * 1e-6 + np.round(wl_logscale[0], 6)
        delta_lr = np.roll(wl_logscale, -1) - wl_logscale
        delta_lr[-1] = delta_lr[-2]

        xx, yy = np.meshgrid(np.arange(aperture_mask.shape[1]), np.arange(aperture_mask.shape[0]))
        pix_in_apertures = aperture_mask > 0
        xstart = np.clip(np.min(xx[pix_in_apertures]) - 1, 0, None)
        ystart = np.clip(np.min(yy[pix_in_apertures]) - 1, 0, None)
        xfin = np.clip(np.max(xx[pix_in_apertures]) + 1, None, aperture_mask.shape[1]-1)
        yfin = np.clip(np.max(yy[pix_in_apertures]) + 1, None, aperture_mask.shape[0]-1)
        aperture_mask_sub = aperture_mask[ystart: yfin + 1, xstart: xfin + 1]
        xx_sub, yy_sub = np.meshgrid(np.arange(xfin - xstart + 1), np.arange(yfin - ystart + 1))
        n_apertures = np.max(aperture_mask)
        # aperture_centers = np.round(fibers_coords).astype(int)
        spectrum = np.zeros(shape=(n_apertures, len(wl_grid)), dtype=np.float32)

        fiber_radius = np.median(fibers_coords[:, 2])
        kern_mask = calc_circular_mask(fiber_radius)
        kern = kernels.CustomKernel(kern_mask.reshape((1, kern_mask.shape[0], kern_mask.shape[1])))

        bar = progressbar.ProgressBar(max_value=len(all_extensions_brt)).start()
        for neb_index, cur_ext in enumerate(all_extensions_brt):
            cur_neb_in_mask = np.zeros_like(xx_sub)
            y0 = self.content[cur_ext].header.get("Y0")
            x0 = self.content[cur_ext].header.get("X0")
            nx = self.content[cur_ext].header.get('NAXIS1')
            ny = self.content[cur_ext].header.get('NAXIS2')
            cur_neb_in_mask[(xx_sub >= (x0 - xstart)) & (xx_sub <= (x0 + nx - xstart)) &
                            (yy_sub >= (y0 - ystart)) & (yy_sub <= (y0 + ny - ystart))] = True
            cur_neb_in_mask_ap = cur_neb_in_mask * (aperture_mask_sub > 0)
            if not np.sum(cur_neb_in_mask_ap):
                bar.update(neb_index + 1)
                continue
            cur_mask_in_neb = np.zeros(shape=self.content[cur_ext].data.shape, dtype=bool)
            xx_neb, yy_neb = np.meshgrid(np.arange(self.content[cur_ext].data.shape[1]),
                                         np.arange(self.content[cur_ext].data.shape[0]))
            cur_mask_in_neb[(xx_neb >= (xstart - x0)) & (xx_neb <= (xfin - x0)) &
                            (yy_neb >= (ystart - y0)) & (yy_neb <= (yfin - y0))] = True
            xstart_neb = np.min(xx_neb[cur_mask_in_neb])
            ystart_neb = np.min(yy_neb[cur_mask_in_neb])
            xfin_neb = np.max(xx_neb[cur_mask_in_neb])
            yfin_neb = np.max(yy_neb[cur_mask_in_neb])
            selected_apertures = np.flatnonzero(((fibers_coords[:, 0] + fibers_coords[:, 2]) >= x0) &
                                                ((fibers_coords[:, 0] - fibers_coords[:, 2]) <= (x0 + nx - 1)) &
                                                ((fibers_coords[:, 1] + fibers_coords[:, 2]) >= y0) &
                                                ((fibers_coords[:, 1] - fibers_coords[:, 2]) <= (y0 + ny - 1)))
            selected_apertures = np.array([sa for sa in selected_apertures if (sa+1) in aperture_mask_sub], dtype=int)
            if len(selected_apertures) == 0:
                bar.update(neb_index + 1)
                continue

            # Here I check if it is necessary to extend all involved arrays to account for fibers at the edges of neb.
            dx0 = abs(np.clip(np.round(np.min(fibers_coords[selected_apertures, 0] -
                                              fibers_coords[selected_apertures, 2])).astype(int) - 1 - x0, None, 0))
            dx1 = np.clip(np.round(np.max(fibers_coords[selected_apertures, 0] +
                                          fibers_coords[selected_apertures, 2])).astype(int) + 2 - x0 - nx, 0, None)
            dy0 = abs(np.clip(np.round(np.min(fibers_coords[selected_apertures, 1] -
                                              fibers_coords[selected_apertures, 2])).astype(int) - 1 - y0, None, 0))
            dy1 = np.clip(np.round(np.max(fibers_coords[selected_apertures, 1] +
                                          fibers_coords[selected_apertures, 2])).astype(int) + 2 - y0 - ny, 0, None)
            if (dx0 > 0) or (dy0 > 0) or (dx1 > 0) or (dy1 > 0):
                npad = ((0, 0), (dy0, dy1), (dx0, dx1))
            else:
                npad = None

            selected_y = fibers_coords[selected_apertures, 1] - ystart_neb - y0 + dy0
            selected_x = fibers_coords[selected_apertures, 0] - xstart_neb - x0 + dx0

            if self.content[cur_ext].header.get("DARK"):
                extinction_map = self.content[cur_ext].data[cur_mask_in_neb].reshape((1, yfin_neb - ystart_neb + 1,
                                                                                      xfin_neb - xstart_neb + 1))
                if npad is not None:
                    extinction_map = np.pad(extinction_map, pad_width=npad, mode='constant', constant_values=0)

                if self.content[cur_ext].header.get('EXT_LAW'):
                    cur_extinction_law = self.content[cur_ext].header.get('EXT_LAW')
                else:
                    cur_extinction_law = self.ext_law
                if self.content[cur_ext].header.get('EXT_RV'):
                    cur_r_v = self.content[cur_ext].header.get('EXT_RV')
                else:
                    cur_r_v = self.R_V

                data_in_apertures = convolve_array(extinction_map,
                                                   kern, selected_x, selected_y, pix_size,
                                                   nchunks=self.content[cur_ext].header.get("NCHUNKS"),
                                                   normalize_kernel=True)

                # data_in_apertures = convolve_fft(extinction_map, kern,
                #                                  allow_huge=True, normalize_kernel=True)[
                #                     :, aperture_centers[selected_apertures, 1] - ystart_neb - y0 + dy0,
                #                     aperture_centers[selected_apertures, 0] - xstart_neb - x0 + dx0]
                data_in_apertures = data_in_apertures.reshape((data_in_apertures.shape[0] * data_in_apertures.shape[1],
                                                               1))
                spectrum[selected_apertures, :] = \
                    spectrum[selected_apertures, :] * ism_extinction(av=data_in_apertures, r_v=cur_r_v,
                                                                     ext_law=cur_extinction_law, wavelength=wl_grid).T

                bar.update(neb_index + 1)
                continue

            my_comp = "_".join(cur_ext.split("_")[:2])

            if self.content[cur_ext].header.get('MAXBRT'):
                if self.content[cur_ext].header.get("LINERAT") == 'Variable':
                    all_wavelength = np.array([extname.split("_")[-1] for extname in all_extensions
                                               if extname is not None and (my_comp + "_FLUX_" in extname)])
                else:
                    all_wavelength = self.content[my_comp + "_FLUXRATIOS"].data[0, :]

                if self.content[cur_ext].header.get("LINERAT") == 'Variable':
                    all_fluxes = np.array([self.content[my_comp + "_FLUX_" + wl].data[
                                               cur_mask_in_neb].reshape((yfin_neb - ystart_neb + 1,
                                                                         xfin_neb - xstart_neb + 1))
                                           for wl in all_wavelength], dtype=np.float32)
                    all_wavelength = all_wavelength.astype(np.float32)
                else:
                    all_fluxes = self.content[cur_ext].data[cur_mask_in_neb].reshape((1, yfin_neb - ystart_neb + 1,
                                                                                      xfin_neb - xstart_neb + 1)
                                                                                     ).astype(np.float32)

                if my_comp + "_LINEPROFILE" in self.content:
                    lsf = self.content[my_comp + "_LINEPROFILE"].data[
                          :, cur_mask_in_neb].reshape((len(self.vel_grid), yfin_neb - ystart_neb + 1,
                                                       xfin_neb - xstart_neb + 1)).astype(np.float32)
                else:
                    if my_comp + "_VEL" in self.content:
                        vel = self.content[my_comp + "_VEL"].data[cur_mask_in_neb].reshape((yfin_neb - ystart_neb + 1,
                                                                                            xfin_neb - xstart_neb + 1))
                    else:
                        vel = np.zeros(shape=(yfin_neb - ystart_neb + 1, xfin_neb - xstart_neb + 1),
                                       dtype=np.float32)  # + self.content[cur_ext].header.get('SysVel')
                        # self.sys_velocity.value
                    if my_comp + "_DISP" in self.content:
                        disp = self.content[my_comp + "_DISP"].data[
                            cur_mask_in_neb].reshape((yfin_neb - ystart_neb + 1, xfin_neb - xstart_neb + 1))
                    else:
                        disp = np.zeros(shape=(yfin_neb - ystart_neb + 1, xfin_neb - xstart_neb + 1),
                                        dtype=np.float32) + self.content[cur_ext].header.get('TurbVel')
                        # self.turbulent_sigma.value
                    lsf = np.exp(-np.power(
                        (self.vel_grid.value[:, None, None] - vel[None, :, :]) / disp[None, :, :], 2.) / 2)
                    lsf = (lsf / np.nansum(lsf, axis=0)).astype(np.float32)

                if npad is not None:
                    lsf = np.pad(lsf, pad_width=npad, mode='constant', constant_values=0)
                    all_fluxes = np.pad(all_fluxes, pad_width=npad, mode='constant', constant_values=0)

                data_in_apertures = [convolve_array(lsf * line_data[None, :, :],
                                                    kern, selected_x, selected_y, pix_size,
                                                    nchunks=self.content[cur_ext].header.get("NCHUNKS"))
                                     for line_data in all_fluxes]
                data_in_apertures = np.array(data_in_apertures)
                if data_in_apertures.shape == 2:
                    data_in_apertures = data_in_apertures.reshape((1, data_in_apertures.shape[0],
                                                                   data_in_apertures.shape[1]))

                # if all_fluxes.shape[0] == 1:
                #     # selected_y = fibers_coords[selected_apertures, 1] - ystart_neb - y0 + dy0
                #     # selected_x = fibers_coords[selected_apertures, 0] - xstart_neb - x0 + dx0
                #     data_in_apertures = convolve_array(lsf * all_fluxes[0][None, :, :],
                #                                        kern, selected_x, selected_y, pix_size,
                #                                        nchunks=self.content[cur_ext].header.get("NCHUNKS"))
                #     data_in_apertures = data_in_apertures.reshape((1, data_in_apertures.shape[0],
                #                                                   data_in_apertures.shape[1]))
                #     # hdu = fits.PrimaryHDU(data=data_in_apertures)
                #     # hdu.writeto(f'./data_{neb_index}.fits', overwrite=True)
                #
                # else:
                #
                #     # selected_y = fibers_coords[selected_apertures, 1] - ystart_neb - y0 + dy0
                #     # selected_x = fibers_coords[selected_apertures, 0] - xstart_neb - x0 + dx0
                #     # data_in_apertures = Parallel(n_jobs=lvmdatasimulator.n_process)(
                #     #     delayed(convolve_array)(lsf * line_data[None, :, :], kern,
                #     #                            selected_y, selected_x, pix_size)
                #     #     for line_data in all_fluxes)
                #     data_in_apertures = [convolve_array(lsf * line_data[None, :, :], kern,
                #                                         selected_y, selected_x, pix_size,
                #                                         nchunks=self.content[cur_ext].header.get("NCHUNKS"))
                #                          for line_data in all_fluxes]
                #     data_in_apertures = np.array(data_in_apertures)
                #     # hdu = fits.PrimaryHDU(data=data_in_apertures)
                #     # hdu.writeto(f'./data_{neb_index}.fits', overwrite=True)

                data_in_apertures = np.moveaxis(data_in_apertures, 2, 0)
                if data_in_apertures.shape[1] > 1:
                    prf_index = np.flatnonzero(all_wavelength == 6562.81)
                else:
                    prf_index = 0
                flux_norm_in_apertures = data_in_apertures.sum(axis=2)
                line_prf_in_apertures = data_in_apertures[:, prf_index, :].reshape(
                    (data_in_apertures.shape[0], data_in_apertures.shape[2])) / \
                    flux_norm_in_apertures[:, prf_index].reshape(data_in_apertures.shape[0], 1)

                wl_logscale_lsf = np.log((self.vel_grid.value + self.content[cur_ext].header.get('SysVel')
                                          ) / 2.9979e5 + 1)
                highres_factor = 1e6
                wl_logscale_lsf_highres = np.arange(np.round((wl_logscale_lsf[-1] - wl_logscale_lsf[0]
                                                              ) / 2. * highres_factor).astype(int) * 2 + 1
                                                    ) / highres_factor + wl_logscale_lsf[0]
                p = interp1d(wl_logscale_lsf, line_prf_in_apertures, axis=1, assume_sorted=True)
                line_highres_log = p(wl_logscale_lsf_highres)
                line_highres_log = line_highres_log / np.sum(line_highres_log, axis=1)[:, None]
                if flux_norm_in_apertures.shape[1] == 1:
                    flux_norm_in_apertures = flux_norm_in_apertures * \
                                        self.content[my_comp + "_FLUXRATIOS"].data[1, None, :]
                wl_indexes = np.round((np.log(all_wavelength * (1 + self.content[cur_ext].header.get('SysVel') /
                                                                2.9979e5)) - wl_logscale_highres[0]
                                       ) * highres_factor).astype(int)
                rec = (wl_indexes > 0) & (wl_indexes < len(wl_logscale_highres))
                spectrum_highres_log = np.zeros(shape=(len(selected_apertures), len(wl_logscale_highres)),
                                                dtype=np.float32)
                win = (len(wl_logscale_lsf_highres) - 1) // 2
                for ind, r in enumerate(rec):
                    if r:
                        spectrum_highres_log[:, wl_indexes[ind] - win: wl_indexes[ind] + win + 1] += \
                            line_highres_log[:, :] * flux_norm_in_apertures[:, ind].reshape(
                                (len(selected_apertures), 1))
                p = interp1d(wl_logscale_highres, spectrum_highres_log, axis=1, assume_sorted=True, bounds_error=False,
                             fill_value='extrapolate')
                spectrum[selected_apertures, :] += (p(wl_logscale) * delta_lr * highres_factor)

            if my_comp + "_CONTINUUM" in self.content:
                brt_max = self.content[cur_ext].header.get('MAXBRT')
                if not brt_max:
                    brt_max = self.content[cur_ext].header.get('MAXEXT')
                if not brt_max:
                    brt_max = 1
                continuum = self._get_continuum(my_comp, wl_grid.value)

                brt = self.content[cur_ext].data[cur_mask_in_neb].reshape((1, yfin_neb - ystart_neb + 1,
                                                                           xfin_neb - xstart_neb + 1)) / brt_max
                if brt.shape[0] != 1:
                    brt = brt.reshape((1, brt.shape[0], brt.shape[1]))
                if npad is not None:
                    brt = np.pad(brt, pad_width=npad, mode='constant', constant_values=0)
                # selected_y = fibers_coords[selected_apertures, 1] - ystart_neb - y0 + dy0
                # selected_x = fibers_coords[selected_apertures, 0] - xstart_neb - x0 + dx0
                data_in_apertures = convolve_array(brt,
                                                   kern, selected_x, selected_y, pix_size,
                                                   nchunks=self.content[cur_ext].header.get("NCHUNKS"))
                # data_in_apertures = convolve_fft(brt, kern,
                #                                  allow_huge=True, normalize_kernel=False)[
                #                     :, aperture_centers[selected_apertures, 1] - ystart_neb - y0 + dy0,
                #                     aperture_centers[selected_apertures, 0] - xstart_neb - x0 + dx0]
                data_in_apertures = data_in_apertures.reshape((data_in_apertures.shape[0] * data_in_apertures.shape[1],
                                                               1))
                spectrum[selected_apertures, :] += continuum[None, :] * data_in_apertures

            bar.update(neb_index + 1)
        bar.finish()
        return spectrum * (proj_plane_pixel_scales(self.wcs)[0] * 3600) ** 2 * fluxunit * u.arcsec ** 2

# encoding: utf-8
#
# @Author: <NAME>, <NAME>
# @Date: Nov 15, 2021
# @Filename: ism.py
# @License: BSD 3-Clause
# @Copyright: <NAME>, <NAME>

import os.path
from astropy import units as u
from astropy import constants as c
import numpy as np
from astropy.io import fits, ascii
from astropy.table import Table
from scipy.special import sph_harm
from astropy.wcs import WCS
from astropy.wcs.utils import proj_plane_pixel_scales
from astropy.coordinates import SkyCoord
from astropy.modeling.models import Sersic2D
from dataclasses import dataclass
import sys
if (sys.version_info[0]+sys.version_info[1]/10.) < 3.8:
    from backports.cached_property import cached_property
else:
    from functools import cached_property
from scipy.ndimage.interpolation import map_coordinates
from scipy.interpolate import interp1d, interp2d
import lvmdatasimulator
from lvmdatasimulator import log
import progressbar
from joblib import Parallel, delayed
from astropy.convolution import convolve_fft, kernels
from lvmdatasimulator.utils import calc_circular_mask, convolve_array, set_default_dict_values, \
    ism_extinction, check_overlap, assign_units
fluxunit = u.erg / (u.cm ** 2 * u.s * u.arcsec ** 2)
velunit = u.km / u.s


def brightness_inhomogeneities_sphere(harm_amplitudes, ll, phi_cur, theta_cur, rho, med, radius, thickness):
    """
    Auxiliary function producing the inhomogeneities on the brightness distribution for the Cloud of Bubble objects
    using the spherical harmonics.
    """
    brt = theta_cur * 0
    for m in np.arange(-ll, ll + 1):
        brt += (harm_amplitudes[m + ll * (ll + 1) - 1] * sph_harm(m, ll, phi_cur, theta_cur).real * med *
                (1 - np.sqrt(abs(rho.value ** 2 / radius.value ** 2 - (1 - thickness / 2) ** 2))))
    return brt


def sphere_brt_in_line(brt_3d, rad_3d, rad_model, flux_model):
    """
    Auxiliary function computing the brightness of the Cloud or Bubble at given radii and in given line
    according to the Cloudy models
    """
    p = interp1d(rad_model, flux_model, fill_value='extrapolate', assume_sorted=True)
    return p(rad_3d) * brt_3d


def interpolate_sphere_to_cartesian(spherical_array, x_grid=None, y_grid=None, z_grid=None,
                                    rad_grid=None, theta_grid=None, phi_grid=None, pxscale=1. * u.pc):
    """
    Auxiliary function to project the brightness or velocities from the spherical to cartesian coordinates
    """
    x, y, z = np.meshgrid(x_grid, y_grid, z_grid, indexing='ij')
    phi_c, theta_c, rad_c = xyz_to_sphere(x, y, z, pxscale=pxscale)
    ir = interp1d(rad_grid, np.arange(len(rad_grid)), bounds_error=False)
    ith = interp1d(theta_grid, np.arange(len(theta_grid)))
    iphi = interp1d(phi_grid, np.arange(len(phi_grid)))
    new_ir = ir(rad_c.ravel())
    new_ith = ith(theta_c.ravel())
    new_iphi = iphi(phi_c.ravel())
    cart_data = map_coordinates(spherical_array, np.vstack([new_ir, new_ith, new_iphi]),
                                order=1, mode='constant', cval=0)
    return cart_data.reshape([len(x_grid), len(y_grid), len(z_grid)]).T


def limit_angle(value, bottom_limit=0, top_limit=np.pi):
    """
    Auxiliary function to limit the angle values to the range of [0, pi]
    """
    value[value < bottom_limit] += (top_limit - bottom_limit)
    value[value > top_limit] -= (top_limit - bottom_limit)
    return value


def xyz_to_sphere(x, y, z, pxscale=1. * u.pc):
    """
    Auxiliary function to map the coordinates from cartesian to spherical system
    """
    phi_c = np.arctan2(y, x)
    rad_c = (np.sqrt(x ** 2 + y ** 2 + z ** 2))
    rad_c[rad_c == 0 * u.pc] = 1e-3 * pxscale
    theta_c = (np.arccos(z / rad_c))
    phi_c = limit_angle(phi_c, 0 * u.radian, 2 * np.pi * u.radian)
    theta_c = limit_angle(theta_c, 0 * u.radian, np.pi * u.radian)
    return phi_c, theta_c, rad_c


def find_model_id(file=lvmdatasimulator.CLOUDY_MODELS,
                  check_id=None, params=lvmdatasimulator.CLOUDY_SPEC_DEFAULTS['id']):
    """
    Checks the input parameters of the pre-computed Cloudy model and return corresponding index in the grid
    """
    with fits.open(file) as hdu:
        if check_id is None:
            if params is None:
                check_id = lvmdatasimulator.CLOUDY_SPEC_DEFAULTS['id']
                log.warning(f'Default Cloudy model will be used (id = {check_id})')
            else:
                summary_table = Table(hdu['Summary'].data)
                indexes = np.arange(len(summary_table)).astype(int)
                rec_table = np.ones(shape=len(summary_table), dtype=bool)

                def closest(rec, prop, val):
                    unique_col = np.unique(summary_table[prop][rec])
                    if isinstance(val, str):
                        res = unique_col[unique_col == val]
                        if len(res) == 0:
                            return ""
                        return res
                    else:
                        return unique_col[np.argsort(np.abs(unique_col - val))[0]]

                for p in params:
                    if p not in summary_table.colnames or params[p] is None or \
                            ((isinstance(params[p], float) or isinstance(params[p], int)) and ~np.isfinite(params[p])):
                        continue
                    rec_table = rec_table & (summary_table[p] == closest(indexes, p, params[p]))
                    indexes = np.flatnonzero(rec_table)
                    if len(indexes) == 0:
                        break
                if len(indexes) == 0 or len(indexes) == len(summary_table):
                    check_id = lvmdatasimulator.CLOUDY_SPEC_DEFAULTS['id']
                    log.warning('Input parameters do not correspond to any pre-computed Cloudy model.'
                                'Default Cloudy model will be used (id = {0})'.format(check_id))
                elif len(indexes) == 1:
                    check_id = summary_table['Model_ID'][indexes[0]]
                    for p in params:
                        if p not in summary_table.colnames or params[p] is None or \
                                ((isinstance(params[p], float) or
                                  isinstance(params[p], int)) and ~np.isfinite(params[p])):
                            continue
                        if params[p] != summary_table[p][indexes[0]]:
                            log.warning(f'Use the closest pre-computed Cloudy model with id = {check_id}')
                            break
                else:
                    check_id = summary_table['Model_ID'][indexes[0]]
                    log.warning(f'Select one of the closest pre-computed Cloudy model with id = {check_id}')
                #
                # for cur_ext in range(len(hdu)):
                #     if cur_ext == 0:
                #         continue
                #     found = False
                #     for p in params:
                #         if p == 'id':
                #             continue
                #         precision = 1
                #         if p == 'Z':
                #             precision = 2
                #         if np.round(params[p], precision) != np.round(hdu[cur_ext].header[p], precision):
                #             break
                #     else:
                #         found = True
                #     if found:
                #         return cur_ext, check_id
                # check_id = lvmdatasimulator.CLOUDY_SPEC_DEFAULTS['id']
                # log.warning('Input parameters do not correspond to any pre-computed Cloudy model.'
                #             'Default Cloudy model will be used (id = {0})'.format(check_id))

        extension_index = None
        while extension_index is None:
            extension_index = [cur_ext for cur_ext in range(len(hdu)) if (
                    check_id == hdu[cur_ext].header.get('MODEL_ID'))]
            if len(extension_index) == 0:
                if check_id == lvmdatasimulator.CLOUDY_SPEC_DEFAULTS['id']:
                    log.warning('Model_ID = {0} is not found in the Cloudy models grid. '
                                'Use the first one in the grid instead'.format(check_id))
                    extension_index = 1
                else:
                    log.warning('Model_ID = {0} is not found in the Cloudy models grid. '
                                'Use default ({1}) instead'.format(check_id,
                                                                   lvmdatasimulator.CLOUDY_SPEC_DEFAULTS['id']))
                    check_id = lvmdatasimulator.CLOUDY_SPEC_DEFAULTS['id']
                    extension_index = None
            else:
                extension_index = extension_index[0]
    return extension_index, check_id


@dataclass
class Nebula:
    """
    Base class defining properties of every nebula type.
    By itself it describes the rectangular nebula (e.g. DIG)
    Constructed nebula has 4 dimensions, where 4th derive its appearance in different lines
        (if spectrum_id is None, or if it is dark nebula => only one line)
    """
    xc: int = None  # Center of the region in the field of view, pix
    yc: int = None  # Center of the region in the field of view, pix
    x0: int = 0  # Coordinates of the bottom-left corner in the field of view, pix
    y0: int = 0  # Coordinates of the bottom-left corner in the field of view, pix
    pix_width: int = None  # full width of cartesian grid, pix (should be odd)
    pix_height: int = None  # full height of cartesian grid, pix (should be odd)
    width: u.pc = 0 * u.pc  # width of the nebula in pc (not used if pix_width is set up)
    height: u.pc = 0 * u.pc  # height of the nebula in pc (not used if pix_height is set up)
    pxscale: u.pc = 0.01 * u.pc  # pixel size in pc
    spectrum_id: int = None  # ID of a template Cloudy emission spectrum for this nebula
    n_brightest_lines: int = None  # limit the number of the lines to the first N brightest
    sys_velocity: velunit = 0 * velunit  # Systemic velocity
    turbulent_sigma: velunit = 10 * velunit  # Velocity dispersion due to turbulence; included in calculations of LSF
    max_brightness: fluxunit = 1e-15 * fluxunit
    max_extinction: u.mag = 0 * u.mag
    perturb_scale: int = 0 * u.pc  # Spatial scale of correlated perturbations
    perturb_amplitude: float = 0.1  # Maximal amplitude of perturbations
    _npix_los: int = 1  # full size along line of sight in pixels
    nchunks: int = -1  # number of chuncks to use for the convolution. If negative, select automatically
    vel_gradient: (velunit / u.pc) = 0  # velocity gradient along the nebula
    vel_pa: u.degree = 0  # Position angle of the kinematical axis (for the velocity gradient or rotation velocity)

    def __post_init__(self):
        self._assign_all_units()
        self._assign_position_params()
        self._ref_line_id = 0
        self.linerat_constant = True  # True if the ratio of line fluxes shouldn't change across the nebula

    def _assign_all_units(self):
        whole_list_properties = ['pxscale', 'sys_velocity', 'turbulent_sigma', 'max_brightness', 'max_extinction',
                                 'perturb_scale', 'radius', 'PA', 'length', 'width', 'vel_gradient', 'r_eff',
                                 'vel_rot', 'expansion_velocity', 'spectral_axis', 'vel_pa']
        whole_list_units = [u.pc, velunit, velunit, fluxunit, u.mag, u.pc, u.pc, u.degree, u.pc, u.pc,
                            (velunit / u.pc), u.kpc, velunit, velunit, velunit, u.degree]
        cur_list_properties = []
        cur_list_units = []
        for prp, unit in zip(whole_list_properties, whole_list_units):
            if hasattr(self, prp):
                cur_list_properties.append(prp)
                cur_list_units.append(unit)
        assign_units(self, cur_list_properties, cur_list_units)

    def _assign_position_params(self, conversion_type='rect'):
        if conversion_type == 'rect':
            for v in ['height', 'width']:
                if self.__getattribute__(f'pix_{v}') is None:
                    val = np.round((self.__getattribute__(v) / self.pxscale).value / 2.).astype(int) * 2 + 1
                else:
                    val = np.round(self.__getattribute__(f'pix_{v}') / 2.).astype(int) * 2 + 1
                setattr(self, f'pix_{v}', val)
        elif conversion_type == 'ellipse':
            self.pix_width = (np.round(np.abs(self.radius / self.pxscale * np.sin(self.PA)) +
                              np.abs(self.radius / self.pxscale *
                                     self.ax_ratio * np.cos(self.PA))).astype(int) * 2 + 1).value
            self.pix_height = (np.round(np.abs(self.radius / self.pxscale * np.cos(self.PA)) +
                               np.abs(self.radius / self.pxscale *
                                      self.ax_ratio * np.sin(self.PA))).astype(int) * 2 + 1).value
        elif conversion_type == 'galaxy':
            self.pix_width = (np.round(np.abs(self.r_max * np.sin(self.PA)) +
                                       np.abs(self.r_max * self.ax_ratio * np.cos(self.PA))).astype(int) * 2 + 1).value
            self.pix_height = (np.round(np.abs(self.r_max * np.cos(self.PA)) +
                                        np.abs(self.r_max * self.ax_ratio * np.sin(self.PA))).astype(int) * 2 + 1).value
        elif conversion_type == 'cylinder':
            self.pix_width = (np.ceil((self.length * np.abs(np.sin(self.PA)) +
                                       self.width * np.abs(np.cos(self.PA))) / self.pxscale / 2.
                                      ).astype(int) * 2 + 1).value
            self.pix_height = (np.ceil((self.length * np.abs(np.cos(self.PA)) +
                                       self.width * np.abs(np.sin(self.PA))) / self.pxscale / 2.
                                       ).astype(int) * 2 + 1).value

        if (self.xc is not None) and (self.yc is not None):
            self.x0 = self.xc - np.round((self.pix_width - 1) / 2).astype(int)
            self.y0 = self.yc - np.round((self.pix_height - 1) / 2).astype(int)
        elif (self.x0 is not None) and (self.y0 is not None):
            self.xc = self.x0 + np.round((self.pix_width - 1) / 2).astype(int)
            self.yc = self.y0 + np.round((self.pix_height - 1) / 2).astype(int)

    @cached_property
    def _cartesian_x_grid(self):
        return np.arange(self.pix_width) * self.pxscale

    @cached_property
    def _cartesian_y_grid(self):
        return np.arange(self.pix_height) * self.pxscale

    @cached_property
    def _cartesian_z_grid(self):
        return np.arange(self._npix_los) * self.pxscale

    @cached_property
    def _max_density(self):
        return self.max_extinction * (1.8e21 / (u.cm ** 2 * u.mag))

    @cached_property
    def _brightness_3d_cartesian(self):
        """
        Method to obtain the brightness (or density) distribution of the nebula in cartesian coordinates
        """
        brt = np.ones(shape=(self.pix_height, self.pix_width, self._npix_los), dtype=float) / self._npix_los
        if (self.perturb_scale > 0) and (self.perturb_amplitude > 0):
            pertscale = (self.perturb_scale / self.pxscale).value
            perturb = np.random.uniform(-1, 1, (self.pix_height, self.pix_width)
                                        ) * self.perturb_amplitude / self._npix_los
            xx, yy = np.meshgrid(np.arange(self.pix_width), np.arange(self.pix_height))
            f = np.exp(-2 * (xx ** 2 + yy ** 2) / pertscale)
            perturb = 4 / np.sqrt(np.pi) / pertscale * np.fft.ifft2(np.fft.fft2(perturb) * np.fft.fft2(f)).real
            brt += (perturb[:, :, None] - np.median(perturb))
        return brt

    @cached_property
    def _brightness_4d_cartesian(self):
        """
        Derive the brightness (or density) distribution of the nebula for each emission line in cartesian coordinates
        """
        if self.spectrum_id is None or self.linerat_constant:
            flux_ratios = np.array([1.])
        else:
            with fits.open(lvmdatasimulator.CLOUDY_MODELS) as hdu:
                flux_ratios = hdu[self.spectrum_id].data[1:, 1]
                index_ha = np.flatnonzero(hdu[self.spectrum_id].data[1:, 0] == 6562.81)
                if self.n_brightest_lines is not None and \
                        (self.n_brightest_lines > 0) and (self.n_brightest_lines < len(flux_ratios)):
                    indexes_sorted = np.argsort(flux_ratios)[::-1]
                    flux_ratios = flux_ratios[indexes_sorted[: self.n_brightest_lines]]
                    index_ha = np.flatnonzero(hdu[self.spectrum_id].data[1:, 0][indexes_sorted] == 6562.81)
                if len(index_ha) == 1:
                    self._ref_line_id = index_ha[0]
        return self._brightness_3d_cartesian[None, :, :, :] * flux_ratios[:, None, None, None]

    @cached_property
    def brightness_skyplane(self):
        """
        Project the 3D nebula onto sky plane (for emission or continuum sources)
        """
        if self.max_brightness > 0:
            norm_max = self.max_brightness
        else:
            norm_max = 1
        map2d = np.nansum(self._brightness_3d_cartesian, 2)
        return map2d / np.nanmax(map2d) * norm_max

    @cached_property
    def brightness_skyplane_lines(self):
        """
        Project the 3D emission nebula line onto sky plane (return images in each emission line)
        """
        if self.max_brightness > 0:
            map2d = np.nansum(self._brightness_4d_cartesian, 3)
            return map2d / np.nanmax(map2d[self._ref_line_id, :, :]) * self.max_brightness
        else:
            return None

    @cached_property
    def extinction_skyplane(self):
        """
        Project the 3D nebula onto sky plane (for dark clouds)
        """
        if self.max_extinction > 0:
            map2d = np.nansum(self._brightness_3d_cartesian, 2)
            return map2d / np.nanmax(map2d) * self._max_density / (1.8e21 / (u.cm ** 2 * u.mag))
        else:
            return None

    @cached_property
    def vel_field(self):
        return self._get_2d_velocity()
        # if vel_field is None:
        #     return np.atleast_1d(self.sys_velocity)
        # else:
        #     return vel_field + self.sys_velocity

    def _get_2d_velocity(self):
        if hasattr(self, 'vel_gradient') and (self.vel_gradient is not None) and (self.vel_gradient != 0):
            xx, yy = np.meshgrid(np.arange(self.pix_width), np.arange(self.pix_height))
            vel_field = (- (xx - (self.pix_width - 1) / 2) * np.sin(self.vel_pa) +
                         (yy - (self.pix_height - 1) / 2) * np.cos(self.vel_pa)) * self.pxscale * self.vel_gradient
            return vel_field
        else:
            return None

    # @cached_property
    # def line_profile(self):
    #     lprf = np.zeros(shape=len(self.los_velocity), dtype=float)
    #     lprf[np.floor(len(lprf) / 2.).astype(int)] = 1.
    #     return lprf


@dataclass
class Rectangle(Nebula):
    """
    Class defining a simple rectangular component.
    This is equal to Nebula, but no perturbations and turbulence by default
    """
    perturb_amplitude: float = 0.0  # Maximal amplitude of perturbations
    turbulent_sigma: velunit = 0 * velunit  # Velocity dispersion due to turbulence; included in calculations of LSF

    def __post_init__(self):
        self._assign_all_units()
        self._assign_position_params()
        self._ref_line_id = 0
        self.linerat_constant = True  # True if the ratio of line fluxes shouldn't change across the nebula


@dataclass
class Ellipse(Nebula):
    """
    Class defining a simple elliptical component.
    No perturbations and turbulence by default
    """
    perturb_amplitude: float = 0.0  # Maximal amplitude of perturbations
    turbulent_sigma: velunit = 0 * velunit  # Velocity dispersion due to turbulence; included in calculations of LSF
    radius: u.pc = 1.0 * u.pc  # Radius along the major axis of the ellipse (or radius of the circle)
    PA: u.degree = 90 * u.degree  # position angle of the major axis
    ax_ratio: float = 1.  # ratio of minor/major axes

    def __post_init__(self):
        self._assign_all_units()
        self._npix_los = 1
        self._assign_position_params(conversion_type='ellipse')
        self._ref_line_id = 0
        self.linerat_constant = True  # True if the ratio of line fluxes shouldn't change across the nebula

    @cached_property
    def _brightness_3d_cartesian(self):
        """
        Method to obtain the brightness (or density) distribution of the nebula in cartesian coordinates
        """
        xx, yy = np.meshgrid(np.arange(self.pix_width), np.arange(self.pix_height))
        brt = np.ones(shape=(self.pix_height, self.pix_width), dtype=np.float32)
        angle = (self.PA + 90 * u.degree).to(u.radian).value
        xct = (xx - (self.pix_width - 1) / 2) * np.cos(angle) + \
              (yy - (self.pix_height - 1) / 2) * np.sin(angle)
        yct = (xx - (self.pix_width - 1) / 2) * np.sin(angle) - \
              (yy - (self.pix_height - 1) / 2) * np.cos(angle)
        rmaj = (self.radius.to(u.pc) / self.pxscale.to(u.pc)).value
        rmin = (self.radius.to(u.pc) / self.pxscale.to(u.pc)).value * self.ax_ratio
        rec = (xct ** 2 / rmaj ** 2) + (yct ** 2 / rmin ** 2) >= 1
        brt[rec] = 0
        brt = brt.reshape((self.pix_height, self.pix_width, 1))
        return brt


@dataclass
class Circle(Ellipse):
    """
    Class defining a simple circular component.
    """
    def __post_init__(self):
        self._assign_all_units()
        self.ax_ratio = 1.
        self._npix_los = 1
        self._assign_position_params(conversion_type='ellipse')
        self._ref_line_id = 0
        self.linerat_constant = True  # True if the ratio of line fluxes shouldn't change across the nebula


@dataclass
class Filament(Nebula):
    """
    Class of an isotropic cylindrical shape filament.
    Defined by its position, length, PA, radius, maximal optical depth.
    If it is emission-type filament, then also maximal brightness is required.
    Velocity gradient also can be set up
    """
    PA: u.degree = 90 * u.degree  # position angle of the filament
    length: u.pc = 10 * u.pc  # full length of the filament
    width: u.pc = 0.1 * u.pc  # full width (diameter) of the filament

    def __post_init__(self):
        self._assign_all_units()
        self._assign_position_params(conversion_type='cylinder')
        self._npix_los = 1
        self._ref_line_id = 0
        self.linerat_constant = True  # True if the ratio of line fluxes shouldn't change across the nebula

    @cached_property
    def _brightness_3d_cartesian(self):
        """
        Method to obtain the brightness (or density) distribution of the nebula in cartesian coordinates
        """
        xx, yy = np.meshgrid(np.arange(self.pix_width), np.arange(self.pix_height))
        brt = np.zeros_like(xx, dtype=np.float32)

        xct = (xx - (self.pix_width - 1) / 2) * np.cos(self.PA + 90 * u.degree) + \
              (yy - (self.pix_height - 1) / 2) * np.sin(self.PA + 90 * u.degree)
        yct = (xx - (self.pix_width - 1) / 2) * np.sin(self.PA + 90 * u.degree) - \
              (yy - (self.pix_height - 1) / 2) * np.cos(self.PA + 90 * u.degree)

        rad = ((self.width / self.pxscale).value / 2.)
        len_px = ((self.length / self.pxscale).value / 2.)
        rec = (np.abs(yct) <= rad) & (np.abs(xct) <= len_px)
        brt[rec] = np.sqrt(1. - (yct[rec] / rad) ** 2)
        brt = brt.reshape((self.pix_height, self.pix_width, 1))
        return brt


@dataclass
class _ObsoleteFilament(Nebula):
    """
    Class of an isotropic cylindrical shape filament.
    Defined by its position, length, PA, radius, maximal optical depth
    if it is emission-type filament, then maximal brightness
    NB: this class is obsolete, but might be considered later in case of implementation of varying line ratios
    """
    PA: u.degree = 90 * u.degree  # position angle of the filament
    length: u.pc = 10 * u.pc  # full length of the filament
    width: u.pc = 0.1 * u.pc  # full width (diameter) of the filament
    vel_gradient: (velunit / u.pc) = 0  # velocity gradient along the filament (to be added)
    _theta_bins: int = 50
    _rad_bins: int = 0
    _h_bins: int = 2
    _npix_los: int = 101

    def __post_init__(self):
        self._assign_all_units()
        if self._rad_bins == 0:
            self._rad_bins = np.ceil(self.width.to(u.pc).value / self.pxscale.to(u.pc).value * 5).astype(int)
        if (self.xc is not None) and (self.yc is not None):
            self.x0 = self.xc - np.round((len(self._cartesian_y_grid) - 1) / 2).astype(int)
            self.y0 = self.yc - np.round((len(self._cartesian_z_grid) - 1) / 2).astype(int)
        elif (self.x0 is not None) and (self.y0 is not None):
            self.xc = self.x0 + np.round((len(self._cartesian_y_grid) - 1) / 2).astype(int)
            self.yc = self.y0 + np.round((len(self._cartesian_z_grid) - 1) / 2).astype(int)
        self._ref_line_id = 0
        self.linerat_constant = True  # True if the ratio of line fluxes shouldn't change across the nebula

    @cached_property
    def _theta_grid(self):
        return np.linspace(0, 2 * np.pi, self._theta_bins)

    @cached_property
    def _h_grid(self):
        return np.linspace(0, self.length, self._h_bins)

    @cached_property
    def _rad_grid(self):
        return np.linspace(0, self.width / 2, self._rad_bins)

    @cached_property
    def _cartesian_y_grid(self):
        npix = np.ceil(1.01 * (self.length * np.abs(np.sin(self.PA)) +
                               self.width * np.abs(np.cos(self.PA))) / self.pxscale).astype(int)
        npix_l = npix / 2 - np.ceil(self.length / 2 * np.sin(-self.PA) / self.pxscale).astype(int)
        return (np.linspace(0, npix, npix + 1) - npix_l) * self.pxscale

    @cached_property
    def _cartesian_z_grid(self):
        npix = np.ceil(1.01 * (self.length * np.abs(np.cos(self.PA)) +
                               self.width * np.abs(np.sin(self.PA))) / self.pxscale).astype(int)
        npix_l = npix / 2 - np.ceil(self.length / 2 * np.cos(-self.PA) / self.pxscale).astype(int)
        return (np.linspace(0, npix, npix + 1) - npix_l) * self.pxscale

    @cached_property
    def _cartesian_x_grid(self):
        return np.linspace(-1.01, 1.01, self._npix_los) * self.width / 2

    @cached_property
    def _brightness_3d_cylindrical(self):
        """
        Method to calculate brightness (or opacity) of the cloud at given theta, phi and radii

        theta: float -- azimuthal angle [0, 2 * np.pi]
        rad: float -- radius [0, self.width / 2]
        h: float -- height [0, self.length]
        Returns:
            3D cube of normalized brightness in theta-rad-h grid; total brightness = 1
        """
        rho, theta, h = np.meshgrid(self._rad_grid, self._theta_grid, self._h_grid, indexing='ij')
        brt = np.ones_like(theta)
        brt[rho > (self.width / 2)] = 0
        norm = np.sum(brt)
        if norm > 0:
            brt = brt / np.sum(brt)
        return brt

    @cached_property
    def _brightness_3d_cartesian(self):
        x, y, z = np.meshgrid(self._cartesian_x_grid, self._cartesian_y_grid,
                              self._cartesian_z_grid, indexing='ij')

        h_c = -y * np.sin(self.PA) + z * np.cos(self.PA)
        theta_c = np.arctan2(y * np.cos(self.PA) + z * np.sin(self.PA), x)
        rad_c = np.sqrt(x ** 2 + (y * np.cos(self.PA) + z * np.sin(self.PA)) ** 2)
        rad_c[rad_c == 0 * u.pc] = 1e-3 * self.pxscale

        theta_c = limit_angle(theta_c, 0 * u.radian, 2 * np.pi * u.radian)

        ir = interp1d(self._rad_grid, np.arange(self._rad_bins), bounds_error=False)
        ith = interp1d(self._theta_grid, np.arange(self._theta_bins))
        ih = interp1d(self._h_grid, np.arange(self._h_bins), bounds_error=False)
        new_ir = ir(rad_c.ravel())
        new_ith = ith(theta_c.ravel())
        new_ih = ih(h_c.ravel())

        cart_data = map_coordinates(self._brightness_3d_cylindrical,
                                    np.vstack([new_ir, new_ith, new_ih]),
                                    order=1, mode='constant', cval=0)

        return cart_data.reshape([len(self._cartesian_x_grid),
                                  len(self._cartesian_y_grid),
                                  len(self._cartesian_z_grid)]).T


@dataclass
class Galaxy(Nebula):
    """
    Class defining the galaxy object (set up it as Sersic2D profile assuming it has continuum and emission components)
    """
    PA: u.degree = 90 * u.degree  # position angle of the major axis
    ax_ratio: float = 0.7  # ratio of minor/major axes
    r_eff: u.kpc = 1 * u.kpc  # Effective radius in kpc
    rad_lim: float = 3.  # Maximum radius for calculations (in R_eff)
    n: float = 1.  # Sersic index
    vel_rot: velunit = 0 * velunit  # Rotational velocity (not implemented yet)

    def __post_init__(self):
        self._assign_all_units()
        self._npix_los = 1
        self.r_max = self.r_eff.to(u.pc).value / self.pxscale.to(u.pc).value * self.rad_lim
        self._assign_position_params(conversion_type='galaxy')
        self._ref_line_id = 0
        self.linerat_constant = True  # True if the ratio of line fluxes shouldn't change across the nebula

    @cached_property
    def _brightness_3d_cartesian(self):
        """
        Method to obtain the brightness (or density) distribution of the nebula in cartesian coordinates
        """
        xx, yy = np.meshgrid(np.arange(self.pix_width), np.arange(self.pix_height))
        angle = (self.PA + 90 * u.degree).to(u.radian).value
        mod = Sersic2D(amplitude=1, r_eff=(self.r_eff.to(u.pc) / self.pxscale.to(u.pc)).value,
                       n=self.n, x_0=(self.pix_width - 1) / 2, y_0=(self.pix_height - 1) / 2,
                       ellip=1 - self.ax_ratio, theta=angle)
        brt = mod(xx, yy)

        xct = (xx - (self.pix_width - 1) / 2) * np.cos(angle) + \
              (yy - (self.pix_height - 1) / 2) * np.sin(angle)
        yct = (xx - (self.pix_width - 1) / 2) * np.sin(angle) - \
              (yy - (self.pix_height - 1) / 2) * np.cos(angle)
        rmaj = self.rad_lim * (self.r_eff.to(u.pc) / self.pxscale.to(u.pc)).value
        rmin = self.rad_lim * (self.r_eff.to(u.pc) / self.pxscale.to(u.pc)).value * self.ax_ratio
        mask = np.ones_like(brt, dtype=np.float32)
        rec = (xct ** 2 / rmaj ** 2) + (yct ** 2 / rmin ** 2) >= 1
        mask[rec] = 0
        mask = convolve_fft(mask, kernels.Gaussian2DKernel(3.), fill_value=0, allow_huge=True)
        brt = brt * mask
        brt = brt.reshape(self.pix_height, self.pix_width, 1)
        return brt

    def _get_2d_velocity(self):
        if hasattr(self, 'vel_rot') and (self.vel_rot is not None) and (self.vel_rot != 0):
            xx, yy = np.meshgrid(np.arange(self.pix_width), np.arange(self.pix_height))
            angle = (self.PA + 90 * u.degree).to(u.radian).value
            xct = (xx - (self.pix_width - 1) / 2) * np.cos(angle) + \
                  (yy - (self.pix_height - 1) / 2) * np.sin(angle)
            yct = (xx - (self.pix_width - 1) / 2) * np.sin(angle) - \
                  (yy - (self.pix_height - 1) / 2) * np.cos(angle)
            rad = np.sqrt(xct ** 2 + yct ** 2)
            vel_field = np.zeros_like(xx, dtype=np.float32) * velunit
            rec = rad > 0
            vel_field[rec] = self.vel_rot * np.sqrt(1 - self.ax_ratio ** 2) * xct[rec] / rad[rec]
            return vel_field
        else:
            return None



@dataclass
class DIG(Nebula):
    """
    Class defining the DIG component. For now it is defined just by its brightness (constant)
    """
    max_brightness: fluxunit = 1e-17 * fluxunit
    vel_gradient: (velunit / u.pc) = 0


@dataclass
class Cloud(Nebula):
    """Class of an isotropic spherical gas cloud without any ionization source.
    Defined by its position, radius, density, maximal optical depth"""
    radius: u.pc = 1.0 * u.pc
    max_brightness: fluxunit = 0 * fluxunit
    max_extinction: u.mag = 2.0 * u.mag
    thickness: float = 1.0
    perturb_degree: int = 0  # Degree of perturbations (max. degree of spherical harmonics for cloud)
    linerat_constant: bool = False  # True if the ratio of line fluxes shouldn't change across the nebula
    _phi_bins: int = 90
    _theta_bins: int = 90
    _rad_bins: int = 0
    _npix_los: int = 100

    def __post_init__(self):
        self._assign_all_units()
        if self._rad_bins == 0:
            self._rad_bins = np.ceil(self.radius.to(u.pc).value / self.pxscale.to(u.pc).value * 3).astype(int)
        delta = np.round((len(self._cartesian_y_grid) - 1) / 2).astype(int)
        if (self.xc is not None) and (self.yc is not None):
            self.x0 = self.xc - delta
            self.y0 = self.yc - delta
        elif (self.x0 is not None) and (self.y0 is not None):
            self.xc = self.x0 + delta
            self.yc = self.y0 + delta
        self._ref_line_id = 0

    @cached_property
    def _theta_grid(self):
        return np.linspace(0, np.pi, self._theta_bins)

    @cached_property
    def _phi_grid(self):
        return np.linspace(0, 2 * np.pi, self._phi_bins)

    @cached_property
    def _rad_grid(self):
        return np.linspace(0, self.radius, self._rad_bins)

    @cached_property
    def _cartesian_z_grid(self):
        npix = np.ceil(1.02 * self.radius / self.pxscale).astype(int)
        return np.linspace(-npix, npix, 2 * npix + 1) * self.pxscale

    @cached_property
    def _cartesian_y_grid(self):
        return self._cartesian_z_grid.copy()

    @cached_property
    def _cartesian_x_grid(self):
        return np.linspace(-1.02, 1.02, self._npix_los) * self.radius

    @cached_property
    def _brightness_3d_spherical(self):
        """
        Method to calculate brightness (or opacity) of the cloud at given theta, phi and radii

        theta: float -- polar angle [0, np.pi]
        phi: float -- azimuthal angle [0, 2 * np.pi]
        rad: float -- radius [0, self.radius]
        Returns:
            3D cube of normalized brightness in theta-phi-rad grid; total brightness = 1
        """
        rho, theta, phi = np.meshgrid(self._rad_grid, self._theta_grid, self._phi_grid, indexing='ij')
        brt = np.ones_like(theta)
        brt[rho < (self.radius * (1 - self.thickness))] = 0
        brt[rho > self.radius] = 0
        med = np.median(brt[brt > 0])
        if self.perturb_degree > 0:
            phi_cur = limit_angle(phi + np.random.uniform(0, 2 * np.pi, 1), 0, 2 * np.pi)
            theta_cur = limit_angle(theta + np.random.uniform(0, np.pi, 1), 0, np.pi)
            harm_amplitudes = self.perturb_amplitude * np.random.randn(self.perturb_degree * (self.perturb_degree + 2))

            brt += np.nansum(Parallel(n_jobs=lvmdatasimulator.n_process)(delayed(brightness_inhomogeneities_sphere)
                                                                         (harm_amplitudes, ll, phi_cur, theta_cur,
                                                                         rho, med, self.radius, self.thickness)
                                                                         for ll in np.arange(1,
                                                                                             self.perturb_degree + 1)),
                             axis=0)
            brt[brt < 0] = 0
        if med > 0:
            brt = brt / np.nansum(brt)
        return brt

    @cached_property
    def _brightness_4d_spherical(self):
        """
        Method to calculate brightness of the cloud at given theta, phi and radii for each line

        theta: float -- polar angle [0, np.pi]
        phi: float -- azimuthal angle [0, 2 * np.pi]
        rad: float -- radius [0, self.radius]
        Returns:
            4D cube of brightness in line-theta-phi-rad grid; normalized to the total brightness in Halpha
        """
        s = self._brightness_3d_spherical.shape
        if self.spectrum_id is None or self.linerat_constant:
            return self._brightness_3d_spherical.reshape((1, s[0], s[1], s[2]))
        rho, _, _ = np.meshgrid(self._rad_grid, self._theta_grid, self._phi_grid, indexing='ij')
        with fits.open(lvmdatasimulator.CLOUDY_MODELS) as hdu:
            radius = hdu[self.spectrum_id].data[0, 2:] * (self.thickness * self.radius) + \
                     self.radius * (1 - self.thickness)
            fluxes = hdu[self.spectrum_id].data[1:, 2:]
            radius = np.insert(radius, 0, self.radius * (1 - self.thickness))
            fluxes = np.insert(fluxes, 0, fluxes[:, 0], axis=1)
            index_ha = np.flatnonzero(hdu[self.spectrum_id].data[1:, 0] == 6562.81)
            if self.n_brightest_lines is not None and \
                    (self.n_brightest_lines > 0) and (self.n_brightest_lines < len(fluxes)):
                indexes_sorted = np.argsort(hdu[self.spectrum_id].data[1:, 1])[::-1]
                fluxes = fluxes[indexes_sorted[:self.n_brightest_lines], :]
                index_ha = np.flatnonzero(hdu[self.spectrum_id].data[1:, 0][indexes_sorted] == 6562.81)
            if len(index_ha) == 1:
                self._ref_line_id = index_ha[0]

            brt = np.array(Parallel(n_jobs=lvmdatasimulator.n_process)(delayed(sphere_brt_in_line)
                                                                       (self._brightness_3d_spherical, rho,
                                                                        radius, flux)
                                                                       for flux in fluxes)).reshape((fluxes.shape[0],
                                                                                                     s[0], s[1], s[2]))
            return brt / np.nansum(brt[self._ref_line_id])

    @cached_property
    def _brightness_3d_cartesian(self):
        return interpolate_sphere_to_cartesian(self._brightness_3d_spherical, x_grid=self._cartesian_x_grid,
                                               y_grid=self._cartesian_y_grid, z_grid=self._cartesian_z_grid,
                                               rad_grid=self._rad_grid, theta_grid=self._theta_grid,
                                               phi_grid=self._phi_grid, pxscale=self.pxscale)

    @cached_property
    def _brightness_4d_cartesian(self):
        s = self._brightness_4d_spherical.shape
        return np.array(Parallel(n_jobs=lvmdatasimulator.n_process)(delayed(interpolate_sphere_to_cartesian)
                                                                    (cur_line_array,
                                                                     self._cartesian_x_grid, self._cartesian_y_grid,
                                                                     self._cartesian_z_grid, self._rad_grid,
                                                                     self._theta_grid, self._phi_grid, self.pxscale)
                                                                    for cur_line_array in self._brightness_4d_spherical)
                        ).reshape((s[0], len(self._cartesian_z_grid), len(self._cartesian_y_grid),
                                   len(self._cartesian_x_grid)))


@dataclass
class Bubble(Cloud):
    """Class of an isotropic thin expanding bubble."""
    spectral_axis: velunit = np.arange(-20, 20, 10) * velunit
    expansion_velocity: velunit = 20 * velunit
    max_brightness: fluxunit = 1e-15 * fluxunit
    max_extinction: u.mag = 0 * u.mag
    thickness: float = 0.2

    @cached_property
    def _velocity_3d_spherical(self) -> velunit:
        """
        Calculate line of sight velocity at given radius, phi, theta

        V ~ 1/brightness (given that v~1/n_e^2 and brightness~ne^2)
        """
        rho, theta, phi = np.meshgrid(self._rad_grid, self._theta_grid, self._phi_grid, indexing='ij')
        vel_cube = np.zeros_like(self._brightness_3d_spherical)
        rec = (rho <= self.radius) & (rho >= (self.radius * (1 - self.thickness)))
        vel_cube[rec] = \
            np.sin(theta[rec]) * \
            np.cos(phi[rec]) * \
            self.expansion_velocity / self._brightness_3d_spherical[rec] * \
            np.median(self._brightness_3d_spherical[self._brightness_3d_spherical > 0])
        return vel_cube

    @cached_property
    def _velocity_3d_cartesian(self) -> velunit:
        return interpolate_sphere_to_cartesian(self._velocity_3d_spherical, x_grid=self._cartesian_x_grid,
                                               y_grid=self._cartesian_y_grid, z_grid=self._cartesian_z_grid,
                                               rad_grid=self._rad_grid, theta_grid=self._theta_grid,
                                               phi_grid=self._phi_grid, pxscale=self.pxscale)

    def _turbulent_lsf(self, velocity):
        """Line spread function as a function of coorinates, including the velocity center shift"""
        # mu = self.velocity(theta, phi)
        mu = self._velocity_3d_cartesian[:, :, :, None] * velunit + self.sys_velocity
        sig = self.turbulent_sigma
        return 1. / (np.sqrt(2. * np.pi) * sig) * np.exp(-np.power((velocity - mu) / sig, 2.) / 2)

    def _d_spectrum_cartesian(self, velocity: velunit):
        """Returns local spectrum, per pc**3 of area"""
        return (self._brightness_3d_cartesian[:, :, :, None] * (
                fluxunit / u.pc ** 3) * self._turbulent_lsf(velocity)).to(fluxunit / velunit / u.pc ** 3)

    @cached_property
    def line_profile(self) -> (fluxunit / velunit):
        """
        Produces the distribution of the observed line profiles in each pixels of the sky plane
        """
        vel_axis = self.spectral_axis.to(velunit, equivalencies=u.spectral())
        _, _, _, vels = np.meshgrid(self._cartesian_z_grid,
                                    self._cartesian_y_grid,
                                    self._cartesian_x_grid,
                                    vel_axis, indexing='ij')
        spectrum = (
                np.sum(self._d_spectrum_cartesian(vels), axis=2
                       ).T * (self._cartesian_x_grid[1] - self._cartesian_x_grid[0]
                              ) * (self._cartesian_y_grid[1] - self._cartesian_y_grid[0]
                                   ) * (self._cartesian_z_grid[1] - self._cartesian_z_grid[0])
        )
        return spectrum / np.sum(spectrum, axis=0)


@dataclass
class CustomNebula(Nebula):
    """
        Class defining the custom nebulae with the user-defined distribution of the brighntess, continuum
        and line shapes in different lines.
    """
    pass


@dataclass
class ISM:
    """
    Class defining the ISM contribution to the field of view
    """
    wcs: WCS
    width: int = 401  # Width of field of view in pixels
    height: int = 401  # Width of field of view in pixels
    spec_resolution: u.Angstrom = 0.06 * u.Angstrom  # Spectral resolution of the simulation
    npix_line: int = 1  # Minimal number of pixels for a resolution element at wl = 10000A for construction of vel.grid
    distance: u.kpc = 50 * u.kpc  # Distance to the object for further conversion between arcsec and pc
    sys_velocity: velunit = 0 * velunit  # Systemic velocity to center the vel.grid on
    vel_amplitude: velunit = 150 * velunit  # Maximal deviation from the systemic velocity to setup vel.grid
    turbulent_sigma: velunit = 20. * velunit  # turbulence vel. disp. to be used for every nebula unless other specified
    R_V: float = 3.1  # R_V value defining the reddening curve (to be used unless other value is provided for a nebula)
    ext_law: str = 'F99'  # Reddening law (to be used unless other value is provided for a nebula)

    def __post_init__(self):
        assign_units(self, ['vel_amplitude', 'turbulent_sigma', 'sys_velocity', 'distance', 'spec_resolution'],
                     [velunit, velunit, velunit, u.kpc, u.Angstrom])
        self.content = fits.HDUList()
        self.content.append(fits.PrimaryHDU(header=self.wcs.to_header(), data=np.zeros(shape=(2, 2), dtype=int)))
        self.vel_grid = np.linspace(-self.vel_amplitude + self.sys_velocity,
                                    self.vel_amplitude + self.sys_velocity,
                                    np.ceil(self.vel_amplitude / self.vel_resolution).astype(int) * 2 + 1)
        self.pxscale = proj_plane_pixel_scales(self.wcs)[0] * 3600 * self.distance.to(u.pc) / 206265.
        self.content[0].header['width'] = (self.width, "Width of field of view, px")
        self.content[0].header['height'] = (self.height, "Height of field of view, px")
        self.content[0].header['PhysRes'] = (self.pxscale.value, "Physical resolution, pc/px")
        self.content[0].header['Dist'] = (self.distance.value, "Distance, kpc")
        self.content[0].header['Vsys'] = (self.distance.value, "Systemic Velocity, km/s")
        self.content[0].header['VelRes'] = (self.vel_resolution.value, "Velocity resolution, km/s/px")
        self.content[0].header['TurbSig'] = (self.turbulent_sigma.value, "Default turbulent velocity dispersion, km/s")
        self.content[0].header['Nobj'] = (0, "Number of generated nebulae")

    @cached_property
    def vel_resolution(self):
        return (self.spec_resolution / self.npix_line / (10000 * u.Angstrom) * c.c).to(velunit)

    def _get_continuum(self, my_comp, wl_grid):
        """
        Properly extracts continuum for current nebula taking into account its shape and surface brightness
        :param my_comp:
        :param wl_grid:
        :return: continuum
        """
        cont_type = self.content[my_comp + "_CONTINUUM"].header.get("CONTTYPE")
        continuum = self.content[my_comp + "_CONTINUUM"].data
        cont_norm = self.content[my_comp + "_CONTINUUM"].header.get("CONTFLUX")
        cont_norm_wl = self.content[my_comp + "_CONTINUUM"].header.get("CONTWL")
        if cont_type.lower() == 'model':
            cont_wl_fullrange = continuum[0, :]
            cont_fullrange = continuum[1, :]
            p = interp1d(continuum[0, :], continuum[1, :], assume_sorted=True,
                         bounds_error=False, fill_value='extrapolate')
            continuum = p(wl_grid)
        elif cont_type.lower() == 'poly':
            p = np.poly1d(continuum)
            cont_wl_fullrange = np.linspace(3500, 10000, 501)
            cont_fullrange = p(cont_wl_fullrange)
            continuum = p(wl_grid)
        elif cont_type.lower() == 'bb':
            cont_wl_fullrange = np.linspace(3500, 10000, 501)
            cont_fullrange = 1 / cont_wl_fullrange ** 5 / (
                    np.exp(6.63e-27 * 3e10 / cont_wl_fullrange / 1e-8 / continuum / 1.38e-16) - 1)
            continuum = 1 / wl_grid ** 5 / (np.exp(6.63e-27 * 3e10 / wl_grid / 1e-8 / continuum / 1.38e-16) - 1)

        t_filter = None
        if type(cont_norm_wl) == str:
            file_filter = os.path.join(lvmdatasimulator.ROOT_DIR, 'data', 'instrument', 'filters', cont_norm_wl+".dat")
            if not os.path.isfile(file_filter):
                log.warning("Cannot find filter {0}. "
                            "Default Wavelength = 5500A will be used for continuum normalization".format(cont_norm_wl))
                cont_norm_wl = 5500.
                t_filter = None
            else:
                t_filter = ascii.read(file_filter, names=['lambda', 'transmission'])
                cont_norm_wl = np.sum(t_filter['lambda'] * t_filter['transmission']) / np.sum(t_filter['transmission'])
        if t_filter is None:
            cont_model_max = cont_fullrange[np.argmin(abs(cont_wl_fullrange - cont_norm_wl))]
        else:
            dl = np.roll(t_filter['lambda'], -1) - t_filter['lambda']
            dl[-1] = dl[-2]
            w_filter = np.sum(dl * t_filter['transmission'])/np.max(t_filter['transmission'])
            p = interp1d(t_filter['lambda'], t_filter['transmission'], assume_sorted=True,
                         fill_value=0, bounds_error=False)
            cont_model_max = np.sum(cont_fullrange * p(cont_wl_fullrange)) / w_filter
        if ~np.isfinite(cont_norm) or cont_norm <= 0:
            cont_norm = self.content[my_comp + "_CONTINUUM"].header.get("CONTMAG") * u.ABmag
            cont_norm = cont_norm.to(u.STmag, u.spectral_density(cont_norm_wl * u.AA)).to(u.erg/u.s/u.cm**2/u.AA).value
        return continuum / cont_model_max * cont_norm * (wl_grid[1] - wl_grid[0])

    def _add_fits_extension(self, name, value, obj_to_add, zorder=0, cur_wavelength=0, add_fits_kw=None,
                            add_counter=False):
        self.content.append(fits.ImageHDU(np.atleast_1d(value), name=name))
        self.content[-1].header['Nebtype'] = (type(obj_to_add).__name__, "Type of the nebula")
        is_dark = ((obj_to_add.max_brightness <= 0) and (obj_to_add.max_extinction > 0))
        self.content[-1].header['Dark'] = (is_dark, " Emitting or absorbing nebula?")
        self.content[-1].header['X0'] = (obj_to_add.x0, "Position in the field of view")
        self.content[-1].header['Y0'] = (obj_to_add.y0, "Position in the field of view")
        self.content[-1].header['Zorder'] = (zorder, "Z-order in the field of view")
        self.content[-1].header['NCHUNKS'] = (obj_to_add.nchunks, "N of chunks for convolution")
        if type(obj_to_add) in [Bubble, Cloud]:
            self.content[-1].header['Radius'] = (obj_to_add.radius.to_value(u.pc), "Radius of the nebula, pc")
            self.content[-1].header['PertOrd'] = (obj_to_add.perturb_degree, "Degree to produce random perturbations")
            self.content[-1].header['PertAmp'] = (obj_to_add.perturb_amplitude, "Max amplitude of random perturb.")
        if type(obj_to_add) in [DIG, Nebula, Rectangle, Ellipse, Circle]:
            self.content[-1].header['PertScl'] = (obj_to_add.perturb_scale.to_value(u.pc),
                                                  "Scale of the random perturbations, pc")
            self.content[-1].header['PertAmp'] = (obj_to_add.perturb_amplitude, "Max amplitude of random perturb.")
        if type(obj_to_add) in [Filament]:
            self.content[-1].header['Width'] = (obj_to_add.width.to_value(u.pc), 'Width of the filament, pc')
            self.content[-1].header['Length'] = (obj_to_add.length.to_value(u.pc), 'Length of the filament, pc')

        if type(obj_to_add) not in [Bubble, Cloud, Galaxy]:
            if obj_to_add.vel_gradient is not None:
                self.content[-1].header['Vgrad'] = (obj_to_add.vel_gradient.to_value(velunit / u.pc),
                                                    'Velocity gradient, km/s per pc')
                if obj_to_add.vel_pa is not None:
                    self.content[-1].header['PAkin'] = (obj_to_add.vel_pa.to_value(u.degree),
                                                        'Kinematical PA, degree')
        if type(obj_to_add) in [Nebula, Rectangle]:
            self.content[-1].header['Width'] = (obj_to_add.width.to_value(u.pc), 'Width of the nebula, pc')
            self.content[-1].header['Height'] = (obj_to_add.height.to_value(u.pc), 'Height of the nebula, pc')
        if type(obj_to_add) in [Ellipse, Filament, Galaxy]:
            self.content[-1].header['PA'] = (obj_to_add.PA.to_value(u.degree), 'Position angle, degree')
        if type(obj_to_add) in [Galaxy]:
            self.content[-1].header['Reff'] = (obj_to_add.r_eff.to_value(u.kpc), 'Effective radius of the galaxy, kpc')
            self.content[-1].header['NSersic'] = (obj_to_add.n, 'Sersic index')
            self.content[-1].header['Rlim'] = (obj_to_add.rad_lim, 'Limiting distance in Reff')
            self.content[-1].header['Vrot'] = (obj_to_add.vel_rot.to_value(velunit),
                                               'Rotational velocity, km/s')
            self.content[-1].header['PAkin'] = (obj_to_add.vel_pa.to_value(u.degree),
                                                'Kinematical PA, degree')
        if type(obj_to_add) in [Ellipse, Circle]:
            self.content[-1].header['Radius'] = (obj_to_add.radius.to_value(u.kpc), 'Radius (major axis), pc')
        if type(obj_to_add) in [Ellipse, Galaxy]:
            self.content[-1].header['AxRat'] = (obj_to_add.ax_ratio, "Axis ratio")
        if (obj_to_add.max_brightness <= 0) and (obj_to_add.max_extinction > 0):
            self.content[-1].header['MaxExt'] = (obj_to_add.max_extinction.value, "Max extinction, mag/pix")
        elif obj_to_add.max_brightness > 0:
            self.content[-1].header['MaxBrt'] = (obj_to_add.max_brightness.value, "Max brightness, erg/s/cm^2/arcsec^2")
        if type(obj_to_add) == Bubble:
            self.content[-1].header['Vexp'] = (obj_to_add.expansion_velocity.to_value(velunit),
                                               'Expansion velocity, km/s')
        self.content[-1].header['SysVel'] = (obj_to_add.sys_velocity.to_value(velunit), "Systemic velocity, km/s")
        self.content[-1].header['TurbVel'] = (obj_to_add.turbulent_sigma.to_value(velunit),
                                              "ISM Velocity dispersion, km/s")
        self.content[-1].header['SpecID'] = (obj_to_add.spectrum_id, "Ref. spectrum ID in model grid")
        self.content[-1].header['NLines'] = (obj_to_add.n_brightest_lines, "Maximal number of lines to use")
        if cur_wavelength:
            self.content[-1].header['Lambda'] = (cur_wavelength, "Current line wavelength")
        if add_fits_kw is not None:
            for kw in add_fits_kw:
                self.content[-1].header[kw] = add_fits_kw[kw]
        if add_counter:
            self.content[0].header['Nobj'] = (self.content[0].header['Nobj'] + 1, "Total number of nebulae")

    def add_nebula(self, obj_to_add, obj_id=0, zorder=0, add_fits_kw=None, continuum=None):
        """
        Method to add the particular nebula to the ISM object and to the output multi-extensions fits file
        """
        if type(obj_to_add) not in [Nebula, Bubble, Filament, DIG, Cloud, Galaxy, Ellipse, Circle,
                                    Rectangle, CustomNebula]:
            log.warning('Skip nebula of wrong type ({0})'.format(type(obj_to_add)))
            return

        if (obj_to_add.max_brightness <= 0) and (obj_to_add.max_extinction <= 0) and (continuum is None):
            log.warning('Skip nebula with zero extinction and brightness')
            return
        if obj_to_add.max_brightness > 0:
            brt = obj_to_add.brightness_skyplane.value
            if obj_to_add.spectrum_id is not None and not obj_to_add.linerat_constant:
                brt_4d = obj_to_add.brightness_skyplane_lines.value
            else:
                brt_4d = None
        elif obj_to_add.max_extinction > 0:
            brt = obj_to_add.extinction_skyplane.value
            brt_4d = None
        elif continuum is not None:
            brt = obj_to_add.brightness_skyplane
            brt_4d = None

        self._add_fits_extension(name="Comp_{0}_Brightness".format(obj_id), value=brt,
                                 obj_to_add=obj_to_add, zorder=zorder, add_fits_kw=add_fits_kw, add_counter=True)
        if obj_to_add.max_brightness > 0:
            if brt_4d is not None:
                with fits.open(lvmdatasimulator.CLOUDY_MODELS) as hdu:
                    wl_list = hdu[obj_to_add.spectrum_id].data[1:, 0]
                    if obj_to_add.n_brightest_lines is not None and \
                            (obj_to_add.n_brightest_lines > 0) and (obj_to_add.n_brightest_lines < len(wl_list)):
                        wl_list = wl_list[np.argsort(hdu[obj_to_add.spectrum_id].data[1:, 1]
                                                     )[::-1][: obj_to_add.n_brightest_lines]]
                for line_ind in range(brt_4d.shape[0]):
                    self._add_fits_extension(name="Comp_{0}_Flux_{1}".format(obj_id, wl_list[line_ind]),
                                             value=brt_4d[line_ind],
                                             obj_to_add=obj_to_add, zorder=zorder, add_fits_kw=add_fits_kw,
                                             cur_wavelength=wl_list[line_ind])
            elif obj_to_add.spectrum_id is not None and obj_to_add.linerat_constant:
                with fits.open(lvmdatasimulator.CLOUDY_MODELS) as hdu:
                    data_save = hdu[obj_to_add.spectrum_id].data[1:, :2]
                    if obj_to_add.n_brightest_lines is not None and \
                            (obj_to_add.n_brightest_lines > 0) and \
                            (obj_to_add.n_brightest_lines < len(hdu[obj_to_add.spectrum_id].data[1:, 0])):
                            data_save = data_save[np.argsort(data_save[:, 1])[::-1][: obj_to_add.n_brightest_lines]]
                    self._add_fits_extension(name="Comp_{0}_FluxRatios".format(obj_id),
                                             value=data_save.T,
                                             obj_to_add=obj_to_add, zorder=zorder, add_fits_kw=add_fits_kw)

        if type(obj_to_add) == Bubble:
            self._add_fits_extension(name="Comp_{0}_LineProfile".format(obj_id), value=obj_to_add.line_profile.value,
                                     obj_to_add=obj_to_add, zorder=zorder, add_fits_kw=add_fits_kw)

        if obj_to_add.vel_field is not None:
            self._add_fits_extension(name="Comp_{0}_Vel".format(obj_id), value=obj_to_add.vel_field.value,
                                     obj_to_add=obj_to_add, zorder=zorder, add_fits_kw=add_fits_kw)

        if continuum is not None:
            self._add_fits_extension(name="Comp_{0}_Continuum".format(obj_id), value=continuum,
                                     obj_to_add=obj_to_add, zorder=zorder, add_fits_kw=add_fits_kw)
        return self.content

    def save_ism(self, filename):
        self.content.writeto(filename, overwrite=True)
        log.info("Generated ISM saved to {0}".format(filename))

    def generate(self, all_objects):
        """
        Generate all the Nebulae from the input list

        Args:
            all_objects: list -- contains a dictionary describing the nebula to add:
                example:
                all_objects = [{type: "Bubble, Filament, DIG, ....",
                                sys_velocity: 0 * u.km/u.s,
                                expansion_velocity: 30 * u.km/u.s,
                                turbulent_sigma: 10 * u.km/u.s,
                                radius: 5 * u.pc,
                                max_brightness: 1e-16 * u.erg / u.cm**2 / u.s / u.arcsec ** 2,
                                RA: "08h12m13s",
                                DEC: "-20d14m13s",
                                'perturb_degree': 8, # max. order of spherical harmonics to generate inhomogeneities
                                'perturb_amplitude': 0.1, # relative max. amplitude of inhomogeneities,
                                'perturb_scale': 200 * u.pc,  # spatial scale to generate inhomogeneities (DIG only)
                                'distance': 50 * u.kpc,  # distance to the nebula (default is from ISM)
                                'cloudy_id': None,  # id of pre-generated Cloudy model
                                'cloudy_params': {'Z': 0.5, 'qH': 49., 'nH': 10, 'Teff': 30000, 'Geometry': 'Sphere'}, #
                                    parameters defining the pre-generated Cloudy model (used if spectrum_id is None)
                                'n_brightest_lines': 10,  # only this number of the brightest lines will be evaluated
                                'linerat_constant': False #  if True -> lines ratios don't vary across Cloud/Bubble
                                'continuum_type': 'BB' or 'Model' or 'Poly' # type of the continuum model
                                'continuum_data': model_id or [poly_coefficients] or Teff or dict with "wl" and "flux"
                                                # value defining cont. shape
                                'continuum_flux': 1e-16 * u.erg / u.cm ** 2 / u.s / u.arcsec **2 / u.AA,
                                'continuum_mag': 22 * u.mag,
                                'continuum_wl': 5500, # could be also R, V, B,
                                'ext_law': 'F99',  # Extinction law, one of those used by pyneb (used for dark nebulae)
                                'ext_rv': 3.1,  # Value of R_V for extinction curve calculation (used for dark nebulae)
                                'vel_gradient: 12. * u.km / u.s / u.pc  # Line-of-sight velocity gradient
                                'vel_pa: 30. * u.degree  # PA of kinematical axis (for vel_gradient or vel_rot)
                                }]
        """
        if type(all_objects) is dict:
            all_objects = [all_objects]
        if type(all_objects) not in [list, tuple]:
            log.warning('Cannot generate nebulae as the input is not a list or tuple')
            return None
        all_objects = [cobj for cobj in all_objects if cobj.get('type') in ['Nebula', 'Bubble', 'Galaxy',
                                                                            'Filament', 'DIG', 'Cloud',
                                                                            'Rectangle', 'Circle', 'Ellipse',
                                                                            'CustomNebula']]
        n_objects = len(all_objects)
        log.info("Start generating {} nebulae".format(n_objects))
        bar = progressbar.ProgressBar(max_value=n_objects).start()
        obj_id = self.content[0].header['Nobj']
        obj_id_ini = self.content[0].header['Nobj']
        for ind_obj, cur_obj in enumerate(all_objects):
            bar.update(ind_obj)
            # Setup default parameters for missing keywords
            kin_pa_default = 0
            if 'PA' in cur_obj:
                kin_pa_default = cur_obj['PA']
            for k, v, unit in zip(['max_brightness', 'max_extinction', 'thickness',
                                   'expansion_velocity', 'sys_velocity',
                                   'turbulent_sigma', 'perturb_degree',
                                   'perturb_amplitude', 'perturb_scale', 'radius', 'distance',
                                   'continuum_type', 'continuum_data', 'continuum_flux', 'continuum_mag',
                                   'continuum_wl', 'ext_law', 'ext_rv', 'vel_gradient', 'vel_rot', 'vel_pa',
                                   'n_brightest_lines', 'offset_RA', 'offset_DEC', 'RA', 'DEC'],
                                  [0, 0, 1., 0, self.sys_velocity, self.turbulent_sigma, 0, 0.1, 0, 0, self.distance,
                                   None, None, 0, None, 5500., self.ext_law, self.R_V, 0, 0, kin_pa_default, None,
                                   None, None, None, None],
                                  [fluxunit, u.mag, None, velunit, velunit, velunit, None, None,
                                   u.pc, u.pc, u.kpc, None, None, fluxunit/u.AA, None, u.Angstrom,
                                   None, None, velunit / u.pc, velunit, u.degree, None, u.arcsec, u.arcsec,
                                   u.degree, u.degree]):
                set_default_dict_values(cur_obj, k, v, unit=unit)

            for k in ['max_brightness', 'max_extinction', 'radius', 'continuum_flux']:
                if cur_obj[k] < 0:
                    cur_obj[k] = 0

            if (cur_obj['max_brightness'] == 0) and (cur_obj['max_extinction'] == 0) and \
                    (((cur_obj['continuum_mag'] is None) and (cur_obj['continuum_flux'] == 0)) or
                     (cur_obj['continuum_data'] is None) or (cur_obj['continuum_type'] is None)):
                log.warning("Wrong set of parameters for the nebula #{0}: skip this one".format(ind_obj))
                continue

            if lvmdatasimulator.CLOUDY_MODELS is None or (cur_obj['max_brightness'] <= 0):
                cloudy_model_index = None
                cloudy_model_id = None
            else:
                if cur_obj.get('cloudy_id') is None:
                    if cur_obj.get('cloudy_params') is None or (type(cur_obj.get('cloudy_params')) is not dict):
                        log.warning("Neither of 'cloudy_id' or 'cloudy_params' is set for the nebula #{0}: "
                                    "use default 'cloudy_id={1}'".format(ind_obj,
                                                                         lvmdatasimulator.CLOUDY_SPEC_DEFAULTS['id']))
                        cur_obj['cloudy_id'] = lvmdatasimulator.CLOUDY_SPEC_DEFAULTS['id']
                    else:
                        for p in lvmdatasimulator.CLOUDY_SPEC_DEFAULTS:
                            if p == 'id':
                                continue
                            if cur_obj['cloudy_params'].get(p) is None:
                                cur_obj['cloudy_params'][p] = lvmdatasimulator.CLOUDY_SPEC_DEFAULTS[p]

                cloudy_model_index, cloudy_model_id = find_model_id(file=lvmdatasimulator.CLOUDY_MODELS,
                                                                    check_id=cur_obj.get('cloudy_id'),
                                                                    params=cur_obj.get('cloudy_params'))

            if cur_obj.get('linerat_constant') is None:
                if cur_obj['type'] in ['Bubble', 'Cloud']:
                    cur_obj['linerat_constant'] = False
                else:
                    cur_obj['linerat_constant'] = True

            if cur_obj['type'] == 'DIG':
                if cur_obj.get('max_brightness') is None or cur_obj.get('max_brightness') <= 0:
                    log.warning("Wrong set of parameters for the nebula #{0}: skip this one".format(ind_obj))
                    continue
                if not cur_obj.get('zorder'):
                    cur_obj['zorder'] = -1
                if cur_obj['perturb_scale'] < 0:
                    cur_obj['perturb_scale'] = 0
                generated_object = DIG(max_brightness=cur_obj.get('max_brightness'),
                                       turbulent_sigma=cur_obj['turbulent_sigma'],
                                       sys_velocity=cur_obj['sys_velocity'],
                                       vel_gradient=cur_obj['vel_gradient'],
                                       vel_pa=cur_obj['vel_pa'],
                                       spectrum_id=cloudy_model_index,
                                       n_brightest_lines=cur_obj['n_brightest_lines'],
                                       pxscale=self.pxscale * (cur_obj['distance'].to(u.pc) / self.distance.to(u.pc)),
                                       perturb_scale=cur_obj['perturb_scale'],
                                       perturb_amplitude=cur_obj['perturb_amplitude'],
                                       pix_width=self.width, pix_height=self.height,
                                       )
            else:
                # ==== Check input parameters and do necessary conversions
                if not cur_obj.get('zorder'):
                    cur_obj['zorder'] = 0
                if not ((cur_obj.get('RA') is not None and cur_obj.get('DEC') is not None) or
                        (cur_obj.get('X') is not None and cur_obj.get('Y') is not None) or
                        (cur_obj.get('offset_X') is not None and cur_obj.get('offset_Y') is not None) or
                        (cur_obj.get('offset_RA') is not None and cur_obj.get('offset_DEC') is not None)):
                    log.warning("Wrong set of parameters for the nebula #{0}: skip this one".format(ind_obj))
                    continue
                if cur_obj['type'] in ['Rectangle', 'Nebula'] and not (('width' in cur_obj) and ('height' in cur_obj)):
                    log.warning("Wrong set of parameters for the nebula #{0}: skip this one".format(ind_obj))
                    continue
                if (cur_obj['type'] in ["Bubble", "Cloud", "Ellipse", 'Circle']) and (cur_obj['radius'] == 0):
                    log.warning("Wrong set of parameters for the nebula #{0}: skip this one".format(ind_obj))
                    continue
                if cur_obj['type'] == 'Filament' and not (('length' in cur_obj) and ('PA' in cur_obj)):
                    log.warning("Wrong set of parameters for the nebula #{0}: skip this one".format(ind_obj))
                    continue
                if cur_obj['type'] == 'Galaxy' and not (('r_eff' in cur_obj) and ('PA' in cur_obj) and
                                                        ('ax_ratio' in cur_obj)):
                    log.warning("Wrong set of parameters for the nebula #{0}: skip this one".format(ind_obj))
                    continue
                if cur_obj['type'] == 'Galaxy':
                    if 'n' not in cur_obj:
                        log.info("Set default Sersic index n=1 for the nebula #{0}".format(ind_obj))
                        cur_obj['n'] = 1
                    if 'rad_lim' not in cur_obj:
                        cur_obj['rad_lim'] = 3.

                if not (cur_obj.get('X') is not None and cur_obj.get('Y') is not None):
                    if cur_obj.get('offset_X') is not None and cur_obj.get('offset_Y') is not None:
                        x = (self.width - 1) / 2. + cur_obj.get('offset_X')
                        y = (self.height - 1) / 2. + cur_obj.get('offset_Y')
                    elif cur_obj.get('RA') is not None and cur_obj.get('DEC') is not None:
                        radec = SkyCoord(ra=cur_obj.get('RA'), dec=cur_obj.get('DEC'))
                        x, y = self.wcs.world_to_pixel(radec)
                    elif cur_obj.get('offset_RA') is not None and cur_obj.get('offset_DEC') is not None:
                        x = (self.width - 1) / 2. - (cur_obj.get('offset_RA').to_value(u.degree) /
                                                     proj_plane_pixel_scales(self.wcs)[0])
                        y = (self.height - 1) / 2. + (cur_obj.get('offset_DEC').to_value(u.degree) /
                                                     proj_plane_pixel_scales(self.wcs)[0])
                    x = np.round(x).astype(int)
                    y = np.round(y).astype(int)
                else:
                    x, y = [cur_obj.get('X'), cur_obj.get('Y')]

                if (cur_obj['thickness'] <= 0) or (cur_obj['thickness'] > 1):
                    log.warning("Wrong value of thickness of the nebula #{0}: set it to 1.".format(ind_obj))
                    cur_obj['thickness'] = 1.

                if cur_obj['type'] == "Bubble" and cur_obj.get('expansion_velocity') <= 0:
                    log.warning("Contracting bubbles are not supported (nebula #{0})."
                                " Use non-expanding cloud instead".format(ind_obj))
                    cur_obj['type'] = "Cloud"

                if cur_obj['type'] in ["Bubble", "Cloud"]:
                    if cur_obj['perturb_degree'] < 0:
                        cur_obj['perturb_degree'] = 0

                elif cur_obj['type'] == 'Filament' and ('width' not in cur_obj):
                    log.info("Set default width of the filament 0.1 pc for the nebula #{0}".format(ind_obj))
                    cur_obj['width'] = 0.1 * u.pc

                # ==== Start calculations of different nebulae
                if cur_obj['type'] == "Bubble":
                    generated_object = Bubble(xc=x, yc=y,
                                              max_brightness=cur_obj.get('max_brightness'),
                                              max_extinction=cur_obj.get('max_extinction'),
                                              spectral_axis=self.vel_grid,
                                              expansion_velocity=cur_obj.get('expansion_velocity'),
                                              thickness=cur_obj['thickness'],
                                              radius=cur_obj['radius'],
                                              pxscale=self.pxscale * (cur_obj['distance'].to(u.pc) /
                                                                      self.distance.to(u.pc)),
                                              perturb_degree=cur_obj['perturb_degree'],
                                              perturb_amplitude=cur_obj['perturb_amplitude'],
                                              turbulent_sigma=cur_obj['turbulent_sigma'],
                                              sys_velocity=cur_obj['sys_velocity'],
                                              spectrum_id=cloudy_model_index,
                                              n_brightest_lines=cur_obj['n_brightest_lines'],
                                              linerat_constant=cur_obj['linerat_constant'],
                                              )
                elif cur_obj['type'] == "Cloud":
                    generated_object = Cloud(xc=x, yc=y,
                                             max_brightness=cur_obj.get('max_brightness'),
                                             max_extinction=cur_obj.get('max_extinction'),
                                             thickness=cur_obj['thickness'],
                                             radius=cur_obj['radius'],
                                             pxscale=self.pxscale * (cur_obj['distance'].to(u.pc) /
                                                                     self.distance.to(u.pc)),
                                             perturb_degree=cur_obj['perturb_degree'],
                                             perturb_amplitude=cur_obj['perturb_amplitude'],
                                             spectrum_id=cloudy_model_index,
                                             n_brightest_lines=cur_obj['n_brightest_lines'],
                                             turbulent_sigma=cur_obj['turbulent_sigma'],
                                             sys_velocity=cur_obj['sys_velocity'],
                                             linerat_constant=cur_obj['linerat_constant'],
                                             vel_gradient=cur_obj['vel_gradient'],
                                             vel_pa=cur_obj['vel_pa'],
                                             )

                elif cur_obj['type'] == "Filament":
                    generated_object = Filament(xc=x, yc=y,
                                                max_brightness=cur_obj.get('max_brightness'),
                                                max_extinction=cur_obj.get('max_extinction'),
                                                width=cur_obj['width'],
                                                length=cur_obj['length'],
                                                PA=cur_obj['PA'],
                                                vel_gradient=cur_obj['vel_gradient'],
                                                vel_pa=cur_obj['vel_pa'],
                                                spectrum_id=cloudy_model_index,
                                                n_brightest_lines=cur_obj['n_brightest_lines'],
                                                turbulent_sigma=cur_obj['turbulent_sigma'],
                                                sys_velocity=cur_obj['sys_velocity'],
                                                pxscale=self.pxscale * (cur_obj['distance'].to(u.pc) /
                                                                        self.distance.to(u.pc)),
                                                )

                elif cur_obj['type'] == "Galaxy":
                    generated_object = Galaxy(xc=x, yc=y,
                                              max_brightness=cur_obj.get('max_brightness'),
                                              max_extinction=cur_obj.get('max_extinction'),
                                              r_eff=cur_obj['r_eff'],
                                              rad_lim=cur_obj['rad_lim'],
                                              ax_ratio=cur_obj['ax_ratio'],
                                              PA=cur_obj['PA'],
                                              n=cur_obj['n'],
                                              vel_rot=cur_obj.get('vel_rot'),
                                              vel_pa=cur_obj['vel_pa'],
                                              spectrum_id=cloudy_model_index,
                                              n_brightest_lines=cur_obj['n_brightest_lines'],
                                              turbulent_sigma=cur_obj['turbulent_sigma'],
                                              sys_velocity=cur_obj['sys_velocity'],
                                              pxscale=self.pxscale * (cur_obj['distance'].to(u.pc) /
                                                                      self.distance.to(u.pc)),
                                              )
                elif cur_obj['type'] == "Ellipse":
                    generated_object = Ellipse(xc=x, yc=y,
                                               max_brightness=cur_obj.get('max_brightness'),
                                               max_extinction=cur_obj.get('max_extinction'),
                                               radius=cur_obj['radius'],
                                               ax_ratio=cur_obj['ax_ratio'],
                                               PA=cur_obj['PA'],
                                               spectrum_id=cloudy_model_index,
                                               n_brightest_lines=cur_obj['n_brightest_lines'],
                                               turbulent_sigma=cur_obj['turbulent_sigma'],
                                               sys_velocity=cur_obj['sys_velocity'],
                                               vel_gradient=cur_obj['vel_gradient'],
                                               vel_pa=cur_obj['vel_pa'],
                                               pxscale=self.pxscale * (cur_obj['distance'].to(u.pc) /
                                                                       self.distance.to(u.pc)),
                                               perturb_scale=cur_obj['perturb_scale'],
                                               perturb_amplitude=cur_obj['perturb_amplitude'],
                                               )
                elif cur_obj['type'] == "Circle":
                    generated_object = Circle(xc=x, yc=y,
                                              max_brightness=cur_obj.get('max_brightness'),
                                              max_extinction=cur_obj.get('max_extinction'),
                                              radius=cur_obj['radius'],
                                              spectrum_id=cloudy_model_index,
                                              n_brightest_lines=cur_obj['n_brightest_lines'],
                                              turbulent_sigma=cur_obj['turbulent_sigma'],
                                              sys_velocity=cur_obj['sys_velocity'],
                                              vel_gradient=cur_obj['vel_gradient'],
                                              vel_pa=cur_obj['vel_pa'],
                                              pxscale=self.pxscale * (cur_obj['distance'].to(u.pc) /
                                                                      self.distance.to(u.pc)),
                                              perturb_scale=cur_obj['perturb_scale'],
                                              perturb_amplitude=cur_obj['perturb_amplitude'],
                                              )
                elif cur_obj['type'] == "Rectangle" or (cur_obj['type'] == "Nebula"):
                    generated_object = Rectangle(xc=x, yc=y,
                                                 width=cur_obj.get('width'), height=cur_obj.get('height'),
                                                 max_brightness=cur_obj.get('max_brightness'),
                                                 max_extinction=cur_obj.get('max_extinction'),
                                                 spectrum_id=cloudy_model_index,
                                                 n_brightest_lines=cur_obj['n_brightest_lines'],
                                                 turbulent_sigma=cur_obj['turbulent_sigma'],
                                                 sys_velocity=cur_obj['sys_velocity'],
                                                 vel_gradient=cur_obj['vel_gradient'],
                                                 vel_pa=cur_obj['vel_pa'],
                                                 pxscale=self.pxscale * (cur_obj['distance'].to(u.pc) /
                                                                         self.distance.to(u.pc)),
                                                 perturb_scale=cur_obj['perturb_scale'],
                                                 perturb_amplitude=cur_obj['perturb_amplitude'],
                                                 )
                elif cur_obj['type'] == "CustomNebulae":
                    # generated_object = CustomNebula(xc=x, yc=y,)
                    log.warning("Custom Nebulae will be added soon")
                    continue
                else:
                    log.warning("Unrecognized type of the nebula #{0}: skip this one".format(ind_obj))
                    continue
            if cloudy_model_index is not None:
                if cur_obj['linerat_constant']:
                    lr = "Constant"
                else:
                    lr = "Variable"
                add_fits_kw = {"Model_ID": cloudy_model_id, "LineRat": lr}
            else:
                add_fits_kw = {}

            add_fits_kw["Distance"] = (cur_obj['distance'].to_value(u.kpc), "Distance to the nebula, kpc")

            continuum = None
            if cur_obj['continuum_type'] is not None and cur_obj['continuum_data'] is not None \
                    and cur_obj['continuum_type'].lower() in ['bb', 'poly', 'model']:

                if cur_obj['continuum_type'].lower() == 'model':
                    if isinstance(cur_obj['continuum_data'], dict) and ('wl' in cur_obj['continuum_data']) and \
                            ('flux' in cur_obj['continuum_data']):
                        if len(cur_obj['continuum_data']['wl']) != len(cur_obj['continuum_data']['flux']):
                            log.error("Number of wavelength and flux points is inconsistent for continuum")
                        else:
                            wlscale = cur_obj['continuum_data']['wl']
                            continuum = np.vstack([wlscale, cur_obj['continuum_data']['flux']])
                    elif ~isinstance(cur_obj['continuum_data'], dict) and lvmdatasimulator.CONTINUUM_MODELS is not None:
                        with fits.open(lvmdatasimulator.CONTINUUM_MODELS) as hdu:
                            if cur_obj['continuum_data'] >= hdu[0].data.shape[0]:
                                log.warning("Wrong continuum model ID for nebula #{0}".format(obj_id))
                            else:
                                wlscale = (np.arange(hdu[0].data.shape[1]) - hdu[0].header['CRPIX1'] + 1
                                           ) * hdu[0].header['CDELT1'] + hdu[0].header['CRVAL1']
                                continuum = np.vstack([wlscale, hdu[0].data[cur_obj['continuum_data']]])
                elif cur_obj['continuum_type'].lower() in ['poly', 'bb']:
                    continuum = cur_obj['continuum_data']
            if continuum is not None:
                if add_fits_kw is None:
                    add_fits_kw = {}
                add_fits_kw['CONTTYPE'] = (cur_obj['continuum_type'], "Type of the continuum")
                if cur_obj['continuum_flux'] > 0:
                    contflux = cur_obj['continuum_flux'].to_value(u.erg / u.cm ** 2 / u.s / u.arcsec ** 2 / u.AA)
                else:
                    contflux = 0
                add_fits_kw['CONTFLUX'] = (contflux,
                                           "Continuum brightness (in erg/s/cm^2/asec^2/AA) at ref. wl/Filter")
                if cur_obj['continuum_mag'] is not None:
                    contmag = cur_obj['continuum_mag'].to_value(u.mag)
                else:
                    contmag = None
                add_fits_kw['CONTMAG'] = (contmag,
                                          "Continuum brightness (in mag/asec^2) at ref. wl/Filter")
                if isinstance(cur_obj['continuum_wl'], str):
                    cont_wl = cur_obj['continuum_wl']
                else:
                    cont_wl = cur_obj['continuum_wl'].to_value(u.AA)
                add_fits_kw['CONTWL'] = (cont_wl, 'Reference wavelength/filter for cont. flux/mag')

            if cur_obj.get('max_extinction') > 0:
                if add_fits_kw is None:
                    add_fits_kw = {}
                add_fits_kw['EXT_LAW'] = (cur_obj['ext_law'], "Extinction law according to pyneb list")
                add_fits_kw['EXT_RV'] = (cur_obj['ext_rv'], "R_V value for extinction calculations")

            self.add_nebula(generated_object, obj_id=obj_id, zorder=cur_obj.get('zorder'), add_fits_kw=add_fits_kw,
                            continuum=continuum)
            obj_id += 1
        bar.finish()
        if (obj_id - obj_id_ini) == 0:
            return None
        else:
            return True

    def load_nebulae(self, file):
        """
        Load previously saved fits-file containing the information about all nebulae.
        Note: Grid should be equal to that used when fits-file was generated!
        """
        if not os.path.isfile(file):
            log.warning("ISM doesn't contain any nebula")
            return None
        else:
            with fits.open(file) as hdu:
                wcs = WCS(hdu[0].header)
                cdelt_file = [cdelt.to(u.degree).value for cdelt in wcs.proj_plane_pixel_scales()]
                cdelt_ism = [cdelt.to(u.degree).value for cdelt in self.wcs.proj_plane_pixel_scales()]
                check = ~np.isclose(cdelt_file, cdelt_ism)
                check = np.append(check, ~np.isclose(wcs.wcs.crval, self.wcs.wcs.crval))
                if any(np.array([cur_hdu.header.get('EXTNAME') for cur_hdu in hdu
                                 if (cur_hdu.header.get('EXTNAME') is not None) and
                                    ("_LINEPROFILE" in cur_hdu.header.get('EXTNAME'))])):
                    check = np.append(check, ~np.isclose(hdu[0].header.get("VELRES"), self.vel_resolution.value))
                if any(check):
                    log.warning("Grid of fits file is inconsistent with that defined in ISM")
                    return None
                self.content = fits.HDUList()
                for cur_ind, hh in enumerate(hdu):
                    if cur_ind == 0:
                        self.content.append(fits.PrimaryHDU(header=hh.header))
                    else:
                        self.content.append(fits.ImageHDU(header=hh.header, data=hh.data,
                                                          name=hh.header.get('EXTNAME')))
                return True

    def calc_extinction(self, wavelength=6562.81, x0=0, y0=0, xs=None, ys=None, extension_name=None):
        """
        Calculate coefficient to reduce flux due to extinction at given wavelength(s)

        Args:
            x0: start x-coordinate in the field of view for calculations
            y0: start x-coordinate in the field of view for calculations
            xs: x-size (in pixels) of the area for calculations (if None => then just pixel x0,y0 is considered; xs=1)
            ys: y-size (in pixels) of the area for calculations (if None => then just pixel x0,y0 is considered; ys=1)
            wavelength: in angstrom, particular wavelength (or wavelengths)
                at which the calculations should be performed
            extension_name (str): name of the extension for current dark nebula

        Returns:
            None (if no any dark nebula at particular location) or np.array of (nlines, ys, xs) shape
        """
        if self.content[0].header['Nobj'] == 0 or (x0 > self.width) or (y0 > self.height):
            return None
        if xs is None:
            x1 = x0
        else:
            x1 = np.clip(x0 + xs - 1, 0, self.width - 1)
        if ys is None:
            y1 = y0
        else:
            y1 = np.clip(y0 + ys - 1, 0, self.height - 1)

        def check_in_region(reg_x0, reg_y0, reg_xs, reg_ys):
            if ((reg_x0 + reg_xs) < x0) or (reg_x0 > x1) or ((reg_y0 + reg_ys) < y0) or (reg_y0 > y1):
                return False
            else:
                return True

        if extension_name is None:
            all_dark_nebulae = [hdu.header.get('EXTNAME') for hdu in self.content if
                                hdu.header.get('EXTNAME') is not None and
                                ("BRIGHTNESS" in hdu.header.get('EXTNAME')) and
                                hdu.header.get('DARK') and
                                check_in_region(hdu.header.get('X0'), hdu.header.get('Y0'),
                                                hdu.header.get('NAXIS1'), hdu.header.get('NAXIS2'))]
        else:
            if not check_in_region(self.content[extension_name].header.get('X0'),
                                   self.content[extension_name].header.get('Y0'),
                                   self.content[extension_name].header.get('NAXIS1'),
                                   self.content[extension_name].header.get('NAXIS2')):
                return None
            all_dark_nebulae = [extension_name]
        if len(all_dark_nebulae) == 0:
            return None
        if type(wavelength) in [float, int, np.float64, np.float32]:
            wavelength = np.atleast_1d(wavelength)
        ext_map = np.ones(shape=(len(wavelength), y1 - y0 + 1, x1 - x0 + 1), dtype=np.float32)

        for dark_nebula in all_dark_nebulae:
            cur_neb_av = np.zeros(shape=(y1 - y0 + 1, x1 - x0 + 1), dtype=np.float32)
            cur_neb_x0 = np.clip(x0 - self.content[dark_nebula].header.get('X0'), 0, None)
            cur_neb_y0 = np.clip(y0 - self.content[dark_nebula].header.get('Y0'), 0, None)
            cur_neb_x1 = self.content[dark_nebula].header.get('NAXIS1') - 1 - np.clip(
                self.content[dark_nebula].header.get('X0') + self.content[dark_nebula].header.get('NAXIS1') - 1 - x1,
                0, None)
            cur_neb_y1 = self.content[dark_nebula].header.get('NAXIS2') - 1 - np.clip(
                self.content[dark_nebula].header.get('Y0') + self.content[dark_nebula].header.get('NAXIS2') - 1 - y1,
                0, None)
            cur_neb_av[cur_neb_y0 + self.content[dark_nebula].header.get('Y0') - y0:
                       cur_neb_y1 + self.content[dark_nebula].header.get('Y0') - y0 + 1,
                       cur_neb_x0 + self.content[dark_nebula].header.get('X0') - x0:
                       cur_neb_x1 + self.content[dark_nebula].header.get('X0') - x0 + 1
                       ] = self.content[dark_nebula].data[cur_neb_y0: cur_neb_y1 + 1, cur_neb_x0: cur_neb_x1 + 1]
            cur_extinction_law = self.ext_law
            cur_r_v = self.R_V
            if self.content[dark_nebula].header.get('EXT_LAW'):
                cur_extinction_law = self.content[dark_nebula].header.get('EXT_LAW')
            if self.content[dark_nebula].header.get('EXT_RV'):
                cur_r_v = self.content[dark_nebula].header.get('EXT_RV')
            ext_map = ext_map * ism_extinction(av=cur_neb_av, r_v=cur_r_v,
                                               ext_law=cur_extinction_law, wavelength=wavelength)
        return ext_map

    def get_map(self, wavelength=6562.81, get_continuum=False):
        """
        Method to produce 2D map of all ISM nebulae in the selected line
        Args:
            wavelength (float or iterative): wavelength (in Angstrom) according to the lines list, or wavelength range
            get_continuum (bool): if True, then also counts the flux from the continuum
        """
        wavelength = np.atleast_1d(wavelength)
        if len(wavelength) == 1:
            wavelength = np.array([wavelength[0]-0.01, wavelength[0]+0.01])
        if self.content[0].header['Nobj'] == 0:
            log.warning("ISM doesn't contain any nebula")
            return None
        all_extensions = [hdu.header.get('EXTNAME') for hdu in self.content]
        all_extensions_brt = np.array([extname for extname in all_extensions
                                       if extname is not None and ("BRIGHTNESS" in extname) and
                                       check_overlap(self.content[extname], (self.height, self.width))])

        if all([self.content[cur_ext].header.get("DARK") for cur_ext in all_extensions_brt]):
            # !!!! ADD later accounting of the continuum and extinction from those nebulae
            log.warning("ISM doesn't contain any emission nebula")
            return None

        all_extensions_brt = all_extensions_brt[
            np.argsort([self.content[cur_ext].header.get('ZORDER') for cur_ext in all_extensions_brt])]

        map_2d = np.zeros(shape=(self.height, self.width), dtype=np.float32)
        map_is_empty = True
        for cur_ext in all_extensions_brt:
            my_comp = "_".join(cur_ext.split("_")[:2])

            y0_in_field = np.clip(self.content[cur_ext].header['Y0'], 0, None)
            y1_in_field = np.clip(self.content[cur_ext].header['Y0'] + self.content[cur_ext].header['NAXIS2'] - 1, None,
                                  self.height-1)
            x0_in_field = np.clip(self.content[cur_ext].header['X0'], 0, None)
            x1_in_field = np.clip(self.content[cur_ext].header['X0'] + self.content[cur_ext].header['NAXIS1'] - 1, None,
                                  self.width - 1)

            if get_continuum and (my_comp + "_CONTINUUM" in all_extensions):
                n_wl_bins = int(np.clip(np.ceil((wavelength[1]-wavelength[0])/20.), 10, 200))
                wl_grid = np.linspace(wavelength[0], wavelength[1], n_wl_bins)
                continuum = np.sum(self._get_continuum(my_comp, wl_grid))
                add_continuum = self.content[cur_ext].data / np.max(self.content[cur_ext].data) * continuum
                map_2d[y0_in_field: y1_in_field + 1, x0_in_field: x1_in_field + 1] += \
                    add_continuum[y0_in_field - self.content[cur_ext].header['Y0']:
                                  y1_in_field - self.content[cur_ext].header['Y0'] + 1,
                                  x0_in_field - self.content[cur_ext].header['X0']: x1_in_field -
                                  self.content[cur_ext].header['X0'] + 1
                                  ]
                map_is_empty = False
            if self.content[cur_ext].header.get("DARK"):
                if map_is_empty:
                    continue
                ext_map = self.calc_extinction(wavelength=(wavelength[-1] + wavelength[0])/2., xs=self.width,
                                               ys=self.height,
                                               extension_name=cur_ext)
                if ext_map is not None:
                    map_2d = map_2d * ext_map[0]
                continue

            all_flux_wl = [extname[-7:].strip('_') for extname in all_extensions
                           if extname is not None and (my_comp in extname and "FLUX_" in extname)]
            all_flux_wl = np.array(all_flux_wl)

            # flux_ext = [extname for extname in all_extensions
            #             if extname is not None and (my_comp in extname and
            #                                         "FLUX_{0}".format(np.round(wavelength, 2)) in extname)]
            add_emission = np.zeros(shape=(self.content[cur_ext].header['NAXIS2'],
                                           self.content[cur_ext].header['NAXIS1']), dtype=np.float32)
            if len(all_flux_wl) == 0:
                fluxrat_ext = [extname for extname in all_extensions
                               if extname is not None and (my_comp in extname and "FLUXRATIOS" in extname)]
                if len(fluxrat_ext) == 0:
                    continue
                fluxrat_ext = fluxrat_ext[0]

                wl_indexes = np.flatnonzero((self.content[fluxrat_ext].data[0, :] > (wavelength[0])) &
                                            (self.content[fluxrat_ext].data[0, :] < (wavelength[1])))

                # wl_index = np.flatnonzero(np.isclose(self.content[fluxrat_ext].data[0, :], wavelength))
                if len(wl_indexes) == 0:
                    continue
                for wl_index in wl_indexes:
                    add_emission += (self.content[cur_ext].data * self.content[fluxrat_ext].data[1, wl_index])
            else:
                all_flux_wl_float = np.array(all_flux_wl).astype(np.float32)
                wl_indexes = np.flatnonzero((all_flux_wl_float > (wavelength[0] - 0.01)) &
                                            (all_flux_wl_float < (wavelength[1] + 0.01)))
                flux_ext_wl = all_flux_wl[wl_indexes]
                if len(wl_indexes) == 0:
                    continue

                for cur_wl in flux_ext_wl:
                    add_emission += self.content[my_comp + "_FLUX_" + cur_wl].data

            map_2d[y0_in_field: y1_in_field + 1, x0_in_field: x1_in_field + 1] += \
                add_emission[y0_in_field - self.content[cur_ext].header['Y0']:
                             y1_in_field - self.content[cur_ext].header['Y0'] + 1,
                             x0_in_field - self.content[cur_ext].header['X0']: x1_in_field -
                             self.content[cur_ext].header['X0'] + 1
                             ]

            map_is_empty = False
        return map_2d * (proj_plane_pixel_scales(self.wcs)[0] * 3600) ** 2

    def get_spectrum(self, wl_grid=None, aperture_mask=None, fibers_coords=None, pix_size=1):
        if aperture_mask is None or (np.sum(aperture_mask) == 0) or (self.content[0].header['Nobj'] == 0):
            log.warning("No overlapping detected between the ISM component and the fibers => no spectra extraction")
            return None
        all_extensions = [hdu.header.get('EXTNAME') for hdu in self.content]
        all_extensions_brt = np.array([extname for extname in all_extensions
                                       if extname is not None and ("BRIGHTNESS" in extname) and
                                       check_overlap(self.content[extname], (self.height, self.width))])
        if all([self.content[cur_ext].header.get("DARK") for cur_ext in all_extensions_brt]):
            return None
        all_extensions_brt = all_extensions_brt[np.argsort([self.content[cur_ext].header.get('ZORDER')
                                                            for cur_ext in all_extensions_brt])]

        wl_logscale = np.log(wl_grid.value)
        wl_logscale_highres = np.arange((np.round(wl_logscale[-1] - wl_logscale[0]) * 1e6
                                         ).astype(int)) * 1e-6 + np.round(wl_logscale[0], 6)
        delta_lr = np.roll(wl_logscale, -1) - wl_logscale
        delta_lr[-1] = delta_lr[-2]

        xx, yy = np.meshgrid(np.arange(aperture_mask.shape[1]), np.arange(aperture_mask.shape[0]))
        pix_in_apertures = aperture_mask > 0
        xstart = np.clip(np.min(xx[pix_in_apertures]) - 1, 0, None)
        ystart = np.clip(np.min(yy[pix_in_apertures]) - 1, 0, None)
        xfin = np.clip(np.max(xx[pix_in_apertures]) + 1, None, aperture_mask.shape[1]-1)
        yfin = np.clip(np.max(yy[pix_in_apertures]) + 1, None, aperture_mask.shape[0]-1)
        aperture_mask_sub = aperture_mask[ystart: yfin + 1, xstart: xfin + 1]
        xx_sub, yy_sub = np.meshgrid(np.arange(xfin - xstart + 1), np.arange(yfin - ystart + 1))
        n_apertures = np.max(aperture_mask)
        # aperture_centers = np.round(fibers_coords).astype(int)
        spectrum = np.zeros(shape=(n_apertures, len(wl_grid)), dtype=np.float32)

        fiber_radius = np.median(fibers_coords[:, 2])
        kern_mask = calc_circular_mask(fiber_radius)
        kern = kernels.CustomKernel(kern_mask.reshape((1, kern_mask.shape[0], kern_mask.shape[1])))

        bar = progressbar.ProgressBar(max_value=len(all_extensions_brt)).start()
        for neb_index, cur_ext in enumerate(all_extensions_brt):
            cur_neb_in_mask = np.zeros_like(xx_sub)
            y0 = self.content[cur_ext].header.get("Y0")
            x0 = self.content[cur_ext].header.get("X0")
            nx = self.content[cur_ext].header.get('NAXIS1')
            ny = self.content[cur_ext].header.get('NAXIS2')
            cur_neb_in_mask[(xx_sub >= (x0 - xstart)) & (xx_sub <= (x0 + nx - xstart)) &
                            (yy_sub >= (y0 - ystart)) & (yy_sub <= (y0 + ny - ystart))] = True
            cur_neb_in_mask_ap = cur_neb_in_mask * (aperture_mask_sub > 0)
            if not np.sum(cur_neb_in_mask_ap):
                bar.update(neb_index + 1)
                continue
            cur_mask_in_neb = np.zeros(shape=self.content[cur_ext].data.shape, dtype=bool)
            xx_neb, yy_neb = np.meshgrid(np.arange(self.content[cur_ext].data.shape[1]),
                                         np.arange(self.content[cur_ext].data.shape[0]))
            cur_mask_in_neb[(xx_neb >= (xstart - x0)) & (xx_neb <= (xfin - x0)) &
                            (yy_neb >= (ystart - y0)) & (yy_neb <= (yfin - y0))] = True
            xstart_neb = np.min(xx_neb[cur_mask_in_neb])
            ystart_neb = np.min(yy_neb[cur_mask_in_neb])
            xfin_neb = np.max(xx_neb[cur_mask_in_neb])
            yfin_neb = np.max(yy_neb[cur_mask_in_neb])
            selected_apertures = np.flatnonzero(((fibers_coords[:, 0] + fibers_coords[:, 2]) >= x0) &
                                                ((fibers_coords[:, 0] - fibers_coords[:, 2]) <= (x0 + nx - 1)) &
                                                ((fibers_coords[:, 1] + fibers_coords[:, 2]) >= y0) &
                                                ((fibers_coords[:, 1] - fibers_coords[:, 2]) <= (y0 + ny - 1)))
            selected_apertures = np.array([sa for sa in selected_apertures if (sa+1) in aperture_mask_sub], dtype=int)
            if len(selected_apertures) == 0:
                bar.update(neb_index + 1)
                continue

            # Here I check if it is necessary to extend all involved arrays to account for fibers at the edges of neb.
            dx0 = abs(np.clip(np.round(np.min(fibers_coords[selected_apertures, 0] -
                                              fibers_coords[selected_apertures, 2])).astype(int) - 1 - x0, None, 0))
            dx1 = np.clip(np.round(np.max(fibers_coords[selected_apertures, 0] +
                                          fibers_coords[selected_apertures, 2])).astype(int) + 2 - x0 - nx, 0, None)
            dy0 = abs(np.clip(np.round(np.min(fibers_coords[selected_apertures, 1] -
                                              fibers_coords[selected_apertures, 2])).astype(int) - 1 - y0, None, 0))
            dy1 = np.clip(np.round(np.max(fibers_coords[selected_apertures, 1] +
                                          fibers_coords[selected_apertures, 2])).astype(int) + 2 - y0 - ny, 0, None)
            if (dx0 > 0) or (dy0 > 0) or (dx1 > 0) or (dy1 > 0):
                npad = ((0, 0), (dy0, dy1), (dx0, dx1))
            else:
                npad = None

            selected_y = fibers_coords[selected_apertures, 1] - ystart_neb - y0 + dy0
            selected_x = fibers_coords[selected_apertures, 0] - xstart_neb - x0 + dx0

            if self.content[cur_ext].header.get("DARK"):
                extinction_map = self.content[cur_ext].data[cur_mask_in_neb].reshape((1, yfin_neb - ystart_neb + 1,
                                                                                      xfin_neb - xstart_neb + 1))
                if npad is not None:
                    extinction_map = np.pad(extinction_map, pad_width=npad, mode='constant', constant_values=0)

                if self.content[cur_ext].header.get('EXT_LAW'):
                    cur_extinction_law = self.content[cur_ext].header.get('EXT_LAW')
                else:
                    cur_extinction_law = self.ext_law
                if self.content[cur_ext].header.get('EXT_RV'):
                    cur_r_v = self.content[cur_ext].header.get('EXT_RV')
                else:
                    cur_r_v = self.R_V

                data_in_apertures = convolve_array(extinction_map,
                                                   kern, selected_x, selected_y, pix_size,
                                                   nchunks=self.content[cur_ext].header.get("NCHUNKS"),
                                                   normalize_kernel=True)

                # data_in_apertures = convolve_fft(extinction_map, kern,
                #                                  allow_huge=True, normalize_kernel=True)[
                #                     :, aperture_centers[selected_apertures, 1] - ystart_neb - y0 + dy0,
                #                     aperture_centers[selected_apertures, 0] - xstart_neb - x0 + dx0]
                data_in_apertures = data_in_apertures.reshape((data_in_apertures.shape[0] * data_in_apertures.shape[1],
                                                               1))
                spectrum[selected_apertures, :] = \
                    spectrum[selected_apertures, :] * ism_extinction(av=data_in_apertures, r_v=cur_r_v,
                                                                     ext_law=cur_extinction_law, wavelength=wl_grid).T

                bar.update(neb_index + 1)
                continue

            my_comp = "_".join(cur_ext.split("_")[:2])

            if self.content[cur_ext].header.get('MAXBRT'):
                if self.content[cur_ext].header.get("LINERAT") == 'Variable':
                    all_wavelength = np.array([extname.split("_")[-1] for extname in all_extensions
                                               if extname is not None and (my_comp + "_FLUX_" in extname)])
                else:
                    all_wavelength = self.content[my_comp + "_FLUXRATIOS"].data[0, :]

                if self.content[cur_ext].header.get("LINERAT") == 'Variable':
                    all_fluxes = np.array([self.content[my_comp + "_FLUX_" + wl].data[
                                               cur_mask_in_neb].reshape((yfin_neb - ystart_neb + 1,
                                                                         xfin_neb - xstart_neb + 1))
                                           for wl in all_wavelength], dtype=np.float32)
                    all_wavelength = all_wavelength.astype(np.float32)
                else:
                    all_fluxes = self.content[cur_ext].data[cur_mask_in_neb].reshape((1, yfin_neb - ystart_neb + 1,
                                                                                      xfin_neb - xstart_neb + 1)
                                                                                     ).astype(np.float32)

                if my_comp + "_LINEPROFILE" in self.content:
                    lsf = self.content[my_comp + "_LINEPROFILE"].data[
                          :, cur_mask_in_neb].reshape((len(self.vel_grid), yfin_neb - ystart_neb + 1,
                                                       xfin_neb - xstart_neb + 1)).astype(np.float32)
                else:
                    if my_comp + "_VEL" in self.content:
                        vel = self.content[my_comp + "_VEL"].data[cur_mask_in_neb].reshape((yfin_neb - ystart_neb + 1,
                                                                                            xfin_neb - xstart_neb + 1))
                    else:
                        vel = np.zeros(shape=(yfin_neb - ystart_neb + 1, xfin_neb - xstart_neb + 1),
                                       dtype=np.float32)  # + self.content[cur_ext].header.get('SysVel')
                        # self.sys_velocity.value
                    if my_comp + "_DISP" in self.content:
                        disp = self.content[my_comp + "_DISP"].data[
                            cur_mask_in_neb].reshape((yfin_neb - ystart_neb + 1, xfin_neb - xstart_neb + 1))
                    else:
                        disp = np.zeros(shape=(yfin_neb - ystart_neb + 1, xfin_neb - xstart_neb + 1),
                                        dtype=np.float32) + self.content[cur_ext].header.get('TurbVel')
                        # self.turbulent_sigma.value
                    lsf = np.exp(-np.power(
                        (self.vel_grid.value[:, None, None] - vel[None, :, :]) / disp[None, :, :], 2.) / 2)
                    lsf = (lsf / np.nansum(lsf, axis=0)).astype(np.float32)

                if npad is not None:
                    lsf = np.pad(lsf, pad_width=npad, mode='constant', constant_values=0)
                    all_fluxes = np.pad(all_fluxes, pad_width=npad, mode='constant', constant_values=0)

                data_in_apertures = [convolve_array(lsf * line_data[None, :, :],
                                                    kern, selected_x, selected_y, pix_size,
                                                    nchunks=self.content[cur_ext].header.get("NCHUNKS"))
                                     for line_data in all_fluxes]
                data_in_apertures = np.array(data_in_apertures)
                if data_in_apertures.shape == 2:
                    data_in_apertures = data_in_apertures.reshape((1, data_in_apertures.shape[0],
                                                                   data_in_apertures.shape[1]))

                # if all_fluxes.shape[0] == 1:
                #     # selected_y = fibers_coords[selected_apertures, 1] - ystart_neb - y0 + dy0
                #     # selected_x = fibers_coords[selected_apertures, 0] - xstart_neb - x0 + dx0
                #     data_in_apertures = convolve_array(lsf * all_fluxes[0][None, :, :],
                #                                        kern, selected_x, selected_y, pix_size,
                #                                        nchunks=self.content[cur_ext].header.get("NCHUNKS"))
                #     data_in_apertures = data_in_apertures.reshape((1, data_in_apertures.shape[0],
                #                                                   data_in_apertures.shape[1]))
                #     # hdu = fits.PrimaryHDU(data=data_in_apertures)
                #     # hdu.writeto(f'./data_{neb_index}.fits', overwrite=True)
                #
                # else:
                #
                #     # selected_y = fibers_coords[selected_apertures, 1] - ystart_neb - y0 + dy0
                #     # selected_x = fibers_coords[selected_apertures, 0] - xstart_neb - x0 + dx0
                #     # data_in_apertures = Parallel(n_jobs=lvmdatasimulator.n_process)(
                #     #     delayed(convolve_array)(lsf * line_data[None, :, :], kern,
                #     #                            selected_y, selected_x, pix_size)
                #     #     for line_data in all_fluxes)
                #     data_in_apertures = [convolve_array(lsf * line_data[None, :, :], kern,
                #                                         selected_y, selected_x, pix_size,
                #                                         nchunks=self.content[cur_ext].header.get("NCHUNKS"))
                #                          for line_data in all_fluxes]
                #     data_in_apertures = np.array(data_in_apertures)
                #     # hdu = fits.PrimaryHDU(data=data_in_apertures)
                #     # hdu.writeto(f'./data_{neb_index}.fits', overwrite=True)

                data_in_apertures = np.moveaxis(data_in_apertures, 2, 0)
                if data_in_apertures.shape[1] > 1:
                    prf_index = np.flatnonzero(all_wavelength == 6562.81)
                else:
                    prf_index = 0
                flux_norm_in_apertures = data_in_apertures.sum(axis=2)
                line_prf_in_apertures = data_in_apertures[:, prf_index, :].reshape(
                    (data_in_apertures.shape[0], data_in_apertures.shape[2])) / \
                    flux_norm_in_apertures[:, prf_index].reshape(data_in_apertures.shape[0], 1)

                wl_logscale_lsf = np.log((self.vel_grid.value + self.content[cur_ext].header.get('SysVel')
                                          ) / 2.9979e5 + 1)
                highres_factor = 1e6
                wl_logscale_lsf_highres = np.arange(np.round((wl_logscale_lsf[-1] - wl_logscale_lsf[0]
                                                              ) / 2. * highres_factor).astype(int) * 2 + 1
                                                    ) / highres_factor + wl_logscale_lsf[0]
                p = interp1d(wl_logscale_lsf, line_prf_in_apertures, axis=1, assume_sorted=True)
                line_highres_log = p(wl_logscale_lsf_highres)
                line_highres_log = line_highres_log / np.sum(line_highres_log, axis=1)[:, None]
                if flux_norm_in_apertures.shape[1] == 1:
                    flux_norm_in_apertures = flux_norm_in_apertures * \
                                        self.content[my_comp + "_FLUXRATIOS"].data[1, None, :]
                wl_indexes = np.round((np.log(all_wavelength * (1 + self.content[cur_ext].header.get('SysVel') /
                                                                2.9979e5)) - wl_logscale_highres[0]
                                       ) * highres_factor).astype(int)
                rec = (wl_indexes > 0) & (wl_indexes < len(wl_logscale_highres))
                spectrum_highres_log = np.zeros(shape=(len(selected_apertures), len(wl_logscale_highres)),
                                                dtype=np.float32)
                win = (len(wl_logscale_lsf_highres) - 1) // 2
                for ind, r in enumerate(rec):
                    if r:
                        spectrum_highres_log[:, wl_indexes[ind] - win: wl_indexes[ind] + win + 1] += \
                            line_highres_log[:, :] * flux_norm_in_apertures[:, ind].reshape(
                                (len(selected_apertures), 1))
                p = interp1d(wl_logscale_highres, spectrum_highres_log, axis=1, assume_sorted=True, bounds_error=False,
                             fill_value='extrapolate')
                spectrum[selected_apertures, :] += (p(wl_logscale) * delta_lr * highres_factor)

            if my_comp + "_CONTINUUM" in self.content:
                brt_max = self.content[cur_ext].header.get('MAXBRT')
                if not brt_max:
                    brt_max = self.content[cur_ext].header.get('MAXEXT')
                if not brt_max:
                    brt_max = 1
                continuum = self._get_continuum(my_comp, wl_grid.value)

                brt = self.content[cur_ext].data[cur_mask_in_neb].reshape((1, yfin_neb - ystart_neb + 1,
                                                                           xfin_neb - xstart_neb + 1)) / brt_max
                if brt.shape[0] != 1:
                    brt = brt.reshape((1, brt.shape[0], brt.shape[1]))
                if npad is not None:
                    brt = np.pad(brt, pad_width=npad, mode='constant', constant_values=0)
                # selected_y = fibers_coords[selected_apertures, 1] - ystart_neb - y0 + dy0
                # selected_x = fibers_coords[selected_apertures, 0] - xstart_neb - x0 + dx0
                data_in_apertures = convolve_array(brt,
                                                   kern, selected_x, selected_y, pix_size,
                                                   nchunks=self.content[cur_ext].header.get("NCHUNKS"))
                # data_in_apertures = convolve_fft(brt, kern,
                #                                  allow_huge=True, normalize_kernel=False)[
                #                     :, aperture_centers[selected_apertures, 1] - ystart_neb - y0 + dy0,
                #                     aperture_centers[selected_apertures, 0] - xstart_neb - x0 + dx0]
                data_in_apertures = data_in_apertures.reshape((data_in_apertures.shape[0] * data_in_apertures.shape[1],
                                                               1))
                spectrum[selected_apertures, :] += continuum[None, :] * data_in_apertures

            bar.update(neb_index + 1)
        bar.finish()
        return spectrum * (proj_plane_pixel_scales(self.wcs)[0] * 3600) ** 2 * fluxunit * u.arcsec ** 2

# encoding: utf-8 # # @Author: <NAME>, <NAME> # @Date: Nov 15, 2021 # @Filename: ism.py # @License: BSD 3-Clause # @Copyright: <NAME>, <NAME> Auxiliary function producing the inhomogeneities on the brightness distribution for the Cloud of Bubble objects     using the spherical harmonics. Auxiliary function computing the brightness of the Cloud or Bubble at given radii and in given line     according to the Cloudy models Auxiliary function to project the brightness or velocities from the spherical to cartesian coordinates Auxiliary function to limit the angle values to the range of [0, pi] Auxiliary function to map the coordinates from cartesian to spherical system Checks the input parameters of the pre-computed Cloudy model and return corresponding index in the grid # # for cur_ext in range(len(hdu)): #     if cur_ext == 0: #         continue #     found = False #     for p in params: #         if p == 'id': #             continue #         precision = 1 #         if p == 'Z': #             precision = 2 #         if np.round(params[p], precision) != np.round(hdu[cur_ext].header[p], precision): #             break #     else: #         found = True #     if found: #         return cur_ext, check_id # check_id = lvmdatasimulator.CLOUDY_SPEC_DEFAULTS['id'] # log.warning('Input parameters do not correspond to any pre-computed Cloudy model.' #             'Default Cloudy model will be used (id = {0})'.format(check_id)) Base class defining properties of every nebula type.     By itself it describes the rectangular nebula (e.g. DIG)     Constructed nebula has 4 dimensions, where 4th derive its appearance in different lines         (if spectrum_id is None, or if it is dark nebula => only one line) # Center of the region in the field of view, pix # Center of the region in the field of view, pix # Coordinates of the bottom-left corner in the field of view, pix # Coordinates of the bottom-left corner in the field of view, pix # full width of cartesian grid, pix (should be odd) # full height of cartesian grid, pix (should be odd) # width of the nebula in pc (not used if pix_width is set up) # height of the nebula in pc (not used if pix_height is set up) # pixel size in pc # ID of a template Cloudy emission spectrum for this nebula # limit the number of the lines to the first N brightest # Systemic velocity # Velocity dispersion due to turbulence; included in calculations of LSF # Spatial scale of correlated perturbations # Maximal amplitude of perturbations # full size along line of sight in pixels # number of chuncks to use for the convolution. If negative, select automatically # velocity gradient along the nebula # Position angle of the kinematical axis (for the velocity gradient or rotation velocity) # True if the ratio of line fluxes shouldn't change across the nebula Method to obtain the brightness (or density) distribution of the nebula in cartesian coordinates Derive the brightness (or density) distribution of the nebula for each emission line in cartesian coordinates Project the 3D nebula onto sky plane (for emission or continuum sources) Project the 3D emission nebula line onto sky plane (return images in each emission line) Project the 3D nebula onto sky plane (for dark clouds) # if vel_field is None: #     return np.atleast_1d(self.sys_velocity) # else: #     return vel_field + self.sys_velocity # @cached_property # def line_profile(self): #     lprf = np.zeros(shape=len(self.los_velocity), dtype=float) #     lprf[np.floor(len(lprf) / 2.).astype(int)] = 1. #     return lprf Class defining a simple rectangular component.     This is equal to Nebula, but no perturbations and turbulence by default # Maximal amplitude of perturbations # Velocity dispersion due to turbulence; included in calculations of LSF # True if the ratio of line fluxes shouldn't change across the nebula Class defining a simple elliptical component.     No perturbations and turbulence by default # Maximal amplitude of perturbations # Velocity dispersion due to turbulence; included in calculations of LSF # Radius along the major axis of the ellipse (or radius of the circle) # position angle of the major axis # ratio of minor/major axes # True if the ratio of line fluxes shouldn't change across the nebula Method to obtain the brightness (or density) distribution of the nebula in cartesian coordinates Class defining a simple circular component. # True if the ratio of line fluxes shouldn't change across the nebula Class of an isotropic cylindrical shape filament.     Defined by its position, length, PA, radius, maximal optical depth.     If it is emission-type filament, then also maximal brightness is required.     Velocity gradient also can be set up # position angle of the filament # full length of the filament # full width (diameter) of the filament # True if the ratio of line fluxes shouldn't change across the nebula Method to obtain the brightness (or density) distribution of the nebula in cartesian coordinates Class of an isotropic cylindrical shape filament.     Defined by its position, length, PA, radius, maximal optical depth     if it is emission-type filament, then maximal brightness     NB: this class is obsolete, but might be considered later in case of implementation of varying line ratios # position angle of the filament # full length of the filament # full width (diameter) of the filament # velocity gradient along the filament (to be added) # True if the ratio of line fluxes shouldn't change across the nebula Method to calculate brightness (or opacity) of the cloud at given theta, phi and radii          theta: float -- azimuthal angle [0, 2 * np.pi]         rad: float -- radius [0, self.width / 2]         h: float -- height [0, self.length]         Returns:             3D cube of normalized brightness in theta-rad-h grid; total brightness = 1 Class defining the galaxy object (set up it as Sersic2D profile assuming it has continuum and emission components) # position angle of the major axis # ratio of minor/major axes # Effective radius in kpc # Maximum radius for calculations (in R_eff) # Sersic index # Rotational velocity (not implemented yet) # True if the ratio of line fluxes shouldn't change across the nebula Method to obtain the brightness (or density) distribution of the nebula in cartesian coordinates Class defining the DIG component. For now it is defined just by its brightness (constant) Class of an isotropic spherical gas cloud without any ionization source.     Defined by its position, radius, density, maximal optical depth # Degree of perturbations (max. degree of spherical harmonics for cloud) # True if the ratio of line fluxes shouldn't change across the nebula Method to calculate brightness (or opacity) of the cloud at given theta, phi and radii          theta: float -- polar angle [0, np.pi]         phi: float -- azimuthal angle [0, 2 * np.pi]         rad: float -- radius [0, self.radius]         Returns:             3D cube of normalized brightness in theta-phi-rad grid; total brightness = 1 Method to calculate brightness of the cloud at given theta, phi and radii for each line          theta: float -- polar angle [0, np.pi]         phi: float -- azimuthal angle [0, 2 * np.pi]         rad: float -- radius [0, self.radius]         Returns:             4D cube of brightness in line-theta-phi-rad grid; normalized to the total brightness in Halpha Class of an isotropic thin expanding bubble. Calculate line of sight velocity at given radius, phi, theta          V ~ 1/brightness (given that v~1/n_e^2 and brightness~ne^2) Line spread function as a function of coorinates, including the velocity center shift # mu = self.velocity(theta, phi) Returns local spectrum, per pc**3 of area Produces the distribution of the observed line profiles in each pixels of the sky plane Class defining the custom nebulae with the user-defined distribution of the brighntess, continuum         and line shapes in different lines. Class defining the ISM contribution to the field of view # Width of field of view in pixels # Width of field of view in pixels # Spectral resolution of the simulation # Minimal number of pixels for a resolution element at wl = 10000A for construction of vel.grid # Distance to the object for further conversion between arcsec and pc # Systemic velocity to center the vel.grid on # Maximal deviation from the systemic velocity to setup vel.grid # turbulence vel. disp. to be used for every nebula unless other specified # R_V value defining the reddening curve (to be used unless other value is provided for a nebula) # Reddening law (to be used unless other value is provided for a nebula) Properly extracts continuum for current nebula taking into account its shape and surface brightness         :param my_comp:         :param wl_grid:         :return: continuum Method to add the particular nebula to the ISM object and to the output multi-extensions fits file Generate all the Nebulae from the input list          Args:             all_objects: list -- contains a dictionary describing the nebula to add:                 example:                 all_objects = [{type: "Bubble, Filament, DIG, ....",                                 sys_velocity: 0 * u.km/u.s,                                 expansion_velocity: 30 * u.km/u.s,                                 turbulent_sigma: 10 * u.km/u.s,                                 radius: 5 * u.pc,                                 max_brightness: 1e-16 * u.erg / u.cm**2 / u.s / u.arcsec ** 2,                                 RA: "08h12m13s",                                 DEC: "-20d14m13s",                                 'perturb_degree': 8, # max. order of spherical harmonics to generate inhomogeneities                                 'perturb_amplitude': 0.1, # relative max. amplitude of inhomogeneities,                                 'perturb_scale': 200 * u.pc,  # spatial scale to generate inhomogeneities (DIG only)                                 'distance': 50 * u.kpc,  # distance to the nebula (default is from ISM)                                 'cloudy_id': None,  # id of pre-generated Cloudy model                                 'cloudy_params': {'Z': 0.5, 'qH': 49., 'nH': 10, 'Teff': 30000, 'Geometry': 'Sphere'}, #                                     parameters defining the pre-generated Cloudy model (used if spectrum_id is None)                                 'n_brightest_lines': 10,  # only this number of the brightest lines will be evaluated                                 'linerat_constant': False #  if True -> lines ratios don't vary across Cloud/Bubble                                 'continuum_type': 'BB' or 'Model' or 'Poly' # type of the continuum model                                 'continuum_data': model_id or [poly_coefficients] or Teff or dict with "wl" and "flux"                                                 # value defining cont. shape                                 'continuum_flux': 1e-16 * u.erg / u.cm ** 2 / u.s / u.arcsec **2 / u.AA,                                 'continuum_mag': 22 * u.mag,                                 'continuum_wl': 5500, # could be also R, V, B,                                 'ext_law': 'F99',  # Extinction law, one of those used by pyneb (used for dark nebulae)                                 'ext_rv': 3.1,  # Value of R_V for extinction curve calculation (used for dark nebulae)                                 'vel_gradient: 12. * u.km / u.s / u.pc  # Line-of-sight velocity gradient                                 'vel_pa: 30. * u.degree  # PA of kinematical axis (for vel_gradient or vel_rot)                                 }] # Setup default parameters for missing keywords #{0}: skip this one".format(ind_obj)) #{0}: " #{0}: skip this one".format(ind_obj)) # ==== Check input parameters and do necessary conversions #{0}: skip this one".format(ind_obj)) #{0}: skip this one".format(ind_obj)) #{0}: skip this one".format(ind_obj)) #{0}: skip this one".format(ind_obj)) #{0}: skip this one".format(ind_obj)) #{0}".format(ind_obj)) #{0}: set it to 1.".format(ind_obj)) #{0})." #{0}".format(ind_obj)) # ==== Start calculations of different nebulae # generated_object = CustomNebula(xc=x, yc=y,) #{0}: skip this one".format(ind_obj)) #{0}".format(obj_id)) Load previously saved fits-file containing the information about all nebulae.         Note: Grid should be equal to that used when fits-file was generated! Calculate coefficient to reduce flux due to extinction at given wavelength(s)          Args:             x0: start x-coordinate in the field of view for calculations             y0: start x-coordinate in the field of view for calculations             xs: x-size (in pixels) of the area for calculations (if None => then just pixel x0,y0 is considered; xs=1)             ys: y-size (in pixels) of the area for calculations (if None => then just pixel x0,y0 is considered; ys=1)             wavelength: in angstrom, particular wavelength (or wavelengths)                 at which the calculations should be performed             extension_name (str): name of the extension for current dark nebula          Returns:             None (if no any dark nebula at particular location) or np.array of (nlines, ys, xs) shape Method to produce 2D map of all ISM nebulae in the selected line         Args:             wavelength (float or iterative): wavelength (in Angstrom) according to the lines list, or wavelength range             get_continuum (bool): if True, then also counts the flux from the continuum # !!!! ADD later accounting of the continuum and extinction from those nebulae # flux_ext = [extname for extname in all_extensions #             if extname is not None and (my_comp in extname and #                                         "FLUX_{0}".format(np.round(wavelength, 2)) in extname)] # wl_index = np.flatnonzero(np.isclose(self.content[fluxrat_ext].data[0, :], wavelength)) # aperture_centers = np.round(fibers_coords).astype(int) # Here I check if it is necessary to extend all involved arrays to account for fibers at the edges of neb. # data_in_apertures = convolve_fft(extinction_map, kern, #                                  allow_huge=True, normalize_kernel=True)[ #                     :, aperture_centers[selected_apertures, 1] - ystart_neb - y0 + dy0, #                     aperture_centers[selected_apertures, 0] - xstart_neb - x0 + dx0] # + self.content[cur_ext].header.get('SysVel') # self.sys_velocity.value # self.turbulent_sigma.value # if all_fluxes.shape[0] == 1: #     # selected_y = fibers_coords[selected_apertures, 1] - ystart_neb - y0 + dy0 #     # selected_x = fibers_coords[selected_apertures, 0] - xstart_neb - x0 + dx0 #     data_in_apertures = convolve_array(lsf * all_fluxes[0][None, :, :], #                                        kern, selected_x, selected_y, pix_size, #                                        nchunks=self.content[cur_ext].header.get("NCHUNKS")) #     data_in_apertures = data_in_apertures.reshape((1, data_in_apertures.shape[0], #                                                   data_in_apertures.shape[1])) #     # hdu = fits.PrimaryHDU(data=data_in_apertures) #     # hdu.writeto(f'./data_{neb_index}.fits', overwrite=True) # # else: # #     # selected_y = fibers_coords[selected_apertures, 1] - ystart_neb - y0 + dy0 #     # selected_x = fibers_coords[selected_apertures, 0] - xstart_neb - x0 + dx0 #     # data_in_apertures = Parallel(n_jobs=lvmdatasimulator.n_process)( #     #     delayed(convolve_array)(lsf * line_data[None, :, :], kern, #     #                            selected_y, selected_x, pix_size) #     #     for line_data in all_fluxes) #     data_in_apertures = [convolve_array(lsf * line_data[None, :, :], kern, #                                         selected_y, selected_x, pix_size, #                                         nchunks=self.content[cur_ext].header.get("NCHUNKS")) #                          for line_data in all_fluxes] #     data_in_apertures = np.array(data_in_apertures) #     # hdu = fits.PrimaryHDU(data=data_in_apertures) #     # hdu.writeto(f'./data_{neb_index}.fits', overwrite=True) # selected_y = fibers_coords[selected_apertures, 1] - ystart_neb - y0 + dy0 # selected_x = fibers_coords[selected_apertures, 0] - xstart_neb - x0 + dx0 # data_in_apertures = convolve_fft(brt, kern, #                                  allow_huge=True, normalize_kernel=False)[ #                     :, aperture_centers[selected_apertures, 1] - ystart_neb - y0 + dy0, #                     aperture_centers[selected_apertures, 0] - xstart_neb - x0 + dx0]