How to use extinction - 10 common examples

wavelength_avail = {'u':3571,'b':4378,'v':5466,'r':6695,'i':8565,'j':12101,
				   'h':16300,'k':21900,'U':3571,'B':4378,'V':5466,'R':6695,'I':8565,'J':12101,
				   'H':16300,'K':21900} #in Angstrom
		# effective wavelength taken from http://svo2.cab.inta-csic.es/svo/theory/fps3/index.php?mode=browse&gname=Generic&gname2=Johnson
		# and from https://en.wikipedia.org/wiki/Photometric_system for h, k
		
		lr,lg,lb = wavelength_avail[r_band],wavelength_avail[g_band],wavelength_avail[b_band] #in Angstrom
		wave = np.array([lb, lg, lr])

		ext_r = np.empty_like(r)
		ext_g = np.empty_like(g)
		ext_b = np.empty_like(b)
	
		for i in range(len(a_v)):
			for j in range(len(a_v[0])):
				ext = extinction.calzetti00(wave.astype(np.float), a_v[i][j].astype(np.float), 3.1, unit='aa', out=None)
				ext_r[i][j] = ext[2]
				ext_g[i][j] = ext[1]
				ext_b[i][j] = ext[0]

		r = r+ext_r
		g = g+ext_g
		b = b+ext_b

	#r,g,b = nw_scale_rgb(r,g,b)
	#r,g,b = nw_arcsinh_fit(r,g,b)

	if mag_range is None:
		rgbim, mag_max = combine(r, g, b, dynamic_range*2.5)
		mag_min = mag_max + 2.5*dynamic_range
	else:
		mag_max, mag_min = mag_range

source_suffix = '.rst'


# The master toctree document.
master_doc = 'index'

# General information about the project.
project = u'extinction'
copyright = u'2016, Kyle Barbary and contributors'

# The version info for the project you're documenting, acts as replacement for
# |version| and |release|, also used in various other places throughout the
# built documents.
#
# THe short X.Y version.
version = '.'.join(extinction.__version__.split('.')[0:2])

# The full version, including alpha/beta/rc tags.
release = extinction.__version__

# List of patterns, relative to source directory, that match files and
# directories to ignore when looking for source files.
exclude_patterns = ['_build']

# The reST default role (used for this markup: `text`) to use for all
# documents.
default_role = 'obj'

# The name of the Pygments (syntax highlighting) style to use.
pygments_style = 'sphinx'

# The master toctree document.
master_doc = 'index'

# General information about the project.
project = u'extinction'
copyright = u'2016, Kyle Barbary and contributors'

# The version info for the project you're documenting, acts as replacement for
# |version| and |release|, also used in various other places throughout the
# built documents.
#
# THe short X.Y version.
version = '.'.join(extinction.__version__.split('.')[0:2])

# The full version, including alpha/beta/rc tags.
release = extinction.__version__

# List of patterns, relative to source directory, that match files and
# directories to ignore when looking for source files.
exclude_patterns = ['_build']

# The reST default role (used for this markup: `text`) to use for all
# documents.
default_role = 'obj'

# The name of the Pygments (syntax highlighting) style to use.
pygments_style = 'sphinx'


# -- Options for HTML output ----------------------------------------------

# The theme to use for HTML and HTML Help pages.  See the documentation for

if self.mwebv:
            extinctions = extinction.fitzpatrick99(wave=self._good_filter_wave, \
                                                   a_v=-3.1 * self.mwebv, r_v=3.1, unit='aa')

        for i, pb in enumerate(self._good_filters):
            mask = (self.passband == pb)

            flux_pb = self.flux[mask]
            fluxerr_pb = self.fluxErr[mask]
            npbobs = len(flux_pb)

            if npbobs < 1:
                return

            if self.mwebv:
                flux_out = extinction.apply(extinctions[i], flux_pb, inplace=False)
                fluxerr_out = extinction.apply(extinctions[i], fluxerr_pb, inplace=False)
            else:
                flux_out = flux_pb
                fluxerr_out = fluxerr_pb
                self.fluxUnred[mask] = flux_out
                self.fluxErrUnred[mask] = fluxerr_out

            if npbobs > 1:
                # there's at least enough observations to find minimum and maximum
                minfluxpb = flux_out.min()
                maxfluxpb = flux_out.max()
                norm = maxfluxpb - minfluxpb
                self.fluxRenorm[mask] = (flux_out - minfluxpb) / norm
                self.fluxErrRenorm[mask] = fluxerr_out / norm
            elif npbobs == 1:
                # deal with the case with one observation in this passband by setting renorm = 0.5

extinctions = extinction.fitzpatrick99(wave=self._good_filter_wave, \
                                                   a_v=-3.1 * self.mwebv, r_v=3.1, unit='aa')

        for i, pb in enumerate(self._good_filters):
            mask = (self.passband == pb)

            flux_pb = self.flux[mask]
            fluxerr_pb = self.fluxErr[mask]
            npbobs = len(flux_pb)

            if npbobs < 1:
                return

            if self.mwebv:
                flux_out = extinction.apply(extinctions[i], flux_pb, inplace=False)
                fluxerr_out = extinction.apply(extinctions[i], fluxerr_pb, inplace=False)
            else:
                flux_out = flux_pb
                fluxerr_out = fluxerr_pb
                self.fluxUnred[mask] = flux_out
                self.fluxErrUnred[mask] = fluxerr_out

            if npbobs > 1:
                # there's at least enough observations to find minimum and maximum
                minfluxpb = flux_out.min()
                maxfluxpb = flux_out.max()
                norm = maxfluxpb - minfluxpb
                self.fluxRenorm[mask] = (flux_out - minfluxpb) / norm
                self.fluxErrRenorm[mask] = fluxerr_out / norm
            elif npbobs == 1:
                # deal with the case with one observation in this passband by setting renorm = 0.5
                norm = self.fluxUnred[mask] / 0.5

def transform(self, wave, flux):
        if self.Rv is not None:
            extinct_mag = self.law(wave, self.Av, self.Rv)
        else:
            extinct_mag = self.law(wave, self.Av)
        extinct_flux = extinction.apply(extinct_mag, flux)
        return wave, extinct_flux

if bi >= 0:
                if bi not in extinct_cache:
                    extinct_cache[bi] = np.zeros_like(
                        self._band_rest_wavelengths[bi])
                    ind = self._ext_indices[bi]
                    if len(ind) > 0:
                        extinct_cache[bi][ind] = odonnell94(
                            self._band_rest_wavelengths[bi][ind],
                            av_host, self._rv_host)
                    ind = self._x_indices[bi]
                    if len(ind) > 0:
                        extinct_cache[bi][ind] = self.mm83(
                            self._nh_host,
                            self._band_rest_wavelengths[bi][ind])
                # Add host and MW contributions
                eapp(
                    self._mw_extinct[bi] + extinct_cache[bi],
                    self._seds[si], inplace=True)
            else:
                # wavelengths = np.array(
                #   [c.c.cgs.value / self._frequencies[si]])
                # Need extinction function for radio
                pass

        # Units of `seds` is ergs / s / Angstrom.
        return {
            'sample_wavelengths': self._sample_wavelengths,
            self.key('seds'): self._seds,
            self.key('avhost'): av_host
        }

:param law: the extinction law, one of {'ccm89', 'odonnell94', 'calzetti00',
        'fitzpatrick99', 'fm07'}
    :type law: str
    :param Av: The scaling total extinction value.
    :type Av: float
    :param Rv: The ratio of total to selective extinction. If using law 'fm07' you do
        not need to provide this (fixed 3.1).
    :type Rv: float

    :raises ValueError: If not using an expected law or ill-specifying Av or Rv
    """
    LAWS = {
        'ccm89': extinction.ccm89,
        'odonnell94': extinction.odonnell94,
        'calzetti00': extinction.calzetti00,
        'fitzpatrick99': extinction.fitzpatrick99,
        'fm07': extinction.fm07,
    }
    def __init__(self, law, Av, Rv=None):
        if not law in self.LAWS:
            raise ValueError('Need to specify a law from {}'.format(self.LAWS.keys()))
        if Av < 0:
            raise ValueError('Cannot have negative extinction')
        if Rv is None or Rv < 0 and law is not 'fm07':
            raise ValueError('Must provide positive r_v for law "{}"'.format(law))
        elif law is 'fm07':
            Rv = None
        self.law = self.LAWS[law]
        self.Av = Av
        self.Rv = Rv

    def transform(self, wave, flux):

def _remove_flux_extinction(self):
        """
        Remove extinction for light curve assuming Fitzpatrick '99 reddening
        law, given some value of E(B-V)
        """
        self.fluxUnred = self.flux.copy()
        self.fluxErrUnred = self.fluxErr.copy()
        self.fluxRenorm = self.flux.copy()
        self.fluxErrRenorm = self.fluxErr.copy()

        # Using negative a_v so that extinction.apply works in reverse and removes the extinction
        if self.mwebv:
            extinctions = extinction.fitzpatrick99(wave=self._good_filter_wave, \
                                                   a_v=-3.1 * self.mwebv, r_v=3.1, unit='aa')

        for i, pb in enumerate(self._good_filters):
            mask = (self.passband == pb)

            flux_pb = self.flux[mask]
            fluxerr_pb = self.fluxErr[mask]
            npbobs = len(flux_pb)

            if npbobs < 1:
                return

            if self.mwebv:
                flux_out = extinction.apply(extinctions[i], flux_pb, inplace=False)
                fluxerr_out = extinction.apply(extinctions[i], fluxerr_pb, inplace=False)
            else:

self._frequencies = kwargs['all_frequencies']
        self._band_rest_wavelengths = self._sample_wavelengths / zp1

        av_host = self._nh_host / 1.8e21

        extinct_cache = OrderedDict()
        for si, cur_band in enumerate(self._bands):
            bi = self._band_indices[si]
            # Extinct out host gal (using rest wavelengths)
            if bi >= 0:
                if bi not in extinct_cache:
                    extinct_cache[bi] = np.zeros_like(
                        self._band_rest_wavelengths[bi])
                    ind = self._ext_indices[bi]
                    if len(ind) > 0:
                        extinct_cache[bi][ind] = odonnell94(
                            self._band_rest_wavelengths[bi][ind],
                            av_host, self._rv_host)
                    ind = self._x_indices[bi]
                    if len(ind) > 0:
                        extinct_cache[bi][ind] = self.mm83(
                            self._nh_host,
                            self._band_rest_wavelengths[bi][ind])
                # Add host and MW contributions
                eapp(
                    self._mw_extinct[bi] + extinct_cache[bi],
                    self._seds[si], inplace=True)
            else:
                # wavelengths = np.array(
                #   [c.c.cgs.value / self._frequencies[si]])
                # Need extinction function for radio
                pass