Cross-correlate the observed spectrum against template fluxes.
772271
def cross_correlate(observed, template_dispersion, template_fluxes,
    rebin="template", wavelength_range=None, continuum_degree=-1,
    apodize=0.10, rescale=False, z_limits=None, full_output=False):
    

    template_fluxes = np.atleast_2d(template_fluxes)
    template_dispersion = np.array(template_dispersion)
    if template_dispersion.shape[0] != template_fluxes.shape[1]:
        raise ValueError("template dispersion must have size (N_pixels,) and "\
            "template fluxes must have size (N_models, N_pixels)")

    try:
        continuum_degree = int(continuum_degree)
    except (TypeError, ValueError):
        raise TypeError("continuum order must be an integer-like object")

    if not (1 > apodize >= 0):
        raise ValueError("apodize fraction must be between 0 and 1")

    if rebin.lower() == "template":
        # Put the template fluxes onto the observed dispersion map.
        dispersion = observed.disp
        template_flux = template_fluxes \
            * resample(template_dispersion, observed.disp)
        observed_flux = observed.flux.copy()
        observed_ivar = observed.ivariance.copy()

    elif rebin.lower() == "observed":
        # Put the observed fluxes onto the template dispersion map.
        dispersion = template_dispersion
        template_flux = template_fluxes
        mat = resample(observed.disp, template_dispersion)
        observed_flux = observed.flux.copy() * mat
        observed_ivar = observed.ivariance.copy() * mat
        
    else:
        raise ValueError("rebin must be either `template` or `observed`")

    if wavelength_range is not None:
        if not isinstance(wavelength_range, (tuple, list, np.ndarray)) \
        or len(wavelength_range) != 2:
            raise TypeError("wavelength range must be a two length tuple")

        indices = np.clip(
            template_dispersion.searchsorted(wavelength_range) + [0, 1],
            0, template_dispersion.size)
        template_dispersion = template_dispersion[indices[0]:indices[1]]
        template_flux = template_flux[indices[0]:indices[1]]
        observed_flux = observed_flux[indices[0]:indices[1]]

    N_templates = template_fluxes.shape[0]

    # Clip out the non-finite edges, if they exist.
    _ = np.where(np.isfinite(template_flux[0] * observed_flux))[0]
    l_idx, u_idx = _[0], _[-1]

    dispersion = dispersion[l_idx:u_idx]
    observed_flux = observed_flux[l_idx:u_idx]
    observed_ivar = observed_ivar[l_idx:u_idx]
    template_flux = template_flux[:, l_idx:u_idx]

    # Ensure an even number of points.
    N = u_idx - l_idx
    N = N - 1 if N % 2 > 0 else N
    dispersion = dispersion[:N]
    observed_flux = observed_flux[:N]
    observed_ivar = observed_ivar[:N]
    template_flux = template_flux[:, :N]

    finite = np.isfinite(observed_flux * observed_ivar) * (observed_flux > 1e-3)
    
    # Continuum.
    if continuum_degree >= 0:
        coefficients = np.polyfit(dispersion[finite], observed_flux[finite],
            continuum_degree)
        observed_flux /= np.polyval(coefficients, dispersion)

    # Interpolate over non-finite pixels.
    observed_flux[~finite] = np.interp(dispersion[~finite], dispersion[finite],
        observed_flux[finite], left=1, right=1)

    # Scale the flux level to that the template intensities
    if rescale:
        observed_flux = (observed_flux * template_flux.ptp()) \
            + template_flux.min()

    # Apodize edges.
    edge_buffer = apodize * (dispersion[-1] - dispersion[0])
    low_w_indices = np.nonzero(dispersion < dispersion[0] + edge_buffer)[0]
    high_w_indices = np.nonzero(dispersion > dispersion[-1] - edge_buffer)[0]

    apod_curve = np.ones(N, dtype='d')
    apod_curve[low_w_indices] = (1.0 + np.cos(np.pi*(
        1.0 - (dispersion[low_w_indices] - dispersion[0])/edge_buffer)))/2.
    apod_curve[high_w_indices] = (1.0 + np.cos(np.pi*(
        1.0 - (dispersion[-1] - dispersion[high_w_indices])/edge_buffer)))/2.

    apod_observed_flux = observed_flux * apod_curve
    apod_template_flux = template_flux * apod_curve

    fft_observed_flux = np.fft.fft(apod_observed_flux)
    fft_template_flux = np.fft.fft(apod_template_flux)
    template_flux_corr = (fft_observed_flux * fft_template_flux.conjugate())
    template_flux_corr /= np.sqrt(
        np.inner(apod_observed_flux, apod_observed_flux))

    z_array = np.array(dispersion.copy())/dispersion[N/2] - 1.0

    z = np.ones(N_templates) * np.nan
    z_err = np.ones(N_templates) * np.nan
    R = np.ones(N_templates) * np.nan
    for i in range(N_templates):

        denominator = np.sqrt(
            np.inner(apod_template_flux[i, :], apod_template_flux[i, :]))
        flux_correlation = template_flux_corr[i, :] / denominator 
        correlation = np.fft.ifft(flux_correlation).real

        # Reflect about zero
        ccf = np.zeros(N)
        ccf[:N/2] = correlation[N/2:]
        ccf[N/2:] = correlation[:N/2]

        # Get height and redshift of best peak
        h = ccf.max()

        # Scale the CCF
        #ccf -= ccf.min()
        #ccf *= (h/ccf.max())
        
        R[i] = h
        z[i] = z_array[ccf.argmax()]
        try:
            z_err[i] = (np.ptp(z_array[np.where(ccf >= 0.5*h)])/2.35482)**2
        except ValueError:
            continue

    # Limits:
    if z_limits is not None:
        R[~((z_limits[1] >= z) * (z >= z_limits[0]))] = np.nan

    # Re-measure the velocity at the best peak.
    index = np.nanargmax(R)

    denominator = np.sqrt(
        np.inner(apod_template_flux[index, :], apod_template_flux[index, :]))

    flux_correlation = template_flux_corr[index, :] / denominator
    correlation = np.fft.ifft(flux_correlation).real

    # FReflect about zero.
    ccf = np.zeros(N)
    ccf[:N/2] = correlation[N/2:]
    ccf[N/2:] = correlation[:N/2]

    h = ccf.max()
    ccf -= ccf.min()
    ccf *= h/ccf.max()

    # Fit +/- 5 pixels
    idx = np.argmax(ccf) - 3, np.argmax(ccf) + 3

    coeffs = np.polyfit(z_array[idx[0]:idx[1]], ccf[idx[0]:idx[1]], 2)
    x_i = np.linspace(z_array[idx[0]], z_array[idx[1]], 1000)
    y_i = np.polyval(coeffs, x_i)

    # Update the value
    z[index] = x_i[y_i.argmax()]

    """
    # Fit the profile peak.
    fwhm = np.ptp(z_array[np.where(ccf >= 0.5)])
    p0 = np.array([z[index], fwhm/2.355, 1.0]) # mu, sigma, peak.
    use = (z_array > (p0[0] - 3 * p0[1])) * ((p0[0] + 3*p0[1]) > z_array)
    x, y = z_array[use], ccf[use]

    # Fit profile:
    f = lambda x, mu, sigma, peak: peak * np.exp(-(mu - x)**2/(2*sigma**2)) 
    
    import matplotlib.pyplot as plt
    plt.plot(z_array, ccf, c='k')
    plt.plot(x, y, c='b')
    plt.plot(x, f(x, *p0), c='r')

    p1 = curve_fit(f, x, y, p0=p0)
    plt.plot(x, f(x, *p1[0]), c='m')

    plt.gca().set_ylim(0,1)
    raise a

    p1 = curve_fit(f, x, y, p0=p0)
    """

    return (z * c, z_err * c, R)
974407
def _check(self):
        """"""
        # Run the super method
        super(Spectrum, self)._check()

        err_str = None
        has_data = self._KEYS.DATA in self
        has_wave = self._KEYS.WAVELENGTHS in self
        has_flux = self._KEYS.FLUXES in self
        has_filename = self._KEYS.FILENAME in self

        if not has_data:
            if (not has_wave or not has_flux) and not has_filename:
                err_str = (
                    "If `{}` not given".format(self._KEYS.DATA) +
                    "; `{}` or `{}` needed".format(
                        self._KEYS.WAVELENGTHS, self._KEYS.FLUXES))

        if err_str is not None:
            raise ValueError(err_str)

        return
437081
def get_spatially_integrated_flux( self, energies):
    
        
       
        if not isinstance(energies, np.ndarray):
            energies = np.array(energies, ndmin=1)

        # Get the differential flux from the spectral components

        results = [self.spatial_shape.get_total_spatial_integral(energies) * component.shape(energies) for component in self.components.values()]

        if isinstance(energies, u.Quantity):

            # Slow version with units

            # We need to sum like this (slower) because using np.sum will not preserve the units
            # (thanks astropy.units)

            differential_flux = sum(results)

        else:

            # Fast version without units, where x is supposed to be in the same units as currently defined in
            # units.get_units()

            differential_flux = np.sum(results, 0)

        return differential_flux
1158644
def get_line_flux(line_wave, wave, flux, **kwargs):
    """"""
    return np.interp(line_wave, wave, flux, **kwargs)
12313
def generate_gaussian_profile(seeing_fwhm):
    """"""
    FWHM_G = 2 * math.sqrt(2 * math.log(2))
    sigma = seeing_fwhm / FWHM_G
    amplitude = 1.0 / (2 * math.pi * sigma * sigma)
    seeing_model = Gaussian2D(amplitude=amplitude,
                              x_mean=0.0,
                              y_mean=0.0,
                              x_stddev=sigma,
                              y_stddev=sigma)
    return seeing_model
856904
def estimate(self, data, full_output=False, **kwargs):
        

        # Number of model comparisons can be specified in the configuration.
        num_model_comparisons = self._configuration.get("estimate", {}).get(
            "num_model_comparisons", self.grid_points.size)
        # If it's a fraction, we need to convert that to an integer.
        if 1 > num_model_comparisons > 0:
            num_model_comparisons *= self.grid_points.size

        # If the num_model_comparison is provided as a keyword argument, use it.
        num_model_comparisons = kwargs.pop("num_model_comparisons",
            int(num_model_comparisons))

        logger.debug("Number of model comparisons to make for initial estimate:"
            " {0}".format(num_model_comparisons))
        
        # Match the data to the model channels.
        matched_channels, missing_channels, ignore_parameters \
            = self._match_channels_to_data(data)

        logger.debug("Matched channels: {0}, missing channels: {1}, ignore "
            "parameters: {2}".format(matched_channels, missing_channels,
                ignore_parameters))

        # Load the intensities
        t = time()
        s = self.grid_points.size/num_model_comparisons # step size
        grid_points = self.grid_points[::s]
        intensities = np.memmap(
            self._configuration["model_grid"]["intensities"], dtype="float32",
            mode="r", shape=(self.grid_points.size, self.wavelengths.size))[::s]
        logger.debug("Took {:.0f} seconds to load and slice intensities".format(
            time() - t))
        # Which matched, data channel has the highest S/N?
        # (This channel will be used to estimate astrophysical parameters)
        data, pixels_affected = self._apply_data_mask(data)
        median_snr = dict(zip(matched_channels,
            [np.nanmedian(spec.flux/(spec.variance**0.5)) for spec in data]))
        median_snr.pop(None, None) # Remove unmatched data spectra

        ccf_channel = self._configuration.get("settings", {}).get("ccf_channel",
            max(median_snr, key=median_snr.get))
        if ccf_channel not in matched_channels:
            logger.warn("Ignoring CCF channel {0} because it was not a matched"
                " channel".format(ccf_channel))
            ccf_channel = max(median_snr, key=median_snr.get)

        logger.debug("Channel with peak SNR is {0}".format(ccf_channel))

        # Are there *any* continuum parameters in any matched channel?
        any_continuum_parameters = any(map(lambda s: s.startswith("continuum_"),
            set(self.parameters).difference(ignore_parameters)))

        # [TODO]: CCF MASK
        # [TODO]: Don't require CCF if we have only continuum parameters.

        z_limits = self._configuration["settings"].get("ccf_z_limits", None)

        theta = {} # Dictionary for the estimated model parameters.
        best_grid_index = None
        c = speed_of_light.to("km/s").value
        for matched_channel, spectrum in zip(matched_channels, data):
            if matched_channel is None: continue

            # Do we need todo cross-correlation for this channel?
            # We do if there are redshift parameters for this channel,
            # or if there is a global redshift or global continuum parameters
            #   and this channel is the highest S/N.
            if "z_{}".format(matched_channel) in self.parameters \
            or ((any_continuum_parameters or "z" in self.parameters) \
            and matched_channel == ccf_channel):

                # Get the continuum degree for this channel.
                continuum_degree = self._configuration["model"].get("continuum",
                    { matched_channel: -1 })[matched_channel]

                logger.debug("Perfoming CCF on {0} channel with a continuum "
                    "degree of {1}".format(matched_channel, continuum_degree))

                # Get model wavelength indices that match the data.
                # get the points that are in the mask, and within the spectrum
                # limits

                # TODO: Make this CCF not model mask.
                idx = np.where(self._model_mask() \
                    * (self.wavelengths >= spectrum.disp[0]) \
                    * (spectrum.disp[-1] >= self.wavelengths))[0]

                v, v_err, R = spectrum.cross_correlate(
                    (self.wavelengths[idx], intensities[:, idx]),
                    #(self.wavelengths, intensities),
                    continuum_degree=continuum_degree, z_limits=z_limits)

                # Identify the best point by the CCF peak.
                best = np.nanargmax(R)

                # Now, why did we do CCF in this channel? Which model parameters
                # should be updated?
                if "z_{}".format(matched_channel) in self.parameters:
                    theta["z_{}".format(matched_channel)] = v[best] / c
                elif "z" in self.parameters: 
                    # If there is a global redshift, update it.
                    theta["z"] = v[best] / c

                    # Continuum parameters will be updated later, so that each
                    # channel is checked to see if it has the highest S/N,
                    # otherwise we might be trying to calculate continuum
                    # parameters when we haven't done CCF on the highest S/N
                    # spectra yet.

                if matched_channel == ccf_channel:
                    # Update astrophysical parameters.
                    theta.update(dict(zip(grid_points.dtype.names,
                        grid_points[best])))
                    best_grid_index = best

        # If there are continuum parameters, calculate them from the best point.
        if any_continuum_parameters:
            for matched_channel, spectrum in zip(matched_channels, data):
                if matched_channel is None: continue

                # The template spectra at the best point needs to be
                # redshifted to the data, and then continuum coefficients
                # calculated from that.

                # Get the continuum degree for this channel.
                continuum_degree = self._configuration["model"].get("continuum",
                    { matched_channel: -1 })[matched_channel]

            
                # Get model wavelength indices that match the data.
                idx = np.clip(self.wavelengths.searchsorted(
                    [spectrum.disp[0], spectrum.disp[-1]]) + [0, 1],
                    0, self.wavelengths.size)

                # Redshift and bin the spectrum.
                z = theta.get("z_{}".format(matched_channel), theta.get("z", 0))

                best_intensities \
                    = np.copy(intensities[best_grid_index, idx[0]:idx[1]]).flatten()

                # Apply model mask.
                model_mask = self._model_mask(self.wavelengths[idx[0]:idx[1]])
                best_intensities[~model_mask] = np.nan

                best_intensities = best_intensities * specutils.sample.resample(
                    self.wavelengths[idx[0]:idx[1]] * (1 + z), spectrum.disp)
                
                # Calculate the continuum coefficients for this channel.
                continuum = spectrum.flux/best_intensities
                finite = np.isfinite(continuum)

                try:
                    coefficients = np.polyfit(
                        spectrum.disp[finite], continuum[finite], continuum_degree,
                        )#w=spectrum.ivariance[finite])
                except np.linalg.linalg.LinAlgError:
                    logger.exception("Exception in initial polynomial fit")
                    coefficients = np.polyfit(spectrum.disp[finite], continuum[finite],
                        continuum_degree)

                # They go into theta backwards. such that coefficients[-1] is
                # continuum_{name}_0
                theta.update(dict(zip(
                    ["continuum_{0}_{1}".format(matched_channel, i) \
                        for i in range(continuum_degree + 1)],
                    coefficients[::-1]
                )))

        # Remaining parameters could be: resolving power, outlier pixels,
        # underestimated variance.
        remaining_parameters = set(self.parameters)\
            .difference(ignore_parameters)\
            .difference(theta)

        if remaining_parameters:
            logger.debug("Remaining parameters to estimate: {0}. For these we "
                "will just assume reasonable initial values.".format(
                remaining_parameters))

            for parameter in remaining_parameters:
                if parameter == "resolution" \
                or parameter.startswith("resolution_"):

                    if parameter.startswith("resolution_"):
                        spectra = [data[matched_channels.index(
                            parameter.split("_")[1])]]
                    else:
                        spectra = [s for s in data if s is not None]

                    R = [s.disp.mean()/np.diff(s.disp).mean() for s in spectra]

                    # Assume oversampling rate of ~5.
                    theta.update({ parameter: np.median(R)/5.})

                elif parameter == "ln_f" or parameter.startswith("ln_f_"):
                    theta.update({ parameter: 0.5 }) # Not overestimated.

                elif parameter in ("Po", "Vo"):
                    theta.update({
                        "Po": 0.01, # 1% outlier pixels.
                        "Vo": np.mean([np.nanmedian(s.variance) for s in data]),
                    })

        logger.info("Initial estimate: {}".format(theta))
        # Having full_output = True means return the best spectra estimate.
        if full_output:

            # Create model fluxes and calculate some metric.
            __intensities = np.copy(intensities[best_grid_index])

            # Apply model masks.
            __intensities[~self._model_mask()] = np.nan

            chi_sq, dof, model_fluxes = self._chi_sq(theta, data,
                __intensities=__intensities, __no_precomputed_binning=True)
            del intensities

            return (theta, chi_sq, dof, model_fluxes)

        # Delete the reference to intensities
        del intensities
        return theta
662168
def _format_data(self, data):
        
        return [spectrum for spectrum in \
            sorted(data if isinstance(data, (list, tuple)) else [data],
                key=lambda x: x.disp[0]) if np.any(np.isfinite(spectrum.flux))]
790723
def estimate_noise(fluxes, contmask):
    
    nstars = fluxes.shape[0]
    scatter = np.zeros(nstars)
    for i,spec in enumerate(fluxes): 
        cont = spec[contmask]
        scatter[i] = stats.funcs.mad_std(cont)
    return scatter
1250196
def _vis_calibrate(self, data, key):
        """"""
        # radiance to reflectance taken as in mipp/xrit/MSG.py
        # again FCI User Guide is not clear on how to do this

        sirr = self.nc[
            '/data/{}/measured/channel_effective_solar_irradiance'
            .format(key.name)][...]

        # reflectance = radiance / sirr * 100
        data.data[:] /= sirr
        data.data[:] *= 100
142622
def stSpectralFlux(X, X_prev):
    
    # compute the spectral flux as the sum of square distances:
    sumX = numpy.sum(X + eps)
    sumPrevX = numpy.sum(X_prev + eps)
    F = numpy.sum((X / sumX - X_prev/sumPrevX) ** 2)

    return F
36256
def get_spectral_index(src, egy):
    """"""
    delta = 1E-5
    f0 = src.spectrum()(pyLike.dArg(egy * (1 - delta)))
    f1 = src.spectrum()(pyLike.dArg(egy * (1 + delta)))

    if f0 > 0 and f1 > 0:
        gamma = np.log10(f0 / f1) / np.log10((1 - delta) / (1 + delta))
    else:
        gamma = np.nan

    return gamma
644642
def compute_vis_noise(msname, sefd, spw_id=0):
    """"""
    from pyrap.tables import table

    tab = table(msname)
    spwtab = table(msname + "/SPECTRAL_WINDOW")

    freq0 = spwtab.getcol("CHAN_FREQ")[spw_id, 0]
    wavelength = 300e+6/freq0
    bw = spwtab.getcol("CHAN_WIDTH")[spw_id, 0]
    dt = tab.getcol("EXPOSURE", 0, 1)[0]
    dtf = (tab.getcol("TIME", tab.nrows()-1, 1)-tab.getcol("TIME", 0, 1))[0]

    # close tables properly, else the calls below will hang waiting for a lock...
    tab.close()
    spwtab.close()

    print(">>> %s freq %.2f MHz (lambda=%.2fm), bandwidth %.2g kHz, %.2fs integrations, %.2fh synthesis"%(msname, freq0*1e-6, wavelength, bw*1e-3, dt, dtf/3600))
    noise = sefd/math.sqrt(abs(2*bw*dt))
    print(">>> SEFD of %.2f Jy gives per-visibility noise of %.2f mJy"%(sefd, noise*1000))

    return noise
591775
protected void computePsiDataPsiGradient(GrayF32 image1, GrayF32 image2,
											 GrayF32 deriv1x, GrayF32 deriv1y,
											 GrayF32 deriv2x, GrayF32 deriv2y,
											 GrayF32 deriv2xx, GrayF32 deriv2yy, GrayF32 deriv2xy,
											 GrayF32 du, GrayF32 dv,
											 GrayF32 psiData, GrayF32 psiGradient ) {
		int N = image1.width * image1.height;

		for( int i = 0; i < N; i++ ) {

			float du_ = du.data[i];
			float dv_ = dv.data[i];

			// compute Psi-data
			float taylor2 = image2.data[i] + deriv2x.data[i]*du_ + deriv2y.data[i]*dv_;
			float v = taylor2 - image1.data[i];

			psiData.data[i] = (float)(1.0/(2.0*Math.sqrt(v*v + EPSILON*EPSILON)));

			// compute Psi-gradient
			float dIx = deriv2x.data[i] + deriv2xx.data[i]*du_ + deriv2xy.data[i]*dv_ - deriv1x.data[i];
			float dIy = deriv2y.data[i] + deriv2xy.data[i]*du_ + deriv2yy.data[i]*dv_ - deriv1y.data[i];
			float dI2 = dIx*dIx +  dIy*dIy;

			psiGradient.data[i] = (float)(1.0/(2.0*Math.sqrt(dI2 + EPSILON*EPSILON)));
		}
	}
216140
def cfn_viz(template, parameters={}, outputs={}, out=sys.stdout):
    
    known_sg, open_sg = _analyze_sg(template['Resources'])
    (graph, edges) = _extract_graph(template.get('Description', ''),
                                    template['Resources'], known_sg, open_sg)
    graph['edges'].extend(edges)
    _handle_terminals(template, graph, 'Parameters', 'source', parameters)
    _handle_terminals(template, graph, 'Outputs', 'sink', outputs)
    graph['subgraphs'].append(_handle_pseudo_params(graph['edges']))

    _render(graph, out=out)
763626
def _counts2radiance(self, counts, coefs, channel):
        """"""
        logger.debug('Converting counts to radiance')

        if self._is_vis(channel):
            # Since the scanline-detector assignment is unknown, use the average
            # coefficients for all scanlines.
            slope = np.array(coefs['slope']).mean()
            offset = np.array(coefs['offset']).mean()
            return self._viscounts2radiance(counts=counts, slope=slope,
                                            offset=offset)

        return self._ircounts2radiance(counts=counts, scale=coefs['scale'],
                                       offset=coefs['offset'])
852919
def overlay_spectra(model, dataset):
    
    best_flux, best_ivar = draw_spectra(model, dataset)
    coeffs_all, covs, scatters, all_chisqs, pivots, label_vector = model.model

    # Overplot original spectra with best-fit spectra
    print("Overplotting spectra for ten random stars")
    res = dataset.test_flux-best_flux
    lambdas = dataset.wl
    npix = len(lambdas)
    nstars = best_flux.shape[0]
    pickstars = []
    for i in range(10):
        pickstars.append(random.randrange(0, nstars-1))
    for i in pickstars:
        print("Star %s" % i)
        ID = dataset.test_ID[i]
        spec_orig = dataset.test_flux[i,:]
        bad = dataset.test_flux[i,:] == 0
        lambdas = np.ma.array(lambdas, mask=bad, dtype=float)
        npix = len(lambdas.compressed())
        spec_orig = np.ma.array(dataset.test_flux[i,:], mask=bad)
        spec_fit = np.ma.array(best_flux[i,:], mask=bad)
        ivars_orig = np.ma.array(dataset.test_ivar[i,:], mask=bad)
        ivars_fit = np.ma.array(best_ivar[i,:], mask=bad)
        red_chisq = np.sum(all_chisqs[:,i], axis=0) / (npix - coeffs_all.shape[1])
        red_chisq = np.round(red_chisq, 2)
        fig,axarr = plt.subplots(2)
        ax1 = axarr[0]
        im = ax1.scatter(lambdas, spec_orig, label="Orig Spec",
                         c=1 / np.sqrt(ivars_orig), s=10)
        ax1.scatter(lambdas, spec_fit, label="Cannon Spec", c='r', s=10)
        ax1.errorbar(lambdas, spec_fit, 
                     yerr=1/np.sqrt(ivars_fit), fmt='ro', ms=1, alpha=0.7)
        ax1.set_xlabel(r"Wavelength $\lambda (\AA)$")
        ax1.set_ylabel("Normalized flux")
        ax1.set_title("Spectrum Fit: %s" % ID)
        ax1.set_title("Spectrum Fit")
        ax1.set_xlim(min(lambdas.compressed())-10, max(lambdas.compressed())+10)
        ax1.legend(loc='lower center', fancybox=True, shadow=True)
        ax2 = axarr[1]
        ax2.scatter(spec_orig, spec_fit, c=1/np.sqrt(ivars_orig), alpha=0.7)
        ax2.errorbar(spec_orig, spec_fit, yerr=1 / np.sqrt(ivars_fit),
                     ecolor='k', fmt="none", ms=1, alpha=0.7)
        #fig.subplots_adjust(right=0.8)
        #cbar_ax = fig.add_axes([0.85, 0.15, 0.05, 0.7])
        fig.colorbar()
        #fig.colorbar(
        #        im, cax=cbar_ax,
        #        label="Uncertainties on the Fluxes from the Original Spectrum")
        xlims = ax2.get_xlim()
        ylims = ax2.get_ylim()
        lims = [np.min([xlims, ylims]), np.max([xlims, ylims])]
        ax2.plot(lims, lims, 'k-', alpha=0.75)
        textstr = "Red Chi Sq: %s" % red_chisq
        props = dict(boxstyle='round', facecolor='palevioletred', alpha=0.5)
        ax2.text(0.05, 0.95, textstr, transform=ax2.transAxes, fontsize=14,
                 verticalalignment='top', bbox=props)
        ax2.set_xlim(xlims)
        ax2.set_ylim(ylims)
        ax2.set_xlabel("Orig Fluxes")
        ax2.set_ylabel("Fitted Fluxes")
        plt.tight_layout()
        filename = "best_fit_spec_Star%s.png" % i
        print("Saved as %s" % filename)
        fig.savefig(filename)
        plt.close(fig)
609462
def _compute_fluxes(self):
        ''' '''
        self.emission = self._compute_emission()
        self.emission_sfc = self._compute_emission_sfc()
        fromspace = self._from_space()
        self.flux_down = self.trans.flux_down(fromspace, self.emission)
        self.flux_reflected_up = self.trans.flux_reflected_up(self.flux_down, self.albedo_sfc)
        # this ensure same dimensions as other fields
        self.flux_to_sfc = self.flux_down[..., -1, np.newaxis]
        self.flux_from_sfc = (self.emission_sfc +
                              self.flux_reflected_up[..., -1, np.newaxis])
        self.flux_up = self.trans.flux_up(self.flux_from_sfc,
                            self.emission + self.flux_reflected_up[...,0:-1])
        self.flux_net = self.flux_up - self.flux_down
        # absorbed radiation (flux convergence) in W / m**2 (per band)
        self.absorbed = np.diff(self.flux_net, axis=-1)
        self.absorbed_total = np.sum(self.absorbed, axis=-1)
        self.flux_to_space = self._compute_flux_top()