python类s()的实例源码

def __call__(self, x, *args, **kwargs):
        """
        Get the spline value and uncertainty at point(s) x. args and kwargs are passed to spline.__call__
        """
        x = np.atleast_1d(x)
        s = np.vstack([curve(x, *args, **kwargs) for curve in self._splines])
        return (np.mean(s, axis=0), np.std(s, axis=0))

def plot_profile(wave, flux, line, name, fosc):
    # define the velocity scale [AA / s]
    vel_aas = (wave - line*(1+zem[i])) / (line*(1+zem[i])) * c
    # convert to [km / s]
    vel_kms = vel_aas.to('km/s')
    # define parameters for the x and y axes
    ax.set_xlim(xmin, xmax)
    ax.xaxis.set_minor_locator(minorLocator)
    ax.tick_params(axis='x', labelsize=xls)
    ymax = 3.e12
    ax.set_ylim(0., ymax)
    # make the plot (using equations 5 and 8 from Savage & Sembach 1991)
    ax.plot(vel_kms, np.log(1/flux) / 2.654e-15 / (fosc * line))
    # include the name of the line
    plt.text(xmin+0.03*(xmax-xmin), ymax*0.6, name)
    # mark the approximate centroid velocity
    plt.axvline(x=vcen[i], ymin=0., ymax = 1.5, linewidth=1, color='k', linestyle='dotted')
    plt.axvline(x=vcen[i]+30., ymin=0., ymax = 1.5, linewidth=0.5, color='k')
    plt.axvline(x=vcen[i]-30., ymin=0., ymax = 1.5, linewidth=0.5, color='k')
    # label other lines for clarity
    for k in range(0, len(lines)):
        vel_off_aas = (lines[k] - line) / line * c
        vel_off_kms = vel_off_aas.to('km/s') / (u.km / u.s)
        plt.axvline(x=vcen[i]+vel_off_kms, ymin=0., ymax = 1.5, linewidth=1, color='k', linestyle='dotted')

# define the data directory

def calculate(wave, flux, line, fosc):
    # define the velocity scale [AA / s]
    vel_aas = (wave - line*(1+zem[i])) / (line*(1+zem[i])) * c
    # convert to [km / s]
    vel_kms = vel_aas.to('km/s')
    # label other lines for clarity

    tau = np.log(1/flux) 

    for k in range(0, len(lines)):
        vel_off_aas = (lines[k] - line) / line * c
        vel_off_kms = vel_off_aas.to('km/s') / (u.km / u.s)
    return(tau)

def scat_plot():
    f = plt.figure()
    # filename = 'MLP5_dap_multi_' + str(plate) + '_quicklook.pdf'
    mpl5_dir = os.environ['MANGADIR_MPL5']
    drp = fits.open(mpl5_dir + 'drpall-v2_0_1.fits')
    drpdata = drp[1].data
    absmag = drpdata.field('nsa_elpetro_absmag')


    plt.xlim(-16,-24)
    plt.ylim(1,7)
    plt.scatter(absmag[:,5], absmag[:,1]-absmag[:,5], marker='.',color=['blue'], s=0.5)
    plt.xlabel('i-band absolute magnitude', fontsize=16)
    plt.ylabel('NUV - i', fontsize=16)
    plt.tick_params(axis='both', labelsize=14)

    ifu_list = drpdata.field('plateifu')
    for i in good_galaxies:
        ithname = str(i[0]) + str(i[1])
        for e in range(0, len(ifu_list)):
            ethname = ifu_list[e]
            ethname = ethname.replace("-","")
            if ithname == ethname:
                plt.scatter(absmag[e, 5], absmag[e, 1] - absmag[e, 5], marker='*',color=['red'])
    f.savefig("scatter.pdf", bbox_inches='tight')
    # pp = PdfPages('scatter.pdf')
    # pp.savefig(plot_1)
    plt.close()
    os.system("open %s &" % 'scatter.pdf')

def luminosity(wavelength,flux,lDistcm):
    return const.c*u.s/u.m/(wavelength*10**-9)*flux*10**-23*(4*math.pi*lDistcm**2)/solarLum

# Given two coefficients and the color will return the mass-to-light ratio
# as defined by Bell and DeJong 2000

def luminosity(wavelength,flux,lDcm):
    return const.c*u.s/u.m/(wavelength*10**-9)*flux*10**-23*(4*math.pi*lDcm**2)/solarLum

def __init__(self, name, z, x,y):
        self.z = z
        self.y = y
        self.x = x
        self.name = name
        #lDcm gives the luminosity distance in centimeters
        #self.lDcm = cosmo.luminosity_distance(0.603)*u.Mpc.to(u.cm)/u.Mpc
        #self.lDcm = cosmo.luminosity_distance(self.z)*u.Mpc.to(u.cm) / u.Mpc
        #pixtoKpc is the conversion betweem pixels to rad to kpc based on z value
        #self.pixtoKpc = conv.arcsec_per_kpc_proper(self.z)*0.05/u.arcsec*u.kpc
        self.lDcm = cosmo.luminosity_distance(self.z)*u.Mpc.to(u.cm) / u.Mpc
        self.radToKpc = conv.arcsec_per_kpc_proper(self.z)*0.05/u.arcsec*u.kpc
# Let's set up some stuff so that I can grab stuff from inside an excel sheet.

def luminosity(wavelength,flux,lDistcm):
    return const.c*u.s/u.m/(wavelength*10**-9)*flux*10**-23*(4*math.pi*lDistcm**2)/solarLum

# Given two coefficients and the color will return the mass-to-light ratio
# as defined by Bell and DeJong 2000

def color_gui(self, cmap):
        """
        Function to change the cmap of the canvas
        :param cmap: string. matplotlib's color map. cmap = 'none' to gray scale again
        :return:
        """
        if cmap == 'none':
            self.color = False
            self.cmap = ""
        else:
            self.color = True
            self.cmap = cmap
        self.reload_canvas()

def load_data(self):
        hdulist = fits.open(self.filename)
        print("MuseCube: Loading the cube fluxes and variances...")

        # import pdb; pdb.set_trace()
        self.cube = hdulist[1].data
        self.stat = hdulist[2].data


        # for ivar weighting ; consider creating it in init ; takes long
        # self.flux_over_ivar = self.cube / self.stat

        self.header_1 = hdulist[1].header  # Necesito el header para crear una buena copia del white.
        self.header_0 = hdulist[0].header

        if self.filename_white is None:
            print("MuseCube: No white image given, creating one.")

            w_data = self.create_white(save=False)

            w_header_0 = copy.deepcopy(self.header_0)
            w_header_1 = copy.deepcopy(self.header_1)

            # These loops remove the third dimension from the header's keywords. This is neccesary in order to
            # create the white image and preserve the cube astrometry
            for i in w_header_0.keys():
                if '3' in i:
                    del w_header_0[i]
            for i in w_header_1.keys():
                if '3' in i:
                    del w_header_1[i]

            # prepare the header
            hdu = fits.HDUList()
            hdu_0 = fits.PrimaryHDU(header=w_header_0)
            hdu_1 = fits.ImageHDU(data=w_data, header=w_header_1)
            hdu.append(hdu_0)
            hdu.append(hdu_1)
            hdu.writeto('new_white.fits', clobber=True)
            self.filename_white = 'new_white.fits'
            print("MuseCube: `new_white.fits` image saved to disk.")

def color_gui(self, cmap):
        """
        Function to change the cmap of the canvas
        :param cmap: string. matplotlib's color map. cmap = 'none' to gray scale again
        :return:
        """
        if cmap == 'none':
            self.color = False
            self.cmap = ""
        else:
            self.color = True
            self.cmap = cmap
        self.reload_canvas()

def test_column_density_small_aperture(self):
        """Test column density for aperture << lengthscale.

        Should be within 1% of ideal value.

        """
        Q = 1e28 / u.s
        v = 1 * u.km / u.s
        rho = 10 * u.km
        parent = 1e4 * u.km
        N_avg = 2 * Haser(Q, v, parent).column_density(rho)
        ideal = Q / v / 2 / rho
        assert np.isclose(N_avg.decompose().value, ideal.decompose().value,
                          rtol=0.01)

def test_total_number_large_aperture(self):
        """Test column density for aperture >> lengthscale."""
        Q = 1 / u.s
        v = 1 * u.km / u.s
        rho = 1000 * u.km
        parent = 10 * u.km
        N = Haser(Q, v, parent).total_number(rho)
        ideal = Q * parent / v
        assert np.isclose(N, ideal.decompose().value)

    # TEST FAILS
    #
    # def test_total_number_rho_NJ78(self):
    #     """Reproduce Newburn and Johnson 1978.

    #     Species, N observed, Q/v (km**-1)
    #     CN, 6.4e26, 5.8e23
    #     C3, 8.3e28, 9.0e23
    #     C2, 7.8e27, 5.9e24

    #     Cannot reproduce C3 quoted in paper.

    #     """

    #     #Nobs = [6.41756750e26, 8.63191842e+28, 7.81278300e27]
    #     #Nobs = [6.4e26, 4.2e27, 7.8e27]
    #     Nobs = [6.4e26, 8.3e28, 7.8e27]
    #     parent = [1.4e4, 0, 1.0e4] * u.km
    #     daughter = [1.7e5, 4.6e4, 7.6e4] * u.km
    #     Q = [5.8e23, 9.0e23, 5.9e24] / u.s
    #     rho = 3300 * u.km

    #     N = np.zeros(3)
    #     for i in range(3):
    #         coma = Haser(Q[i], 1 * u.km / u.s, parent[i], daughter[i])
    #         N[i] = coma.total_number(rho)

    #     assert np.allclose(N, Nobs)

def test_circular_integration_0step(self):
        """Compare total number and integrated column density for circle."""

        from ..core import CircularAperture

        Nobs = 6.41756750e26
        parent = 1.4e4 * u.km
        Q = 5.8e23 / u.s
        v = 1 * u.km / u.s
        aper = CircularAperture(3300 * u.km)

        coma = Haser(Q, v, parent)
        N1 = coma.total_number(aper)
        N2 = coma._integrate_column_density(aper)
        assert np.isclose(N1, N2)

def test_total_number_rectangular_ap(self):
        """

        compare with:

        import astropy.units as u
        from sbpy.imageanalysis.utils import rarray
        from sbpy.activity import Haser

        r = rarray((5000, 3300), subsample=10) * u.km
        parent = 1.4e4 * u.km
        daughter = 1.7e5 * u.km
        Q = 5.8e23 / u.s
        v = 1 * u.km / u.s
        coma = Haser(Q, v, parent, daughter)
        sigma = coma.column_density(r)
        print((sigma * 1 * u.km**2).decompose())
        --> <Quantity 3.449607967230623e+26>

        This differs from the test value below by XXX

        """

        from ..core import RectangularAperture

        parent = 1.4e4 * u.km
        daughter = 1.7e5 * u.km
        Q = 5.8e23 / u.s
        v = 1 * u.km / u.s
        aper = RectangularAperture([5000, 3300] * u.km)

        coma = Haser(Q, v, parent, daughter)
        N = coma.total_number(aper)

        assert np.isclose(N, 3.449607967230623e+26)

def __init__(self, Q, v):
        assert isinstance(Q, u.Quantity)
        assert Q.unit.is_equivalent((u.s**-1, u.mol / u.s))
        self.Q = Q

        assert isinstance(v, u.Quantity)
        assert v.unit.is_equivalent(u.m / u.s)
        self.v = v

def __getattr__(self, field):
        if field in self.table.columns:
            return self.table[field]
        else:
           raise AttributeError ("field '{:s}' does not exist".format(field))

def calculate_counts_simple(channel, loop, *args):
        """
        Calculate the AIA intensity using only the temperature response functions.
        """
        response_function = (splev(np.ravel(loop.electron_temperature), channel['temperature_response_spline'])
                             * u.count*u.cm**5/u.s/u.pixel)
        counts = np.reshape(np.ravel(loop.density**2)*response_function, loop.density.shape)
        return counts

def __init__(self, observing_time: u.s, observing_area=None):
        self.logger = logging.getLogger(name=type(self).__name__)
        self.observing_time = np.arange(observing_time[0].to(u.s).value,
                                        observing_time[1].to(u.s).value,
                                        self.cadence.value)*u.s
        self.observing_area = observing_area

def make_aia_animation(aia, start_time: u.s, stop_time: u.s, root_dir, figsize=None, norm=None, fontsize=14, **kwargs):
    """
    Build animation from a series of synthesized AIA observations
    """
    with h5py.File(aia.counts_file, 'r') as hf:
        reference_time = u.Quantity(hf['time'], hf['time'].attrs['units'])
    start_index = np.where(reference_time == start_time)[0][0]
    stop_index = np.where(reference_time == stop_time)[0][0]
    fig_format = os.path.join(root_dir, f'{aia.name}', '{}', 'map_t{:06d}.fits')
    fig, ims = plot_aia_channels(aia, start_time, root_dir, figsize=figsize, norm=norm, fontsize=fontsize, 
                                 use_with_animation=True)

    def update_fig(i):
        for channel in aia.channels:
            tmp = Map(fig_format.format(channel['name'], i))
            ims[channel['name']].set_array(tmp.data)
        fig.suptitle(r'$t={:.0f}$ {}'.format(reference_time[i].value, reference_time.unit.to_string()), 
                     fontsize=fontsize)
        return [ims[k] for k in ims]

    animator_settings = {'interval': 50, 'blit': True}
    animator_settings.update(kwargs.get('animator_settings', {}))
    animation = matplotlib.animation.FuncAnimation(fig, update_fig, frames=range(start_index, stop_index),
                                                   **animator_settings)

    return animation