Python XSpectrum1D.from_array方法代码示例

# 需要导入模块: from linetools.spectra.xspectrum1d import XSpectrum1D [as 别名]
# 或者: from linetools.spectra.xspectrum1d.XSpectrum1D import from_array [as 别名]
def parse_UVES_popler(hdulist):
    """ Read a spectrum from a UVES_popler-style fits file.
    """
    from linetools.spectra.xspectrum1d import XSpectrum1D

    hd = hdulist[0].header
    uwave = setwave(hd) * u.Angstrom
    co = hdulist[0].data[3]
    fx = hdulist[0].data[0] * co  #  Flux
    sig = hdulist[0].data[1] * co
    xspec1d = XSpectrum1D.from_array(uwave, u.Quantity(fx),
                                     uncertainty=StdDevUncertainty(sig))
    xspec1d.co = co
    return xspec1d

# 需要导入模块: from linetools.spectra.xspectrum1d import XSpectrum1D [as 别名]
# 或者: from linetools.spectra.xspectrum1d.XSpectrum1D import from_array [as 别名]
def voigt_model(spec, line, Npix=None, flg_ret=1):
    """Generates a Voigt model from a line or list of lines

    Parameters:
    ------------
        spec: wave array or Spectrum
        line: Abs_Line, List of Abs_line, or array of parameters
        flg_ret : int (1)  Byte-wise Flag for return
          1: vmodel
          2: tau

    ToDo:
        1.  May need to more finely sample the wavelength array

    JXP 01 Nov 2014
    """
    # Imports
    import copy
    from linetools.spectra.xspectrum1d import XSpectrum1D
    from linetools.spectralline import AbsLine
    from astropy.nddata import StdDevUncertainty

    # Spectrum input
    if isinstance(spec,np.ndarray):  # Standard wavelength array
        vmodel = XSpectrum1D.from_array(spec, np.ones(len(spec)), 
        	uncertainty=StdDevUncertainty(np.zeros(len(spec))))
    elif isinstance(spec,Spectrum1D):
        vmodel = copy.deepcopy(spec)
    else:
        raise ValueError('voigt_model: Unknown spectrum input')
    #xdb.set_trace()
        
    # Line input
    if isinstance(line,AbsLine):  # Single line as a Abs_Line Class
        par = [0*i for i in range(6)] # Dummy list
        par[0] = line.attrib['N'] # logN; Won't have units
        par[1] = line.attrib['z']
        par[2] = line.attrib['b'] # Should have units
        par[3] = line.wrest # Should have units
        par[4] = line.data['f']
        par[5] = line.data['gamma'] # Should have units
    elif isinstance(line,list):
        if isinstance(line[0],AbsLine):  # List of Abs_Line
            tau = np.zeros(len(vmodel.flux))
            for iline in line:
                tau += voigt_model(vmodel.dispersion, iline, Npix=None, flg_ret=2) 
        else:
            par = line  # Single line as a vector
    else: 
        raise ValueError('voigt: Unknown type for voigt line')

    # tau
    if 'tau' not in locals():
        #xdb.set_trace()
        cold = 10.0**par[0] / u.cm / u.cm
        zp1=par[1]+1.0

        wv=par[3].to(u.cm) #*1.0e-8
        nujk = (const.c / wv).to(u.Hz)
        dnu = (par[2].to(u.km/u.s) / wv).to('Hz')
        avoigt = (par[5]/( 4 * np.pi * dnu)).to(u.dimensionless_unscaled)

        uvoigt = ( ((const.c / (vmodel.dispersion/zp1)) - nujk) / dnu).to(u.dimensionless_unscaled)
        '''
        wv=par[3].to(u.cm) #*1.0e-8
        bl=((par[2]*u.km/u.s)*wv/const.c).to(u.cm) # 2.99792458E5
        a= ((par[5]/u.s)*wv*wv/( 4 * np.pi * const.c.to('cm/s') * bl)).to(u.dimensionless_unscaled)
        cns=wv*wv*par[4]/(bl*2.002134602291006E12) * u.cm

        # Converting to Voigt units
        cne=cold*cns
        ww=(vmodel.dispersion.to('cm'))/zp1
        v=wv*ww*((1.0/ww)-(1.0/wv))/bl
        '''

        # Voigt
        cne = 0.01497 * cold * par[4] * u.cm * u.cm * u.Hz
        tau = cne * voigtking(uvoigt,avoigt) / (np.sqrt(np.pi) * dnu)

    # Flux
    if vmodel.unit is None:
        try:
            vmodel.flux = (np.exp(-1.0*tau)).value
        except AttributeError:
            vmodel.flux = np.exp(-1.0*tau)
    else:
        vmodel.flux = np.exp(-1.0*tau).to(vmodel.unit)

    # Convolve
    if Npix is not None:
        vmodel.gauss_smooth(npix=Npix)
    
    # Return
    ret_val = []
    if flg_ret % 2 == 1: ret_val.append(vmodel)
    if flg_ret % 4 >= 2: ret_val.append(tau)
    #xdb.set_trace()
    if len(ret_val) == 1: ret_val = ret_val[0]
    return ret_val

# 需要导入模块: from linetools.spectra.xspectrum1d import XSpectrum1D [as 别名]
# 或者: from linetools.spectra.xspectrum1d.XSpectrum1D import from_array [as 别名]

#.........这里部分代码省略.........
            else:
                uncertainty = None

            #Log-Linear?
            try:
                dc_flag = head0['DC-FLAG']
            except KeyError:
                # The following is necessary for Becker's XShooter output
                cdelt1, dc_flag = get_cdelt_dcflag(head0)

            # Read
            if dc_flag == 0:
                # Read FITS file
                spec1d = spec_read_fits.read_fits_spectrum1d(os.path.expanduser(datfil),
                                                             dispersion_unit='AA')
                spec1d.uncertainty = uncertainty
                xspec1d = XSpectrum1D.from_spec1d(spec1d)
            elif dc_flag == 1: # Generate wavelengths and use array approach
                fx = hdulist[0].data
                # Generate wave
                wave = setwave(head0)
            else:
                raise ValueError('DC-FLAG has unusual value {:d}'.format(dc_flag))

        elif hdulist[0].name == 'FLUX':
            # NEW SCHOOL (one file for flux and error)
            if 'WAVELENGTH' not in hdulist:
                spec1d = spec_read_fits.read_fits_spectrum1d(
                    os.path.expanduser(datfil), dispersion_unit='AA')
                xspec1d = XSpectrum1D.from_spec1d(spec1d)
            else:
                wave = hdulist['WAVELENGTH'].data * u.AA
                fx = hdulist['FLUX'].data
                xspec1d = XSpectrum1D.from_array(wave, u.Quantity(fx))

            # Error array
            if 'ERROR' in hdulist:
                sig = hdulist['ERROR'].data
                xspec1d.uncertainty = StdDevUncertainty(sig)
            else:
                sig = None

            if 'CONTINUUM' in hdulist:
                xspec1d.co = hdulist['CONTINUUM'].data

            if 'METADATA' in head0:
                xspec1d.meta.update(json.loads(head0['METADATA']))

        else:  # ASSUMING MULTI-EXTENSION
            if len(hdulist) <= 2:
                raise RuntimeError('No wavelength info but only 2 extensions!')
            fx = hdulist[0].data.flatten()
            try:
                sig = hdulist[1].data.flatten()
            except AttributeError:  # Error array is "None"
                sig = None
            wave = hdulist[2].data.flatten()
            # BOSS/SDSS?
            try:
                multi_ivar = head0['TELESCOP'][0:4] in ['SDSS']
            except KeyError:
                pass
            #
            if multi_ivar is True:
                tmpsig = np.zeros(len(sig))
                gdp = np.where(sig > 0.)[0]