# Path: helita/utils/utilsmath.py
# def planck(wavelength, temp, dist='wavelength'):
#     """
#     Calculates the Planck function, either per wavelength or per
#     frequency.
# 
#     Parameters
#     ----------
#     wavelength : `Quantity` object (number or sequence, length units)
#         Wavelength(s) to calculate.
#     temp : `Quantity` object (number or sequence, temperature units)
#         Temperature(s) to calculate
#     dist: str, optional
#         How to calculate the distribution. Options are 'wavelength' (default),
#         or 'frequency'.
# 
#     Returns
#     -------
#     planck : `Quantity` object (number or sequence)
#         Planck distribution. Units are energy per time per area
#         per frequency (or wavelength) per solid angle.
# 
#     Notes
#     -----
#     For solid angle integrated one must multiply it by pi.
# 
#     """
#     from astropy.constants import c, h, k_B
#     import astropy.units as u
# 
#     wave = wavelength.to('nm')
#     if temp.shape and wave.shape:
#         temp = temp[:, np.newaxis]  # array broadcast when T and wave are arrays
#     if dist.lower() == 'wavelength':
#         iplanck = 2 * h * c**2 / wave**5 / (np.exp(h * c / (wave * k_B * temp)) - 1)
#         return (iplanck / u.sr).to('erg / (s cm2 Angstrom sr)')
#     elif dist.lower() == 'frequency':
#         iplanck = 2 * h * c / wave**3 / (np.exp(h * c / (wave * k_B * temp)) - 1)
#         return (iplanck / u.sr).to('W / (m2 Hz sr)')
#     else:
#         raise ValueError('invalid distribution ' % dist)
# 
# @vectorize([float32(float32, float32), float64(float64, float64)])
# def voigt(a, v):
#     """
#     Returns the Voigt function:
# 
#     H(a,v) = a/pi * \int_{-Inf}^{+Inf} exp(-y**2)/[(v-y)**2 + a**2] dy
# 
#     Based on approximation from old Fortran routine voigtv from Aake Nordlund.
#     Makes use of numba vectorize, can be used as numpy ufunc.
# 
#     Parameters
#     ----------
#     a : scalar or n-D array (float)
#         Parameter 'a' in Voigt function, typically a scalar. If n-D, must
#         have same shape of v.
#     v : scalar or n-D array (float)
#         Velocity or Doppler value for Voigt function, typically a 1D array.
# 
#     Returns
#     -------
#     h : scalar or n-D array (float)
#         Voigt function. Same shape and type as inputs.
#     """
#     a0 = 122.607931777104326
#     a1 = 214.382388694706425
#     a2 = 181.928533092181549
#     a3 = 93.155580458138441
#     a4 = 30.180142196210589
#     a5 = 5.912626209773153
#     a6 = 0.564189583562615
#     b0 = 122.607931773875350
#     b1 = 352.730625110963558
#     b2 = 457.334478783897737
#     b3 = 348.703917719495792
#     b4 = 170.354001821091472
#     b5 = 53.992906912940207
#     b6 = 10.479857114260399
#     if a == 0:
#         return exp(-v ** 2)
#     z = v * 1j + a
#     h = (((((((a6 * z + a5) * z + a4) * z + a3) * z + a2) * z + a1) * z + a0) /
#      (((((((z + b6) * z + b5) * z + b4) * z + b3) * z + b2) * z + b1) * z + b0))
#     return h.real


import os
import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
from pkg_resources import resource_filename
from ipywidgets import interact, fixed, Dropdown, IntSlider, FloatSlider
from scipy.integrate import cumtrapz
from scipy.interpolate import interp1d
from astropy import units as u
from ..utils.utilsmath import planck, voigt
                             dist='frequency').value
            fig, ax = plt.subplots()
            ax.plot(height, sf, 'b-', label=r'S$_\mathrm{total}$', lw=1)
            ax.set_yscale('log')
            ax.plot(height, obj.ray.Jlambda[0, 0, :, 0].dropna('height'),
                    'y-', label='J', lw=1)
            ax.plot(height, bplanck, 'r--', label=r'B$_\mathrm{Planck}$',
                    lw=1)
            ax.set_xlabel("Height (Mm)")
            ax.set_ylabel(r"W m$^{-2}$ Hz$^{-1}$ sr$^{-1}$")
            ax.set_title("%.3f nm" % obj.ray.wavelength_selected[0])
            lg = ax.legend(loc='upper center')
            lg.draw_frame(0)
            # tau annotations
            tau_v, h_idx = __get_tau_levels(0, 0, 0)
            for i, level in enumerate(tau_levels):
                xval = tau_v[i]
                yval = sf[h_idx[i]]
                ax.annotate(r'$\tau$=%s' % level,
                            xy=(xval, yval),
                            xytext=(xval, yval / (0.2 - 0.03 * i)),
                            arrowprops=ARROW, ha='center', va='top')
            return ax

        ax = _sf_plot()

        @interact(wavelength=(0, nwave - 1, 1), y_log=True,
                  x=x_slider, y=y_slider)
        def _sf_update(wavelength=0, y_log=True, x=0, y=0):
            obj = self.rhobj

# Path: helita/utils/utilsmath.py
# def planck(wavelength, temp, dist='wavelength'):
#     """
#     Calculates the Planck function, either per wavelength or per
#     frequency.
# 
#     Parameters
#     ----------
#     wavelength : `Quantity` object (number or sequence, length units)
#         Wavelength(s) to calculate.
#     temp : `Quantity` object (number or sequence, temperature units)
#         Temperature(s) to calculate
#     dist: str, optional
#         How to calculate the distribution. Options are 'wavelength' (default),
#         or 'frequency'.
# 
#     Returns
#     -------
#     planck : `Quantity` object (number or sequence)
#         Planck distribution. Units are energy per time per area
#         per frequency (or wavelength) per solid angle.
# 
#     Notes
#     -----
#     For solid angle integrated one must multiply it by pi.
# 
#     """
#     from astropy.constants import c, h, k_B
#     import astropy.units as u
# 
#     wave = wavelength.to('nm')
#     if temp.shape and wave.shape:
#         temp = temp[:, np.newaxis]  # array broadcast when T and wave are arrays
#     if dist.lower() == 'wavelength':
#         iplanck = 2 * h * c**2 / wave**5 / (np.exp(h * c / (wave * k_B * temp)) - 1)
#         return (iplanck / u.sr).to('erg / (s cm2 Angstrom sr)')
#     elif dist.lower() == 'frequency':
#         iplanck = 2 * h * c / wave**3 / (np.exp(h * c / (wave * k_B * temp)) - 1)
#         return (iplanck / u.sr).to('W / (m2 Hz sr)')
#     else:
#         raise ValueError('invalid distribution ' % dist)
# 
# @vectorize([float32(float32, float32), float64(float64, float64)])
# def voigt(a, v):
#     """
#     Returns the Voigt function:
# 
#     H(a,v) = a/pi * \int_{-Inf}^{+Inf} exp(-y**2)/[(v-y)**2 + a**2] dy
# 
#     Based on approximation from old Fortran routine voigtv from Aake Nordlund.
#     Makes use of numba vectorize, can be used as numpy ufunc.
# 
#     Parameters
#     ----------
#     a : scalar or n-D array (float)
#         Parameter 'a' in Voigt function, typically a scalar. If n-D, must
#         have same shape of v.
#     v : scalar or n-D array (float)
#         Velocity or Doppler value for Voigt function, typically a 1D array.
# 
#     Returns
#     -------
#     h : scalar or n-D array (float)
#         Voigt function. Same shape and type as inputs.
#     """
#     a0 = 122.607931777104326
#     a1 = 214.382388694706425
#     a2 = 181.928533092181549
#     a3 = 93.155580458138441
#     a4 = 30.180142196210589
#     a5 = 5.912626209773153
#     a6 = 0.564189583562615
#     b0 = 122.607931773875350
#     b1 = 352.730625110963558
#     b2 = 457.334478783897737
#     b3 = 348.703917719495792
#     b4 = 170.354001821091472
#     b5 = 53.992906912940207
#     b6 = 10.479857114260399
#     if a == 0:
#         return exp(-v ** 2)
#     z = v * 1j + a
#     h = (((((((a6 * z + a5) * z + a4) * z + a3) * z + a2) * z + a1) * z + a0) /
#      (((((((z + b6) * z + b5) * z + b4) * z + b3) * z + b2) * z + b1) * z + b0))
#     return h.real


# Path: helita/vis/rh15d_vis.py
import os
import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
from pkg_resources import resource_filename
from ipywidgets import interact, fixed, Dropdown, IntSlider, FloatSlider
from scipy.integrate import cumtrapz
from scipy.interpolate import interp1d
from astropy import units as u
from ..utils.utilsmath import planck, voigt

                             dist='frequency').value
            fig, ax = plt.subplots()
            ax.plot(height, sf, 'b-', label=r'S$_\mathrm{total}$', lw=1)
            ax.set_yscale('log')
            ax.plot(height, obj.ray.Jlambda[0, 0, :, 0].dropna('height'),
                    'y-', label='J', lw=1)
            ax.plot(height, bplanck, 'r--', label=r'B$_\mathrm{Planck}$',
                    lw=1)
            ax.set_xlabel("Height (Mm)")
            ax.set_ylabel(r"W m$^{-2}$ Hz$^{-1}$ sr$^{-1}$")
            ax.set_title("%.3f nm" % obj.ray.wavelength_selected[0])
            lg = ax.legend(loc='upper center')
            lg.draw_frame(0)
            # tau annotations
            tau_v, h_idx = __get_tau_levels(0, 0, 0)
            for i, level in enumerate(tau_levels):
                xval = tau_v[i]
                yval = sf[h_idx[i]]
                ax.annotate(r'$\tau$=%s' % level,
                            xy=(xval, yval),
                            xytext=(xval, yval / (0.2 - 0.03 * i)),
                            arrowprops=ARROW, ha='center', va='top')
            return ax

        ax = _sf_plot()

        @interact(wavelength=(0, nwave - 1, 1), y_log=True,
                  x=x_slider, y=y_slider)
        def _sf_update(wavelength=0, y_log=True, x=0, y=0):
            obj = self.rhobj

            bplanck = planck(obj.ray.wavelength_selected[wavelength].values * u.nm,