python类AA的实例源码

def plot_profile(wave, flux, line, name):
    # 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)
    ax.set_ylim(0., 1.5)
    # make the plot
    ax.plot(vel_kms, flux)
    # include the name of the line
    plt.text(xmin+0.03*(xmax-xmin), 0.15, 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 plot_profile(wave, flux, line, name):
    # 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)
    ax.set_ylim(0., 3.)
    # make the plot
    ax.plot(vel_kms, np.log(1/flux))
    # include the name of the line
    plt.text(xmin+0.03*(xmax-xmin), 2.3, 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 get_stellar_mass(self, z=0, rnorm=1., cosmo=Planck15):
        """
        Get M/L, L, stellar mass, as measured in R-band

        Returns: 
            M/Lr, log10(Lr), log10(stellar_mass)

        """
        from sedpy.observate import Filter
        import astropy.units as u
        import astropy.constants as const

        MLr, MLi, MLk = self.get_M2L()
        #plt.plot(self.wave, self.spec); print(MLr)

        rfilt = Filter('sdss_r0')

        norm_spec = self.get_model(in_place=False)*rnorm

        #rband_Jy = rfilt.obj_counts(self.wave, norm_spec)/rfilt.ab_zero_counts*3631
        #rband_flam = rband_Jy/(3.34e4*rfilt.wave_pivot**2)#*u.erg/u.second/u.cm**2/u.AA
        #dL = Planck15.luminosity_distance(z)

        Mr = rfilt.ab_mag(self.wave, norm_spec) - cosmo.distmod(z).value
        Lr = 10**(-0.4*(Mr-rfilt.solar_ab_mag))
        #Lr = (rband_flam*4*np.pi*dL.to(u.cm).value**2)/3.828e+33*rfilt.wave_pivot*(1+z)
        stellar_mass = Lr*MLr
        return MLr, np.log10(Lr), np.log10(stellar_mass)

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 compute_ck2004_response(self, path, verbose=False):
        """
        Computes Castelli & Kurucz (2004) intensities across the entire
        range of model atmospheres.

        @path: path to the directory containing ck2004 SEDs
        @verbose: switch to determine whether computing progress should
        be printed on screen

        Returns: n/a
        """

        models = glob.glob(path+'/*M1.000*')
        Nmodels = len(models)

        # Store the length of the filename extensions for parsing:
        offset = len(models[0])-models[0].rfind('.')

        Teff, logg, abun = np.empty(Nmodels), np.empty(Nmodels), np.empty(Nmodels)
        InormE, InormP = np.empty(Nmodels), np.empty(Nmodels)

        if verbose:
            print('Computing Castelli & Kurucz (2004) passband intensities for %s:%s. This will take a while.' % (self.pbset, self.pbname))

        for i, model in enumerate(models):
            #~ spc = np.loadtxt(model).T -- waaay slower
            spc = np.fromfile(model, sep=' ').reshape(-1,2).T

            Teff[i] = float(model[-17-offset:-12-offset])
            logg[i] = float(model[-11-offset:-9-offset])/10
            sign = 1. if model[-9-offset]=='P' else -1.
            abun[i] = sign*float(model[-8-offset:-6-offset])/10

            spc[0] /= 1e10 # AA -> m
            spc[1] *= 1e7  # erg/s/cm^2/A -> W/m^3
            wl = spc[0][(spc[0] >= self.ptf_table['wl'][0]) & (spc[0] <= self.ptf_table['wl'][-1])]
            fl = spc[1][(spc[0] >= self.ptf_table['wl'][0]) & (spc[0] <= self.ptf_table['wl'][-1])]
            fl *= self.ptf(wl)
            flP = fl*wl
            InormE[i] = np.log10(fl.sum()/self.ptf_area*(wl[1]-wl[0]))             # energy-weighted intensity
            InormP[i] = np.log10(flP.sum()/self.ptf_photon_area*(wl[1]-wl[0]))     # photon-weighted intensity
            if verbose:
                if 100*i % (len(models)) == 0:
                    print('%d%% done.' % (100*i/(len(models)-1)))

        # Store axes (Teff, logg, abun) and the full grid of Inorm, with
        # nans where the grid isn't complete.
        self._ck2004_axes = (np.unique(Teff), np.unique(logg), np.unique(abun))

        self._ck2004_energy_grid = np.nan*np.ones((len(self._ck2004_axes[0]), len(self._ck2004_axes[1]), len(self._ck2004_axes[2]), 1))
        self._ck2004_photon_grid = np.nan*np.ones((len(self._ck2004_axes[0]), len(self._ck2004_axes[1]), len(self._ck2004_axes[2]), 1))
        for i, I0 in enumerate(InormE):
            self._ck2004_energy_grid[Teff[i] == self._ck2004_axes[0], logg[i] == self._ck2004_axes[1], abun[i] == self._ck2004_axes[2], 0] = I0
        for i, I0 in enumerate(InormP):
            self._ck2004_photon_grid[Teff[i] == self._ck2004_axes[0], logg[i] == self._ck2004_axes[1], abun[i] == self._ck2004_axes[2], 0] = I0

        # Tried radial basis functions but they were just terrible.
        #~ self._log10_Inorm_ck2004 = interpolate.Rbf(self._ck2004_Teff, self._ck2004_logg, self._ck2004_met, self._ck2004_Inorm, function='linear')
        self.content.append('ck2004')
        self.atmlist.append('ck2004')

def ingest_lines(self):

        print("Ingesting lines from {}".format(self.data_source.short_name))

        for rdr in self.ion_readers:

            atomic_number = rdr.ion.Z
            ion_charge = rdr.ion.Ion - 1

            ion = Ion.as_unique(self.session, atomic_number=atomic_number, ion_charge=ion_charge)

            try:
                bound_lines = rdr.bound_lines
            except ChiantiIonReaderError:
                print("Lines not found for ion {} {}".format(convert_atomic_number2symbol(atomic_number), ion_charge))
                continue

            print("Ingesting lines for {} {}".format(convert_atomic_number2symbol(atomic_number), ion_charge))

            lvl_index2id = self.get_lvl_index2id(ion)

            for index, row in bound_lines.iterrows():

                # index: (lower_level_index, upper_level_index)
                lower_level_index, upper_level_index = index

                try:
                    lower_level_id = int(lvl_index2id.loc[lower_level_index])
                    upper_level_id = int(lvl_index2id.loc[upper_level_index])
                except KeyError:
                    raise IngesterError("Levels from this source have not been found."
                                        "You must ingest levels before transitions")

                # Create a new line
                line = Line(
                    lower_level_id=lower_level_id,
                    upper_level_id=upper_level_id,
                    data_source=self.data_source,
                    wavelengths=[
                        LineWavelength(quantity=row["wavelength"]*u.AA,
                                       data_source=self.data_source,
                                       medium=MEDIUM_VACUUM,
                                       method=row["method"])
                    ],
                    a_values=[
                        LineAValue(quantity=row["a_value"]*u.Unit("s**-1"),
                                   data_source=self.data_source)
                    ],
                    gf_values=[
                        LineGFValue(quantity=row["gf_value"],
                                    data_source=self.data_source)
                    ]
                )

                self.session.add(line)