exo_k.atm ========= .. py:module:: exo_k.atm .. autoapi-nested-parse:: @author: jeremy leconte This module contain classes to handle the radiative properties of atmospheres. This allows us to compute the transmission and emission spectra of those atmospheres. Physical properties of the atmosphere are handled in atm_profile.py. The nomenclature for layers, levels, etc, is as follows (example with 4 layers; Nlay=4): :: ------------------- Model top or first level (plev[0]) - - - - - - - - - - First atmospheric layer (play[0]=plev[0], tlay[0], xlay[0]) ------------------- plev[1] - - - - - - - - - - play[1], tlay[1], xlay[1] ------------------- plev[2] - - - - - - - - - - play[2], tlay[2], xlay[2] ------------------- plev[3] - - - - - - - - - - bottom layer (play[Nlay-1]=psurf, tlay[3], xlay[3]) ------------------- Surface (plev[Nlev-1=Nlay]=psurf) /////////////////// .. image:: /images/atm_schema.png Temperatures (`tlay`) and volume mixing ratios (`xlay`) are provided at the *mid layer* point (`play`) (Note that this does not mean that this point is the middle of the layer, in general, it is not). If pressure levels are not specified by the user with `logplevel`, they are at the mid point between the pressure of the layers directly above and below. The pressure of the top and bottom levels are confounded with the pressure of the top and bottom mid layer points. For radiative calculations, the source function (temperature) needs to be known at the boundaries of the *radiative* layers but the opacity needs to be known inside the *radiative* layer. For this reason, there are `Nlay-1` radiative layers and they are offset with respect to atmospheric layers. Opacities are computed at the center of those radiative layers, i.e. at the pressure levels. The temperature is interpolated at these levels with an arithmetic average. Volume mixing ratios are interpolated using a geometric average. Module Contents --------------- .. py:class:: Atm(k_database=None, cia_database=None, a_database=None, wn_range=None, wl_range=None, internal_flux = 0.0, rayleigh = False, flux_top_dw = 0.0, Tstar = 5570, albedo_surf = 0.0, wn_albedo_cutoff=5000.0, **kwargs) Bases: :py:obj:`exo_k.atm_profile.Atm_profile` Class based on Atm_profile that handles radiative transfer calculations. Radiative data are accessed through the :any:`gas_mix.Gas_mix` class. Initialization method that calls Atm_Profile().__init__(**kwargs) and links to Kdatabase and other radiative data. .. py:attribute:: known_spectral_axis :value: False .. py:attribute:: number_computation_albedo_surf_nu :type: int :value: 0 .. py:attribute:: number_computation_flux_top_dw_nu :type: int :value: 0 .. py:attribute:: compute_flux_top_dw_nu :value: True .. py:attribute:: Tstar :type: Optional[float] :value: None .. py:attribute:: stellar_spectrum :type: Optional[exo_k.util.spectrum.Spectrum] :value: None .. py:attribute:: compute_albedo_surf_nu :value: True .. py:attribute:: flux_net_nu :value: None .. py:attribute:: kernel :value: None .. py:method:: set_k_database(k_database = None) Change the radiative database used by the :class:`Gas_mix` object handling opacities inside :class:`Atm`. See :any:`gas_mix.Gas_mix.set_k_database` for details. :param k_database: New Kdatabase to use. :type k_database: :class:`Kdatabase` object .. py:method:: set_cia_database(cia_database = None) Change the CIA database used by the :class:`Gas_mix` object handling opacities inside :class:`Atm`. See :any:`gas_mix.Gas_mix.set_cia_database` for details. :param cia_database: New CIAdatabase to use. :type cia_database: :class:`CIAdatabase` object .. py:method:: set_a_database(a_database = None) Change the Aerosol database used by the :class:`Aerosols` object handling aerosol optical properties. See :any:`aerosols.Aerosols.set_a_database` for details. :param a_database: New Adatabase to use. :type a_database: :class:`Adatabase` object .. py:method:: set_spectral_range(wn_range = None, wl_range = None) Sets the spectral range in which computations will be done by specifying either the wavenumber (in cm^-1) or the wavelength (in micron) range. See :any:`gas_mix.Gas_mix.set_spectral_range` for details. .. py:method:: set_incoming_stellar_flux(flux_top_dw = None, Tstar = None, stellar_spectrum = None, **kwargs) Sets the stellar incoming flux integrated in each wavenumber channel. .. important:: The normalization is such that the flux input is exactly flux_top_dw whatever the spectral range used. :param flux_top_dw: Bolometric Incoming flux (in W/m^2). :param Tstar: Stellar temperature (in K) used to compute the spectral distribution of incoming flux using a blackbody. :param stellar_spectrum: Spectrum of the star in units of per wavenumber. The specific units do not matter as the overall flux will be renormalized. .. py:property:: spectral_incoming_flux Returns a :class:`Spectrum` object with the value of the incoming stellar flux at the top of atmosphere. .. py:method:: set_internal_flux(internal_flux) Sets internal flux from the subsurface in W/m^2 .. py:method:: set_surface_albedo(albedo_surf = None, wn_albedo_cutoff = None, **kwargs) Sets the value of the mean surface albedo using an heaviside filter. :param albedo_surf: Effective visible surface albedo. If `albedo_surf` is an array, we suppose that the albedo is already in the form for `albedo_surf_nu`. :param wn_albedo_cutoff: wavenumber value in cm-1 dividing the visible range from the IR range, where the albedo goes from a given value 'albedo_surf' to 0. .. py:method:: _setup_albedo_surf_nu() Compute the value of the mean surface albedo for each wavenumber in an array. If the albedo value is modified, this method needs to be called again. For now, only available with wavenumbers .. py:method:: set_rayleigh(rayleigh = False) Sets whether or not Rayleigh scattering is included. .. py:method:: spectral_integration(spectral_var) Spectrally integrate an array, taking care of whether we are dealing with corr-k or xsec data. :param spectral_var: array to integrate :returns: var: np.ndarray array integrated over wavenumber (and g-space if relevant) .. py:method:: interpolate_atm_profile(logplay, Mgas=None, gas_vmr=None, **kwargs) Re interpolates the current profile on a new log pressure grid. Recomputes the radiative kernel on the new grid. See Atm_profile.interpolate_profile for details and options :param logplay: New log pressure grid :type logplay: :class:`array` .. py:method:: g_integration(spectral_var) Integrate an array along the g direction (only for corrk) :param spectral_var: array to integrate :returns: var: array, np.ndarray array integrated over g-space if relevant .. py:method:: opacity(rayleigh = None, compute_all_opt_prop = False, wn_range = None, wl_range = None, **kwargs) Computes the opacity of each of the radiative layers (m^2/molecule). :param rayleigh: If true, the rayleigh cross section is computed in self.kdata_scat and added to kdata(total extinction cross section) If None, the global attribute self.rayleigh is used. :param compute_all_opt_prop: :param wn_range: :param wl_range: See :any:`gas_mix.Gas_mix.cross_section` for details. .. py:method:: source_function(integral = True, source = True) Compute the blackbody source function (Pi*Bnu) for each layer of the atmosphere. :param integral: * If true, the black body is integrated within each wavenumber bin. * If not, only the central value is used. False is faster and should be ok for small bins, but True is the correct version. :param source: If False, the source function is put to 0 (for solar absorption calculations) :returns: **$\pi B_\nu$** -- Blackbody source function. :rtype: :class:`np.ndarray` .. py:method:: setup_emission_caculation(mu_eff = 0.5, rayleigh = None, integral = True, source = True, gas_vmr=None, Mgas=None, aer_reffs_densities=None, **kwargs) Computes all necessary quantities for emission calculations (opacity, source, etc.) .. py:method:: emission_spectrum(integral=True, mu0=0.5, mu_quad_order=None, dtau_min=1e-13, **kwargs) Returns the emission flux at the top of the atmosphere (in W/m^2/cm^-1) :param integral: * If true, the black body is integrated within each wavenumber bin. * If not, only the central value is used. False is faster and should be ok for small bins, but True is the correct version. :type integral: :class:`boolean`, *optional* :param mu0: Cosine of the quadrature angle use to compute output flux :type mu0: :class:`float` :param mu_quad_order: If an integer is given, the emission intensity is computed for a number of angles and integrated following a gauss legendre quadrature rule of order `mu_quad_order`. :type mu_quad_order: :class:`int` :param dtau_min: If the optical depth in a layer is smaller than dtau_min, dtau_min is used in that layer instead. Important as too transparent layers can cause important numerical rounding errors. :type dtau_min: :class:`float` :returns: Spectrum object A spectrum with the Spectral flux at the top of the atmosphere (in W/m^2/cm^-1) .. py:method:: emission_spectrum_quad(integral=True, mu_quad_order=3, dtau_min=1e-13, **kwargs) Returns the emission flux at the top of the atmosphere (in W/m^2/cm^-1) using gauss legendre qudrature of order `mu_quad_order` :param integral: * If true, the black body is integrated within each wavenumber bin. * If not, only the central value is used. False is faster and should be ok for small bins, but True is the correct version. :type integral: :class:`boolean`, *optional* :param dtau_min: If the optical depth in a layer is smaller than dtau_min, dtau_min is used in that layer instead. Important as too transparent layers can cause important numerical rounding errors. :type dtau_min: :class:`float` :returns: Spectrum object A spectrum with the Spectral flux at the top of the atmosphere (in W/m^2/cm^-1) .. py:method:: emission_spectrum_2stream(integral = True, mu0 = 0.5, method = 'toon', dtau_min = 1e-10, flux_at_level = False, rayleigh = None, source = True, compute_kernel = False, flux_top_dw = None, mu_star = None, stellar_mode = 'diffusive', **kwargs) Returns the emission flux at the top of the atmosphere (in W/m^2/cm^-1) :param integral: * If True, the source function (black body) is integrated within each wavenumber bin. * If False, only the central value is used. False is faster and should be ok for small bins (e.g. used with `Xtable`), but True is the most accurate mode. :param rayleigh: Whether to account for rayleigh scattering or not. If None, the global attribute self.rayleigh is used. :param mu0: Cosine of the quadrature angle used to compute output flux :param flux_top_dw: (W/m^2) allows you to specify the stellar flux impinging on the top of the atmosphere. (Used to rescale the stellar spectrum). :param stellar_mode: When `flux_top_dw` is provided, set the method to be used to take it into account. - ``diffusive``: Incoming diffuse flux at the upper boundary. - ``collimated``: Incoming diffuse flux is treated as a source term. :param mu_star: $\mu_*$ is the incident direction of the solar beam. Used when the incoming diffuse flux is treated as a source term. :param dtau_min: If the optical depth in a layer is smaller than dtau_min, dtau_min is used in that layer instead. Important as too transparent layers can cause important numerical rounding errors. :param flux_at_level: Whether the fluxes are computed at the levels (True) or in the middle of the layers (False). e.g. this needs to be True to be able to use the kernel with the evolution model. :param source: Whether to include the self emission of the gas in the computation. :param compute_kernel: Whether we want to recompute the kernel of the heating rates along with the fluxes. :param method: What method to use to compute the 2 stream fluxes. - 'toon' (default) uses Toon et al. (1989). :returns: Spectrum object : A spectrum with the Spectral flux at the top of the atmosphere (in W/m^2/cm^-1) .. py:method:: _compute_kernel(solve_2stream_specialized, epsilon = 0.01, per_unit_mass = True, integral = True, **kwargs) Computes the Jacobian matrix d Heating[lay=i] / d T[lay=j] The user should not call this method directly. To recompute the kernel, call `emission_spectrum_2stream` with the option `compute_kernel=True`. After the computation, the kernel is stored in `self.kernel` and the temperature array at which the kernel has been computed in `self.tlay_kernel`. :param solve_2stream_specialized: Name of the function that will be used to compute the fluxes. `emission_spectrum_2stream` will create a lambda function that take care of specifying all arguments leaving `source_nu` to be specified. :param epsilon: The temperature increment used to compute the derivative will be `epsilon*self.tlay` :param per_unit_mass: Whether the heating rates will be normalized per unit of mass (i.e. the kernel values will have units of W/kg/K). If False, kernel in W/layer/K. :param integral: Whether the integral mode is used in flux computations. e.g. this needs to be `True` to be able to use the kernel with the evolution model. .. py:method:: flux_divergence(per_unit_mass = True) Computes the divergence of the net flux in the layers (used to compute heating rates). :func:`emission_spectrum_2stream` needs to be run first. :param per_unit_mass: If True, the heating rates are normalized by the mass of each layer (result in W/kg). :returns: H: np.ndarray Heating rate in each layer (Difference of the net fluxes). Positive means heating. The last value is the net flux impinging on the surface + the internal flux. net: np.ndarray Net fluxes at level surfaces .. py:method:: bolometric_fluxes(per_unit_mass = True) Computes the bolometric fluxes at levels and the divergence of the net flux in the layers (used to compute heating rates). :func:`emission_spectrum_2stream` needs to be run first. :param per_unit_mass: If True, the heating rates are normalized by the mass of each layer (result in W/kg). :returns: up: np.ndarray Upward fluxes at level surfaces dw: np.ndarray Downward fluxes at level surfaces net: np.ndarray Net fluxes at level surfaces H: np.ndarray Heating rate in each layer (Difference of the net fluxes). Positive means heating. The last value is the net flux impinging on the surface + the internal flux. .. py:method:: heating_rate(compute_kernel = False, per_unit_mass = True, dTmax_use_kernel = None, **kwargs) Computes the heating rates and net fluxes in the atmosphere. :param compute_kernel: If True, the kernel (jacobian of the heat rates) is recomputed :param per_unit_mass: If True, the heating rates are normalized by the mass of each layer (result in W/kg). :param dTmax_use_kernel: Maximum temperature difference between the current temperature and the temperature of the last kernel computation before a new kernel is recomputed. :param kwargs: Arguments to use when calling `emission_spectrum_2stream`. :returns: H: np.ndarray Heating rate in each layer (Difference of the net fluxes). Positive means heating. The last value is the net flux impinging on the surface + the internal flux. net: np.ndarray Net fluxes at level surfaces .. py:method:: bolometric_fluxes_2band(flux_top_dw=None, per_unit_mass = True, **kwargs) Computes the bolometric fluxes at levels and heating rates using `bolometric_fluxes`. However, the (up, down, net) fluxes are separated in two contributions: - the part emitted directly by the atmosphere (_emis). - the part due to the incoming stellar light (_stell), that can be used to compute the absorbed stellar radiation and the bond_albedo. We also provide the associated heating rates (H in W/kg) .. py:method:: spectral_fluxes_2band(**kwargs) Computes the spectral fluxes at levels. The (up, down, net) fluxes are separated in two contributions: - the part emitted directly by the atmosphere (_emis). - the part due to the incoming stellar light (_stell), that can be used to compute the absorbed stellar radiation and the bond_albedo. .. py:method:: transmittance_profile(**kwargs) Computes the transmittance profile of an atmosphere, i.e. Exp(-tau) for each layer of the model. Real work done in the numbafied function path_integral_corrk/xsec depending on the type of data. .. py:method:: transmission_spectrum(normalized=False, Rstar=None, **kwargs) Computes the transmission spectrum of the atmosphere. In general (see options below), the code returns the transit depth: .. math:: \delta_\nu=(\pi R_p^2+\alpha_\nu)/(\pi R_{star}^2), where .. math:: \alpha_\nu=2 \pi \int_0^{z_{max}} (R_p+z)*(1-e^{-\tau_\nu(z)) d z. :param Rstar: Radius of the host star. If a float is specified, meters are assumed. Does not need to be given here if it has already been specified as an attribute of the :class:`Atm` object. If specified, the result is the transit depth: .. math:: \delta_\nu=(\pi R_p^2+\alpha_\nu)/(\pi R_{star}^2). :type Rstar: :class:`float` or :class:`astropy.unit object`, *optional* :param normalized: Used only if self.Rstar and Rstar are None: * If True, the result is normalized to the planetary radius: .. math:: \delta_\nu=1+\frac{\alpha_\nu}{\pi R_p^2}. * If False, .. math:: \delta_\nu=\pi R_p^2+\alpha_\nu. :type normalized: :class:`boolean`, *optional* :returns: array The transit spectrum (see above for normalization options). .. py:method:: exp_minus_tau() Sums Exp(-tau) over gauss points .. py:method:: exp_minus_tau_g(g_index) Sums Exp(-tau) over gauss point .. py:method:: blackbody(layer_idx = -1, integral = True) Computes the surface black body flux (in W/m^2/cm^-1) for the temperature of layer `layer_idx`. :param layer_idx: Index of layer used for the temperature. :param integral: * If true, the black body is integrated within each wavenumber bin. * If not, only the central value is used. False is faster and should be ok for small bins, but True is the correct version. :returns: Spectrum object Spectral flux in W/m^2/cm^-1 .. py:method:: surf_bb(**kwargs) Computes the surface black body flux (in W/m^2/cm^-1). See :any:`blackbody` for options. :returns: Spectrum object Spectral flux of a bb at the surface (in W/m^2/cm^-1) .. py:method:: top_bb(**kwargs) Computes the top of atmosphere black body flux (in W/m^2/cm^-1). See :any:`blackbody` for options. :returns: Spectrum object Spectral flux of a bb at the temperature at the top of atmosphere (in W/m^2/cm^-1) .. py:method:: contribution_function() Compute the contribution function $Cf(P,\lambda) = B(T,\lambda) \frac{d e^\tau}{d \log P}$. $\tau$ and $\log P$ are taken at the mid-layers, the temperature needs to be taken at the level surfaces. For that, we use the computed temperature `t_opac`. The result is not normalized to allow comparison. :returns: Contribution function [W/m^2/str/cm-1]. Shape: **(Nlay - 1, Nw)**. :rtype: :class:`np.ndarray`