import os import glob import sys from os import path import numpy as np import astropy.units as u from astropy.io import fits from scipy.stats import norm from scipy.interpolate import interp1d data_path = os.path.dirname(__file__) def readcol(filename, **kwargs): """ Tries to reproduce the simplicity of the IDL procedure READCOL. Given a file with some columns of strings and columns of numbers, this function extract the columns from a file and places them in Numpy vectors with the proper type: name, mass = readcol('prova.txt', usecols=(0, 2)) where the file prova.txt contains the following: ################## # name radius mass ################## abc 25. 36. cde 45. 56. rdh 55 57. qtr 75. 46. hdt 47. 56. ################## This function is a wrapper for numpy.genfromtxt() and accepts the same input. See the following website for the full documentation https://docs.scipy.org/doc/numpy/reference/generated/numpy.genfromtxt.html """ f = np.genfromtxt(filename, dtype=None, **kwargs) t = type(f[0]) if t == np.ndarray or t == np.void: # array or structured array f = map(np.array, zip(*f)) # In Python 3.x all strings (e.g. name='NGC1023') are Unicode strings by defauls. # However genfromtxt() returns byte strings b'NGC1023' for non-numeric columns. # To have the same behaviour in Python 3 as in Python 2, I convert the Numpy # byte string 'S' type into Unicode strings, which behaves like normal strings. # With this change I can read the string a='NGC1023' from a text file and the # test a == 'NGC1023' will give True as expected. if sys.version >= '3': f = [v.astype(str) if v.dtype.char=='S' else v for v in f] return f def log_rebin(lamRange, spec, oversample=False, velscale=None, flux=False): """ Logarithmically rebin a spectrum, while rigorously conserving the flux. Basically the photons in the spectrum are simply redistributed according to a new grid of pixels, with non-uniform size in the spectral direction. When the flux keyword is set, this program performs an exact integration of the original spectrum, assumed to be a step function within the linearly-spaced pixels, onto the new logarithmically-spaced pixels. The output was tested to agree with the analytic solution. :param lamRange: two elements vector containing the central wavelength of the first and last pixels in the spectrum, which is assumed to have constant wavelength scale! E.g. from the values in the standard FITS keywords: LAMRANGE = CRVAL1 + [0,CDELT1*(NAXIS1-1)]. It must be LAMRANGE[0] < LAMRANGE[1]. :param spec: input spectrum. :param oversample: Oversampling can be done, not to loose spectral resolution, especally for extended wavelength ranges and to avoid aliasing. Default: OVERSAMPLE=1 ==> Same number of output pixels as input. :param velscale: velocity scale in km/s per pixels. If this variable is not defined, then it will contain in output the velocity scale. If this variable is defined by the user it will be used to set the output number of pixels and wavelength scale. :param flux: (boolean) True to preserve total flux. In this case the log rebinning changes the pixels flux in proportion to their dLam so the following command will show large differences beween the spectral shape before and after LOG_REBIN: plt.plot(exp(logLam), specNew) # Plot log-rebinned spectrum plt.plot(np.linspace(lamRange[0], lamRange[1], spec.size), spec) By defaul, when this is False, the above two lines produce two spectra that almost perfectly overlap each other. :return: [specNew, logLam, velscale] """ lamRange = np.asarray(lamRange) assert len(lamRange) == 2, 'lamRange must contain two elements' assert lamRange[0] < lamRange[1], 'It must be lamRange[0] < lamRange[1]' s = spec.shape assert len(s) == 1, 'input spectrum must be a vector' n = s[0] if oversample: m = int(n*oversample) else: m = int(n) dLam = np.diff(lamRange)/(n - 1.) # Assume constant dLam lim = lamRange/dLam + [-0.5, 0.5] # All in units of dLam borders = np.linspace(*lim, num=n+1) # Linearly logLim = np.log(lim) c = 299792.458 # Speed of light in km/s if velscale is None: # Velocity scale is set by user velscale = np.diff(logLim)/m*c # Only for output else: logScale = velscale/c m = int(np.diff(logLim)/logScale) # Number of output pixels logLim[1] = logLim[0] + m*logScale newBorders = np.exp(np.linspace(*logLim, num=m+1)) # Logarithmically k = (newBorders - lim[0]).clip(0, n-1).astype(int) specNew = np.add.reduceat(spec, k)[:-1] # Do analytic integral specNew *= np.diff(k) > 0 # fix for design flaw of reduceat() specNew += np.diff((newBorders - borders[k])*spec[k]) if not flux: specNew /= np.diff(newBorders) # Output log(wavelength): log of geometric mean logLam = np.log(np.sqrt(newBorders[1:]*newBorders[:-1])*dLam) return specNew, logLam, velscale def gaussian_filter1d(spec, sig): """ Convolve a spectrum by a Gaussian with different sigma for every pixel. If all sigma are the same this routine produces the same output as scipy.ndimage.gaussian_filter1d, except for the border treatment. Here the first/last p pixels are filled with zeros. When creating a template library for SDSS data, this implementation is 60x faster than a naive for loop over pixels. :param spec: vector with the spectrum to convolve :param sig: vector of sigma values (in pixels) for every pixel :return: spec convolved with a Gaussian with dispersion sig """ sig = sig.clip(0.01) # forces zero sigmas to have 0.01 pixels p = int(np.ceil(np.max(3*sig))) m = 2*p + 1 # kernel size x2 = np.linspace(-p, p, m)**2 n = spec.size a = np.zeros((m, n)) for j in range(m): # Loop over the small size of the kernel a[j, p:-p] = spec[j:n-m+j+1] gau = np.exp(-x2[:, None]/(2*sig**2)) gau /= np.sum(gau, 0)[None, :] # Normalize kernel conv_spectrum = np.sum(a*gau, 0) return conv_spectrum def calibrate(wave, flux, mag, filtername = 'SLOAN_SDSS.r'): """ Calibrate the spectra according to the magnitude. :param wave: _description_ :type wave: _type_ :param flux: _description_ :type flux: _type_ :param mag: _description_ :type mag: _type_ :param filtername: _description_, defaults to './data/SLOAN_SDSS.r' :type filtername: str, optional :return: _description_ :rtype: _type_ """ # Loading response curve if filtername == '5100': wave0 = np.linspace(3000,10000,7000) response0 = np.zeros(7000) response0[(wave0 > 5050) & (wave0 < 5150)] = 1. else: filter_file = data_path + '/data/filter/' + filtername+'.filter' wave0, response0 = readcol(filter_file) # Setting the response func = interp1d(wave0, response0) response = np.copy(wave) ind_extra = (wave > max(wave0)) | (wave < min(wave0)) response[ind_extra] = 0 ind_inside = (wave < max(wave0)) & (wave > min(wave0)) response[ind_inside] = func(wave[ind_inside]) # Flux map of datacube for given filter band preflux = np.sum(flux * response * np.mean(np.diff(wave))) / np.sum(response * np.mean(np.diff(wave))) # Real flux from magnitude for given filter realflux = (mag * u.STmag).to(u.erg/u.s/u.cm**2/u.AA).value # Normalization flux_ratio = realflux / preflux flux_calibrate = flux * flux_ratio * 1e17 # Units: 10^-17 erg/s/A/cm^2 return flux_calibrate # ---------------- # Reddening Module def Calzetti_Law(wave, Rv = 4.05): """ Calzetti_Law Dust Extinction Curve of Calzetti et al. (2000) Parameters ---------- wave : float Wavelength Rv : float, optional Extinction coefficient, by default 4.05 Returns ------- float Extinction value corresponding to the input wavelength """ wave_number = 1./(wave * 1e-4) reddening_curve = np.zeros(len(wave)) idx = np.logical_and(wave >= 1200, wave <= 6300) reddening_curve[idx] = 2.659 * ( -2.156 + 1.509 * wave_number[idx] - 0.198 * \ (wave_number[idx] ** 2)) + 0.011 * (wave_number[idx] **3 ) + Rv idx = np.logical_and(wave >= 6300, wave <= 22000) reddening_curve[idx] = 2.659 * ( -1.857 + 1.040 * wave_number[idx]) + Rv return reddening_curve def reddening(wave, flux, ebv = 0.0, law = 'calzetti', Rv = 4.05): """ Reddening an input spectra through a given reddening curve. Parameters ---------- wave : float Wavelength of input spectra flux : float Flux of input spectra ebv : float, optional E(B-V) value, by default 0 law : str, optional Extinction curve, by default 'calzetti' Rv : float, optional Extinction coefficient, by default 4.05 Returns ------- float Flux of spectra after reddening """ curve = Calzetti_Law(wave, Rv = Rv) fluxNew = flux / (10. ** (0.4 * ebv * curve)) return fluxNew ################ # Emission Lines ################ def SingleEmissinoLine(wave, line_wave, FWHM_inst): """ Profile of single emission line (Gaussian profile) Parameters ---------- wave : float Wavelength of line_wave : float Wavelength of emission line at the line center FWHM_inst : float Intrinsic broadening of emission line, units: A Returns ------- float Spectra of single emission line """ sigma = FWHM_inst / 2.355 flux = norm.pdf(wave, line_wave, sigma) return flux class EmissionLineTemplate(): """ Template for the emission lines """ def __init__(self, instrument, lam_range = [500, 15000], dlam = 0.1, model = 'hii'): """ __init__ _summary_ Parameters ---------- lam_range : list, optional _description_, by default [500, 15000] dlam : float, optional _description_, by default 0.1 model : str, optional _description_, by default 'nlr' FWHM_inst : float, optional _description_, by default 0.5 """ self.lam_range = lam_range self.wave = np.arange(lam_range[0], lam_range[1], 0.1) self.FWHM_inst = instrument.inst_fwhm self.model = model # HII region model of fsps-cloudy if model == 'hii': # Loading emission line flux table flux_table_file = data_path + '/data/fsps.nebular.fits' line_table = fits.open(flux_table_file) # Emission line list line_list = line_table[1].data line_wave = line_list['Wave'] line_names = line_list['Name'] w = (line_wave > lam_range[0]) & (line_wave < lam_range[1]) self.line_names = line_names[w] self.line_wave = line_wave[w] # Make parameter grid grid = line_table[2].data self.zh_grid = grid['logZ'] # Narrow line region model of if model == 'nlr': # Loading emission line flux table flux_table_file = data_path + '/data/AGN.NLR.fits' line_table = fits.open(flux_table_file) # Emission line list line_list = line_table[1].data line_wave = line_list['Wave'] line_names = line_list['Name'] w = (line_wave > lam_range[0]) & (line_wave < lam_range[1]) self.line_names = line_names[w] self.line_wave = line_wave[w] # Make parameter grid grid = line_table[2].data self.zh_grid = grid['logZ'] # Flux ratio flux_ratio = line_table[3].data self.flux_ratio = flux_ratio # Make emission line nline = len(line_wave) for i in range(nline): if i==0: emission_line = SingleEmissinoLine(self.wave, line_wave[i], self.FWHM_inst) emission_lines = emission_line else: emission_line = SingleEmissinoLine(self.wave, line_wave[i], self.FWHM_inst) emission_lines = np.vstack((emission_lines, emission_line)) self.emission_lines = emission_lines.T class HII_Region(): """ Class for the spectra of HII region Parameters ---------- wave : _type_ _description_ temp : _type_ _description_ Halpha : int, optional _description_, by default 100 logZ : bool, optional _description_, by default False vel : int, optional _description_, by default 100 vdisp : int, optional _description_, by default 120 Ebv : float, optional _description_, by default 0.1 Raises ------ ValueError _description_ """ def __init__(self, inst, temp, halpha = 100, zh = 0, vel = 100, vdisp = 120, ebv = 0.1): if (zh > -2) & (zh < 0.5): indz = np.argmin(np.abs(zh - temp.zh_grid)) flux_ratio = temp.flux_ratio[indz, :] else: raise ValueError('The range of logZ is not correct!') # Make emission line spectra through adding emission lines emlines = temp.emission_lines * flux_ratio flux_combine = np.sum(emlines, axis = 1) flux_calibrate = flux_combine * halpha # Units: erg/s/A/cm^2 # Dust attenuation if np.isscalar(ebv): flux_dust = reddening(temp.wave, flux_calibrate, ebv = ebv) # Broadening caused by Velocity Dispersion velscale = 10 lam_range = [np.min(temp.wave), np.max(temp.wave)] flux_logwave, logLam = log_rebin(lam_range, flux_dust, velscale=velscale)[:2] sigma_gas = vdisp / velscale # in pixel sigma_LSF = temp.FWHM_inst / (np.exp(logLam)) * 3e5 / velscale # in pixel if sigma_gas>0: sigma_dif = np.zeros(len(flux_logwave)) idx = (sigma_gas > sigma_LSF) sigma_dif[idx] = np.sqrt(sigma_gas ** 2. - sigma_LSF[idx] ** 2.) idx = (sigma_gas <= sigma_LSF) sigma_dif[idx] = 0.1 flux_broad = gaussian_filter1d(flux_logwave, sigma_dif) # Redshift redshift = vel / 3e5 wave_r = np.exp(logLam) * (1 + redshift) flux_red = np.interp(inst.wave, wave_r, flux_broad) self.wave = inst.wave self.flux = flux_red ############# # AGN Spectra ############# class AGN_NLR(): """ Class for narrow line region of AGN Parameters ---------- wave : float Wavelength temp : class Class of emission line template Halpha : float, optional Total flux of Halpha emission line, units 1e-17 erg/s, by default 100 logZ : float, optional Gas-phase metallicity, by default False vel : float, optional Line-of-sight velocity, by default 100 vdisp : float, optional Velocity dispersion, by default 120 Ebv : float, optional Dust extinction, by default 0.1 Raises ------ ValueError _description_ """ def __init__(self, temp, instrument, halpha = 100, zh = False, vel = 100, vdisp = 120, ebv = 0.1): if (zh > -2) & (zh < 0.5): indz = np.argmin(np.abs(zh - temp.zh_grid)) flux_ratio = temp.flux_ratio[indz, :] else: raise ValueError('The value of logZ is not correct!') # Make emission line spectra through adding emission lines emlines = temp.emission_lines * flux_ratio flux_combine = np.sum(emlines, axis = 1) flux_calibrate = flux_combine * halpha # Units: 1e-17 erg/s/A/cm^2 # Dust attenuation if np.isscalar(ebv): flux_dust = reddening(temp.wave, flux_calibrate, ebv = ebv) # Broadening caused by Velocity Dispersion velscale = 10 lam_range = [np.min(temp.wave), np.max(temp.wave)] flux_logwave, logLam = log_rebin(lam_range, flux_dust, velscale=velscale)[:2] sigma_gas = vdisp / velscale # in pixel sigma_LSF = temp.FWHM_inst / (np.exp(logLam)) * 3e5 / velscale # in pixel if sigma_gas>0: sigma_dif = np.zeros(len(flux_logwave)) idx = (sigma_gas > sigma_LSF) sigma_dif[idx] = np.sqrt(sigma_gas ** 2. - sigma_LSF[idx] ** 2.) idx = (sigma_gas <= sigma_LSF) sigma_dif[idx] = 0.1 flux_broad = gaussian_filter1d(flux_logwave, sigma_dif) # Redshift redshift = vel / 3e5 wave_r = np.exp(logLam) * (1 + redshift) flux_red = np.interp(instrument.wave, wave_r, flux_broad) self.wave = instrument.wave self.flux = flux_red class AGN_BLR(): """ Class for the broad line region of AGN """ def __init__(self, instrument, hbeta_flux = 100, hbeta_fwhm = 2000, ebv = None, vel = 0., lam_range = [500, 15000]): wave_rest = np.arange(lam_range[0], lam_range[1], 0.1) line_names = ['Hepsilon', 'Hdelta', 'Hgamma', 'Hbeta', 'Halpha'] line_waves = [3970.079, 4101.742, 4340.471, 4861.333, 6562.819] line_ratio = [0.101, 0.208, 0.405, 1.000, 2.579] # From Ilic et al. (2006) # Make emission lines for i in range(len(line_names)): if i==0: emission_line = SingleEmissinoLine(wave_rest, line_waves[i], hbeta_fwhm / 3e5 * line_waves[i]) emission_lines = emission_line else: emission_line = SingleEmissinoLine(wave_rest, line_waves[i], hbeta_fwhm / 3e5 * line_waves[i]) emission_lines = np.vstack((emission_lines, emission_line)) emlines = emission_lines.T * line_ratio flux_combine = np.sum(emlines, axis = 1) # Flux callibration flux_calibrate = flux_combine * hbeta_flux # Units: 1e-17 erg/s/A/cm^2 # Dust attenuation if np.isscalar(ebv): flux_dust = reddening(wave_rest, flux_calibrate, ebv = ebv) else: flux_dust = flux_calibrate # Redshift redshift = vel / 3e5 wave_r = wave_rest * (1 + redshift) flux_red = np.interp(instrument.wave, wave_r, flux_dust) self.wave = instrument.wave self.flux = flux_red class AGN_FeII(): """ Class for FeII emission lines of AGN """ def __init__(self, instrument, hbeta_broad = 100, r4570 = 0.4, ebv = None, vel = 0.): """ _summary_ """ filename = data_path + '/data/FeII.AGN.fits' # Loading FeII template hdulist = fits.open(filename) data = hdulist[1].data wave_rest = data['WAVE'] flux_model = data['FLUX'] # Determine the flux of FeII Fe4570_temp = 100 Fe4570_model = hbeta_broad * r4570 Ratio_Fe4570 = Fe4570_model / Fe4570_temp # Flux calibration flux_calibrate = flux_model * Ratio_Fe4570 # Dust attenuation if np.isscalar(ebv): flux_dust = reddening(wave_rest, flux_calibrate, ebv = ebv) else: flux_dust = flux_calibrate # Redshift redshift = vel / 3e5 wave_r = wave_rest * (1 + redshift) flux_red = np.interp(instrument.wave, wave_r, flux_dust) self.wave = instrument.wave self.flux = flux_red class AGN_Powerlaw(): """ AGN_Powerlaw _summary_ """ def __init__(self, instrument, m5100 = 1000, alpha = -1.5, ebv = None, vel = 0.,): """ _summary_ """ wave_rest = np.linspace(1000,20000,10000) flux = wave_rest ** alpha # Flux calibration flux_calibrate = calibrate(wave_rest, flux, m5100, filtername='5100') # Dust attenuation if np.isscalar(ebv): flux_dust = reddening(wave_rest, flux_calibrate, ebv = ebv) else: flux_dust = flux_calibrate # Redshift redshift = vel / 3e5 wave_r = wave_rest * (1 + redshift) flux_red = np.interp(instrument.wave, wave_r, flux_dust) self.wave = instrument.wave self.flux = flux_red class AGN(): def __init__(self, instrument, nlr_template, bhmass = 1e6, edd_ratio = 0.05, halpha_broad = 100, halpha_narrow = 100, vdisp_broad = 2000, vdisp_narrow = 500, vel = 1000, zh = 0, ebv = 0.1, dist = 20): NLR = AGN_NLR(nlr_template, instrument, halpha = halpha_narrow, zh = zh, vel = vel, vdisp = vdisp_narrow, ebv = ebv) if halpha_broad > 0: BLR = AGN_BLR(instrument, hbeta_flux = halpha_broad / 2.579, hbeta_fwhm = vdisp_broad / 2.355, ebv = ebv, vel = vel) m5100 = BHmass_to_M5100(bhmass, edd_ratio = edd_ratio, dist = dist) PL = AGN_Powerlaw(instrument, m5100 = m5100, ebv = ebv, vel = vel) Fe = AGN_FeII(instrument, hbeta_broad = halpha_broad / 2.579, ebv = ebv, vel = vel) self.wave = instrument.wave self.flux = NLR.flux + PL.flux + Fe.flux if halpha_broad > 0: self.flux = self.flux + BLR.flux def BHmass_to_M5100(bhmass, edd_ratio = 0.05, dist = 21): """ Caculate luminosity at 5100A according to the black hole mass Parameters ---------- BHmass : float Black hole mass, unit: solar mass Edd_ratio : float, optional Eddtington ratio, by default 0.05 Dist : int, optional Distance to the black hole, by default 21 Returns ------- float Luminosity at 5100A """ # Calculate bolometric luminosity Ledd = 3e4 * bhmass Lbol = Ledd * edd_ratio # Convert bolometric luminosity to 5100A luminosity (Marconi et al. 2004) L5100 = Lbol / 10.9 M5100 = 4.86 - 2.5 * np.log10(L5100) m5100 = M5100 + 5. * np.log10(dist * 1e5) return m5100 ################# # Stellar Spectra ################# def age_metal(filename): """ Extract the age and metallicity from the name of a file of the MILES library of Single Stellar Population models as downloaded from http://miles.iac.es/ as of 2016 :param filename: string possibly including full path (e.g. 'miles_library/Mun1.30Zm0.40T03.9811.fits') :return: age (Gyr), [M/H] """ s = path.basename(filename) age = float(s[s.find("T")+1:s.find("_iPp0.00_baseFe.fits")]) metal = s[s.find("Z")+1:s.find("T")] if "m" in metal: metal = -float(metal[1:]) elif "p" in metal: metal = float(metal[1:]) else: raise ValueError("This is not a standard MILES filename") return age, metal class StellarContinuumTemplate(object): def __init__(self, instrument, velscale = 50, pathname = data_path + '/data/EMILES/Ech*_baseFe.fits', normalize = False): """ Produces an array of logarithmically-binned templates by reading the spectra from the Single Stellar Population (SSP) library by Vazdekis et al. (2010, MNRAS, 404, 1639) http://miles.iac.es/. The code checks that the model specctra form a rectangular grid in age and metallicity and properly sorts them in both parameters. The code also returns the age and metallicity of each template by reading these parameters directly from the file names. The templates are broadened by a Gaussian with dispersion sigma_diff = np.sqrt(、**2 - sigma_tem**2). Thie script relies on the files naming convention adopted by the MILES library, where SSP spectra have the form below Mun1.30Zm0.40T03.9811.fits This code can be easily adapted by the users to deal with other stellar libraries, different IMFs or different abundances. :param pathname: path with wildcards returning the list files to use (e.g. 'miles_models/Mun1.30*.fits'). The files must form a Cartesian grid in age and metallicity and the procedure returns an error if they are not. :param velscale: desired velocity scale for the output templates library in km/s (e.g. 60). This is generally the same or an integer fraction of the velscale of the galaxy spectrum. :param FWHM_gal: vector or scalar of the FWHM of the instrumental resolution of the galaxy spectrum in Angstrom. :param normalize: set to True to normalize each template to mean=1. This is useful to compute light-weighted stellar population quantities. :return: The following variables are stored as attributes of the miles class: .templates: array has dimensions templates[npixels, n_ages, n_metals]; .log_lam_temp: natural np.log() wavelength of every pixel npixels; .age_grid: (Gyr) has dimensions age_grid[n_ages, n_metals]; .metal_grid: [M/H] has dimensions metal_grid[n_ages, n_metals]. .n_ages: number of different ages .n_metal: number of different metallicities """ FWHM_inst = instrument.inst_fwhm files = glob.glob(pathname) assert len(files) > 0, "Files not found %s" % pathname all = [age_metal(f) for f in files] all_ages, all_metals = np.array(all).T ages, metals = np.unique(all_ages), np.unique(all_metals) n_ages, n_metal = len(ages), len(metals) assert set(all) == set([(a, b) for a in ages for b in metals]), \ 'Ages and Metals do not form a Cartesian grid' # Extract the wavelength range and logarithmically rebin one spectrum # to the same velocity scale of the SDSS galaxy spectrum, to determine # the size needed for the array which will contain the template spectra. hdu = fits.open(files[0]) ssp = hdu[0].data h2 = hdu[0].header lam_range_temp = h2['CRVAL1'] + np.array([0, h2['CDELT1']*(h2['NAXIS1']-1)]) sspNew, log_lam_temp = log_rebin(lam_range_temp, ssp, velscale=velscale)[:2] #wave=((np.arange(hdr['NAXIS1'])+1.0)-hdr['CRPIX1'])*hdr['CDELT1']+hdr['CRVAL1'] templates = np.empty((sspNew.size, n_ages, n_metal)) age_grid = np.empty((n_ages, n_metal)) metal_grid = np.empty((n_ages, n_metal)) # Convolve the whole Vazdekis library of spectral templates # with the quadratic difference between the galaxy and the # Vazdekis instrumental resolution. Logarithmically rebin # and store each template as a column in the array TEMPLATES. # Quadratic sigma difference in pixels Vazdekis --> galaxy # The formula below is rigorously valid if the shapes of the # instrumental spectral profiles are well approximated by Gaussians. # FWHM of Emiles templates Emile_wave = np.exp(log_lam_temp) Emile_FWHM = np.zeros(h2['NAXIS1']) Emile_FWHM[np.where(Emile_wave < 3060)] = 3. Emile_FWHM[np.where((Emile_wave >= 3060) & (Emile_wave < 3540))] = 3. Emile_FWHM[np.where((Emile_wave >= 3540) & (Emile_wave < 8950))] = 2.5 Lwave = Emile_wave[np.where(Emile_wave >= 8950)] Emile_FWHM[np.where(Emile_wave >= 8950)]=60*2.35/3.e5*Lwave # sigma=60km/s at lambda > 8950 LSF = Emile_FWHM FWHM_eff=Emile_FWHM.copy() # combined FWHM from stellar library and instrument(input) if np.isscalar(FWHM_inst): FWHM_eff[Emile_FWHM < FWHM_inst] = FWHM_inst LSF[Emile_FWHM < FWHM_inst] = FWHM_inst #else: # FWHM_eff[Emile_FWHM < FWHM_inst] = FWHM_inst[Emile_FWHM < FWHM_inst] # LSF[Emile_FWHM < FWHM_inst] = FWHM_inst[Emile_FWHM < FWHM_inst] FWHM_dif = np.sqrt(FWHM_eff**2 - Emile_FWHM**2) sigma_dif = FWHM_dif/2.355/h2['CDELT1'] # Sigma difference in pixels # Here we make sure the spectra are sorted in both [M/H] and Age # along the two axes of the rectangular grid of templates. for j, age in enumerate(ages): for k, metal in enumerate(metals): p = all.index((age, metal)) hdu = fits.open(files[p]) ssp = hdu[0].data #if np.isscalar(FWHM_dif): # ssp = ndimage.gaussian_filter1d(ssp, sigma_dif) #else: ssp = gaussian_filter1d(ssp, sigma_dif) # convolution with variable sigma sspNew = log_rebin(lam_range_temp, ssp, velscale=velscale)[0] #if normalize: # sspNew /= np.mean(sspNew) templates[:, j, k] = sspNew age_grid[j, k] = age metal_grid[j, k] = metal self.templates = templates/np.median(templates) # Normalize by a scalar self.log_lam_temp = log_lam_temp self.age_grid = age_grid self.metal_grid = metal_grid self.n_ages = n_ages self.n_metal = n_metal self.lsf = log_rebin(lam_range_temp, LSF, velscale=velscale)[0] self.velscale = velscale def fmass_ssp(self): isedpath = data_path + '/data/EMILES/model/' massfile = isedpath + 'out_mass_CH_PADOVA00' n_metal = self.n_metal n_ages = self.n_ages fage = self.age_grid[:,0] fMs = np.zeros((n_ages,n_metal)) Metal, Age, Ms = readcol(massfile, usecols=(2, 3, 6)) for i in range(n_metal): for j in range(self.n_ages): locmin = np.argmin(abs(Metal - self.metal_grid[j, i])) & np.argmin(abs(Age - self.age_grid[j, i])) fMs[j,i] = Ms[locmin] return fMs class StellarContinuum(): def __init__(self, template, instrument, mag = 15, age = 1, feh = 0, vel = 100, vdisp = 100, ebv = 0): """ Modelling the spectra of stellar population Pars: ssp - The class of simple stellar population obswave - Wavelength, 1D-array mag - Magnitude in r-band used for the calibration of spectra, scalar Age - Mean age of stellar population used for determining the spectra profile, scalar (Gyr) FeH - Mean FeH of stellar population used for determining the spectra profile, scalar vel - Line of sight velocity used for determining the doppler effect, scalar (km/s) vdisp - Line of sight velocity dispersion used for broadening the spectra, scalar (km/s) Ebv - Extinction, scalar (mag), default 0 (no dust extinction) dFeH - Range of FeH when searching the best templates, default 0.2 (searching the template with -0.2 < FeH - dFeh < 0.2) Attri: wave - Wavelength is the same with obswave in Pars flux - Flux of best model """ # ----------------- # Stellar Continuum SSP_temp = template.templates # Select metal bins metals = template.metal_grid[0,:] minloc = np.argmin(abs(feh - metals)) tpls = SSP_temp[:, :, minloc] fmass = template.fmass_ssp()[:, minloc] # Select age bins Ages = template.age_grid[:,0] minloc = np.argmin(abs(age-Ages)) Stellar = tpls[:, minloc] wave = np.exp(template.log_lam_temp) # Broadening caused by Velocity Dispersion sigma_gal = vdisp / template.velscale # in pixel sigma_LSF = template.lsf / template.velscale # in pixel if sigma_gal>0: sigma_dif = np.zeros(len(Stellar)) idx = (sigma_gal > sigma_LSF) sigma_dif[idx] = np.sqrt(sigma_gal ** 2. - sigma_LSF[idx] ** 2.) idx = (sigma_gal <= sigma_LSF) sigma_dif[idx] = 0.1 flux0 = gaussian_filter1d(Stellar, sigma_dif) # Dust Reddening if np.isscalar(ebv): flux0 = reddening(wave, flux0, ebv = ebv) # Redshift redshift = vel / 3e5 wave_r = wave * (1 + redshift) flux = np.interp(instrument.wave, wave_r, flux0) # Calibration if np.isscalar(mag): flux = calibrate(instrument.wave, flux, mag, filtername='SLOAN_SDSS.r') # Convert to input wavelength self.wave = instrument.wave self.flux = flux ##################### # Single Star Spectra ##################### class SingleStarTemplate(): def __init__(self, inst, velscale = 20): filename = data_path + '/data/Starlib.XSL.fits' fwhm_inst = inst.inst_fwhm hdulist = fits.open(filename) lam = hdulist[1].data['Wave'] flux = hdulist[2].data par = hdulist[3].data lam_range_temp = np.array([3500, 12000]) TemNew, log_lam_temp = log_rebin(lam_range_temp, flux[1, :], velscale = velscale)[:2] # FWHM of XLS templates Temp_wave = np.exp(log_lam_temp) Temp_FWHM = np.zeros(len(log_lam_temp)) Temp_FWHM[(Temp_wave < 5330)] = 13 * 2.35 / 3e5 * Temp_wave[(Temp_wave < 5330)] # sigma = 13km/s at lambda <5330 Temp_FWHM[(Temp_wave >= 5330) & (Temp_wave < 9440)] = 11 * 2.35 / 3e5 * Temp_wave[(Temp_wave >= 5330) & (Temp_wave < 9440)] # sigma = 13km/s at 5330 < lambda < 9440 Temp_FWHM[(Temp_wave >= 9440)] = 16 * 2.35 / 3e5 * Temp_wave[(Temp_wave >= 9440)] # sigma=16km/s at lambda > 9440 LSF = Temp_FWHM FWHM_eff = Temp_FWHM.copy() # combined FWHM from stellar library and instrument(input) if np.isscalar(fwhm_inst): FWHM_eff[Temp_FWHM < fwhm_inst] = fwhm_inst LSF[Temp_FWHM < fwhm_inst] = fwhm_inst #else: # FWHM_eff[Temp_FWHM < FWHM_inst] = FWHM_inst[Temp_FWHM < FWHM_inst] # LSF[Temp_FWHM < FWHM_inst] = FWHM_inst[Temp_FWHM < FWHM_inst] FWHM_dif = np.sqrt(FWHM_eff ** 2 - Temp_FWHM ** 2) sigma_dif = FWHM_dif / 2.355 / (lam[1] - lam[0]) # Sigma difference in pixels temp = np.empty((TemNew.size, par.size)) for i in range(par.size): temp0 = log_rebin(lam_range_temp, flux[i, :], velscale=velscale)[0] #if np.isscalar(FWHM_dif): # temp1 = ndimage.gaussian_filter1d(temp0, sigma_dif) #else: temp1 = gaussian_filter1d(temp0, sigma_dif) # convolution with variable sigma tempNew = temp1 / np.mean(temp1) temp[:, i] = tempNew self.templates = temp self.log_lam_temp = log_lam_temp self.teff_grid = par['Teff'] self.feh_grid = par['FeH'] self.logg_grid = par['logg'] self.lsf = Temp_FWHM self.velscale = velscale class SingleStar(): def __init__(self, template, instrument, mag = 15, teff = 10000, feh = 0, vel = 100, ebv = 0): StarTemp = template.templates # Select metal bins idx_feh = (np.abs(template.feh_grid - feh) < 0.5) tpls = StarTemp[:, idx_feh] # Select Teff bins teff_feH = template.teff_grid[idx_feh] minloc = np.argmin(abs(teff - teff_feH)) starspec = tpls[:, minloc] wave = np.exp(template.log_lam_temp) # Dust Reddening if np.isscalar(ebv): starspec = reddening(wave, starspec, ebv = ebv) # Redshift redshift = vel / 3e5 wave_r = wave * (1 + redshift) flux = np.interp(instrument.wave, wave_r, starspec) # Calibration if np.isscalar(mag): flux = calibrate(instrument.wave, flux, mag, filtername='SLOAN_SDSS.r') # Convert to input wavelength self.wave = instrument.wave self.flux = flux