From 8afb62f5c642c4f436cfcde199ce15a9bc21187e Mon Sep 17 00:00:00 2001 From: Shuai Feng <451424498@qq.com> Date: Mon, 15 Apr 2024 15:40:29 +0800 Subject: [PATCH] 0415 first upload --- LICENSE | 19 + PKG-INFO | 9 + README.md | 4 +- setup.cfg | 8 + setup.py | 18 + src/csst_ifs_gehong.egg-info/PKG-INFO | 9 + src/csst_ifs_gehong.egg-info/SOURCES.txt | 17 + .../dependency_links.txt | 1 + src/csst_ifs_gehong.egg-info/requires.txt | 1 + src/csst_ifs_gehong.egg-info/top_level.txt | 1 + src/gehong/__init__.py | 0 src/gehong/config.py | 36 + src/gehong/cube3d.py | 94 ++ src/gehong/map2d.py | 362 ++++++ src/gehong/spec1d.py | 1101 +++++++++++++++++ 15 files changed, 1678 insertions(+), 2 deletions(-) create mode 100644 LICENSE create mode 100644 PKG-INFO create mode 100644 setup.cfg create mode 100644 setup.py create mode 100644 src/csst_ifs_gehong.egg-info/PKG-INFO create mode 100644 src/csst_ifs_gehong.egg-info/SOURCES.txt create mode 100644 src/csst_ifs_gehong.egg-info/dependency_links.txt create mode 100644 src/csst_ifs_gehong.egg-info/requires.txt create mode 100644 src/csst_ifs_gehong.egg-info/top_level.txt create mode 100644 src/gehong/__init__.py create mode 100644 src/gehong/config.py create mode 100644 src/gehong/cube3d.py create mode 100644 src/gehong/map2d.py create mode 100644 src/gehong/spec1d.py diff --git a/LICENSE b/LICENSE new file mode 100644 index 0000000..335ea9d --- /dev/null +++ b/LICENSE @@ -0,0 +1,19 @@ +Copyright (c) 2018 The Python Packaging Authority + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. \ No newline at end of file diff --git a/PKG-INFO b/PKG-INFO new file mode 100644 index 0000000..2de492c --- /dev/null +++ b/PKG-INFO @@ -0,0 +1,9 @@ +Metadata-Version: 2.1 +Name: csst-ifs-gehong +Version: 1.0.0 +Home-page: https://csst-ifs-gehong.readthedocs.io/en/latest/ +Author: Shuai Feng +Author-email: sfeng@hebtu.edu.cn +License: MIT +Keywords: CSST-IFS +License-File: LICENSE diff --git a/README.md b/README.md index 4416792..660348c 100644 --- a/README.md +++ b/README.md @@ -1,3 +1,3 @@ -# GEHONG +# csst-ifs-gehong Package -GEnerate tHe data Of iNtegral field spectrograph of Galaxy +This is a Python package for modelling the data of intergral field spectrascopy mounted on the Chinese Space Station Telescopy (CSST-IFS). See more detail at [csst-ifs-gehong](https://csst-ifs-gehong.readthedocs.io/). \ No newline at end of file diff --git a/setup.cfg b/setup.cfg new file mode 100644 index 0000000..bac36ec --- /dev/null +++ b/setup.cfg @@ -0,0 +1,8 @@ +[metadata] +description-file = README.md +license_files = LICENSE + +[egg_info] +tag_build = +tag_date = 0 + diff --git a/setup.py b/setup.py new file mode 100644 index 0000000..ce4b821 --- /dev/null +++ b/setup.py @@ -0,0 +1,18 @@ +from setuptools import setup, find_packages + + +setup( + name='csst-ifs-gehong', + version='1.0.0', + license='MIT', + author="Shuai Feng", + author_email='sfeng@hebtu.edu.cn', + packages=find_packages('src'), + package_dir={'': 'src'}, + url='https://csst-ifs-gehong.readthedocs.io/en/latest/', + keywords='CSST-IFS', + install_requires=[ + 'astropy', + ], + +) \ No newline at end of file diff --git a/src/csst_ifs_gehong.egg-info/PKG-INFO b/src/csst_ifs_gehong.egg-info/PKG-INFO new file mode 100644 index 0000000..2de492c --- /dev/null +++ b/src/csst_ifs_gehong.egg-info/PKG-INFO @@ -0,0 +1,9 @@ +Metadata-Version: 2.1 +Name: csst-ifs-gehong +Version: 1.0.0 +Home-page: https://csst-ifs-gehong.readthedocs.io/en/latest/ +Author: Shuai Feng +Author-email: sfeng@hebtu.edu.cn +License: MIT +Keywords: CSST-IFS +License-File: LICENSE diff --git a/src/csst_ifs_gehong.egg-info/SOURCES.txt b/src/csst_ifs_gehong.egg-info/SOURCES.txt new file mode 100644 index 0000000..46f205c --- /dev/null +++ b/src/csst_ifs_gehong.egg-info/SOURCES.txt @@ -0,0 +1,17 @@ +LICENSE +README.md +setup.cfg +setup.py +src/csst_ifs_gehong.egg-info/PKG-INFO +src/csst_ifs_gehong.egg-info/SOURCES.txt +src/csst_ifs_gehong.egg-info/dependency_links.txt +src/csst_ifs_gehong.egg-info/requires.txt +src/csst_ifs_gehong.egg-info/top_level.txt +src/gehong/__init__.py +src/gehong/config.py +src/gehong/cube3d.py +src/gehong/map2d.py +src/gehong/spec1d.py +test/test_cube3d.py +test/test_map2d.py +test/test_spec1d.py \ No newline at end of file diff --git a/src/csst_ifs_gehong.egg-info/dependency_links.txt b/src/csst_ifs_gehong.egg-info/dependency_links.txt new file mode 100644 index 0000000..8b13789 --- /dev/null +++ b/src/csst_ifs_gehong.egg-info/dependency_links.txt @@ -0,0 +1 @@ + diff --git a/src/csst_ifs_gehong.egg-info/requires.txt b/src/csst_ifs_gehong.egg-info/requires.txt new file mode 100644 index 0000000..ae8bed8 --- /dev/null +++ b/src/csst_ifs_gehong.egg-info/requires.txt @@ -0,0 +1 @@ +astropy diff --git a/src/csst_ifs_gehong.egg-info/top_level.txt b/src/csst_ifs_gehong.egg-info/top_level.txt new file mode 100644 index 0000000..70c83a1 --- /dev/null +++ b/src/csst_ifs_gehong.egg-info/top_level.txt @@ -0,0 +1 @@ +gehong diff --git a/src/gehong/__init__.py b/src/gehong/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/src/gehong/config.py b/src/gehong/config.py new file mode 100644 index 0000000..a0ab871 --- /dev/null +++ b/src/gehong/config.py @@ -0,0 +1,36 @@ +import numpy as np + +class config(): + """ + The configuration of spectral modeling. The default value is used for ETC calculation. + + Parameters + ---------- + wave_min : float, optional + Minimum value of wavelength coverage, by default 3500.0A + wave_max : float, optional + Minimum value of wavelength coverage, by default 10000.0A + dlam : float, optional + Wavelength width of each spaxel, by default 2.0A + inst_fwhm : float, optional + Spectral resolution, by default 0.1A + nx : int, optional + Number of spaxel in a spatial direction, by default 30 + ny : int, optional + Number of spaxel in a spatial direction, by default 30 + dpix : float, optional + Pixel size in the spatial direction, by default 0.2arcsec + """ + def __init__(self, wave_min = 3500.0, wave_max = 10000.0, + dlam = 2.0, inst_fwhm = 0.1, + nx = 30, ny = 30, dpix = 0.2): + self.dlam = dlam + self.wave = np.arange(wave_min, wave_max, dlam) + self.wave_min = wave_min + + self.inst_fwhm = inst_fwhm + self.nx = nx + self.ny = ny + self.dpix = dpix + self.fov_x = nx * dpix + self.fov_y = ny * dpix \ No newline at end of file diff --git a/src/gehong/cube3d.py b/src/gehong/cube3d.py new file mode 100644 index 0000000..1a92edd --- /dev/null +++ b/src/gehong/cube3d.py @@ -0,0 +1,94 @@ +import astropy.units as u +import numpy as np +from scipy.interpolate import interp1d +from .spec1d import * +import astropy.wcs + +class Cube3D(): + """ + Class of 3-dimentional spectral cube + """ + def __init__(self, config, stellar_map = None, gas_map = None): + + self.config= config + + self.nx = config.nx + self.ny = config.ny + self.dpix = config.dpix + self.fov_x = config.fov_x + self.fov_y = config.fov_y + + self.wave = config.wave + self.nz = len(self.wave) + self.wave0 = np.min(self.wave) + self.inst_fwhm = config.inst_fwhm + + self.flux = np.zeros((self.nx,self.ny,self.nz)) + + self.stellar_map = stellar_map + self.gas_map = gas_map + + def make_cube(self, stellar_tem = None, hii_tem = None): + + for i in range(self.nx): + for j in range(self.ny): + if self.stellar_map is not None: + ss = StellarContinuum(self.config, stellar_tem, mag = self.stellar_map.mag[i,j], + age = self.stellar_map.age[i,j], feh = self.stellar_map.feh[i,j], + vel = self.stellar_map.vel[i,j], vdisp = self.stellar_map.vdisp[i,j], + ebv = self.stellar_map.ebv[i,j]) + if self.gas_map is not None: + gg = HII_Region(self.config, hii_tem, halpha = self.gas_map.halpha[i,j], + logz = self.gas_map.zh[i,j], vel = self.gas_map.vel[i,j], + vdisp = self.gas_map.vdisp[i,j], ebv = self.gas_map.ebv[i,j]) + self.flux[i,j,:] = ss.flux + gg.flux + else: + self.flux[i,j,:] = ss.flux + + def wcs_info(self): + wcs = fits.Header() + wcs_dict = {'CTYPE1': 'WAVE ', + 'CUNIT1': 'Angstrom', + 'CDELT1': self.config.dlam, + 'CRPIX1': 1, + 'CRVAL1': np.min(self.wave), + 'CTYPE2': 'RA---TAN', + 'CUNIT2': 'deg', + 'CDELT2': self.dpix / 3600., + 'CRPIX2': np.round(self.ny / 2.), + 'CRVAL2': 0.5, + 'CTYPE3': 'DEC--TAN', + 'CUNIT3': 'deg', + 'CDELT3': self.dpix / 3600., + 'CRPIX3': np.round(self.nx / 2.), + 'CRVAL3': 1, + 'BUNIT' : '10**(-17)*erg/s/cm**2/Angstrom'} + input_wcs = astropy.wcs.WCS(wcs_dict) + self.wcs_header = input_wcs.to_header() + + def insert_spec(self, spec, dx = 0, dy = 0): + x0 = np.int(np.round(self.config.nx / 2.)) + y0 = np.int(np.round(self.config.ny / 2.)) + self.flux[x0 + dx, y0 + dy, :] = self.flux[x0 + dx, y0 + dy, :] + spec.flux + + def savefits(self, filename, path = './'): + + hdr = fits.Header() + + hdr['FILETYPE'] = 'SCICUBE' + hdr['CODE'] = 'CSST-IFS-GEHONG' + hdr['VERSION'] = '0.0.1' + + hdr['OBJECT'] = 'NGC1234' + hdr['RA'] = 0.0 + hdr['DEC'] = 0.0 + + hdu0 = fits.PrimaryHDU(header = hdr) + + self.wcs_info() + hdr = self.wcs_header + hdu1 = fits.ImageHDU(self.flux, header = hdr) + + # Output + hdulist = fits.HDUList([hdu0, hdu1]) + hdulist.writeto(path + filename, overwrite = True) \ No newline at end of file diff --git a/src/gehong/map2d.py b/src/gehong/map2d.py new file mode 100644 index 0000000..12c1e37 --- /dev/null +++ b/src/gehong/map2d.py @@ -0,0 +1,362 @@ +from __future__ import division + +import scipy.special as sp +import numpy as np +from astropy.io import fits +from skimage.transform import resize + +def Sersic2D(x, y, mag = 12.0, r_eff = 1.0, n = 2.0, ellip = 0.5, + theta = 0.0, x_0 = 0.0, y_0 = 0.0, pixelscale = 0.01): + """ + Model of 2D - Sersic profile + + Parameters + ---------- + x : float array + _description_ + y : _type_ + _description_ + mag : float, optional + Integral magnitude of sersic model, by default 12.0 + r_eff : float, optional + Effective radius in pixel, by default 1.0 + n : float, optional + Sersic index, by default 2.0 + ellip : float, optional + Ellipticity, by default 0.5 + theta : float, optional + Position angle in degree, by default 0.0 + x_0 : float, optional + Offset of the center of Sersic model, by default 0.0 + y_0 : float, optional + Offset of the center of Sersic model, by default 0.0 + pixelscale : float, optional + Size of each pixel in arcsec^2, by default 0.01 + + Returns + ------- + + _description_ + """ + # Produce Sersic profile + bn = sp.gammaincinv(2. * n, 0.5) + a, b = r_eff, (1 - ellip) * r_eff + cos_theta, sin_theta = np.cos(theta), np.sin(theta) + x_maj = (x - x_0) * cos_theta + (y - y_0) * sin_theta + x_min = -(x - x_0) * sin_theta + (y - y_0) * cos_theta + z = (abs(x_maj / a) ** 2 + abs(x_min / b) ** 2) ** (1 / 2) + profile = np.exp(-bn * (z ** (1 / n) - 1)) + + # Normalization + integral = a * b * 2 * np.pi * n * np.exp(bn) / (bn ** (2 * n)) * sp.gamma(2 * n) + prof_norm = profile / integral * pixelscale + + # Calibration + total_flux = 10. ** ((22.5 - mag) * 0.4) + sb_mag = 22.5 - 2.5 * np.log10(prof_norm * total_flux / pixelscale) + + return sb_mag + +def VelMap2D(x, y, vmax = 200.0, rt = 1.0, ellip = 0.5, + theta = 0.0, x_0 = 0.0, y_0 = 0.0): + """ + VelMap2D _summary_ + + Parameters + ---------- + x : _type_ + _description_ + y : _type_ + _description_ + vmax : int, optional + _description_, by default 200 + rt : int, optional + _description_, by default 1 + ellip : float, optional + _description_, by default 0.5 + theta : int, optional + _description_, by default 0 + x_0 : int, optional + _description_, by default 0 + y_0 : int, optional + _description_, by default 0 + + Returns + ------- + _type_ + _description_ + """ + # Produce tanh profile + a, b = rt, (1 - ellip) * rt + cos_theta, sin_theta = np.cos(theta), np.sin(theta) + x_maj = (x - x_0) * cos_theta + (y - y_0) * sin_theta + x_min = -(x - x_0) * sin_theta + (y - y_0) * cos_theta + z = (abs(x_maj / a) ** 2 + abs(x_min / b) ** 2) ** (1 / 2) + profile = vmax * np.tanh(z) * ((x_maj / a) / z) + + return profile + +def GradMap2D(x, y, a0 = 10, r_eff = 1, gred = -1, ellip = 0.5, + theta = 0, x_0 = 0, y_0 = 0): + """ + GradMap2D _summary_ + + Parameters + ---------- + x : _type_ + _description_ + y : _type_ + _description_ + a0 : int, optional + _description_, by default 10 + r_eff : int, optional + _description_, by default 1 + gred : int, optional + _description_, by default -1 + ellip : float, optional + _description_, by default 0.5 + theta : int, optional + _description_, by default 0 + x_0 : int, optional + _description_, by default 0 + y_0 : int, optional + _description_, by default 0 + + Returns + ------- + _type_ + _description_ + """ + # Produce gradiant profile + a, b = r_eff, (1 - ellip) * r_eff + cos_theta, sin_theta = np.cos(theta), np.sin(theta) + x_maj = (x - x_0) * cos_theta + (y - y_0) * sin_theta + x_min = -(x - x_0) * sin_theta + (y - y_0) * cos_theta + z = (abs(x_maj / a) ** 2 + abs(x_min / b) ** 2) ** (1 / 2) + profile = a0 + z * gred + + return profile + +class Map2d(object): + + def __init__(self, config): + """ + __init__ _summary_ + + Parameters + ---------- + inst : _type_ + _description_ + """ + self.xsamp = config.dpix + self.ysamp = config.dpix + startx = -(config.nx - 1) / 2.0 * self.xsamp + stopx = (config.nx - 1) / 2.0 * self.xsamp + starty = -(config.ny - 1) / 2.0 * self.ysamp + stopy = (config.ny - 1) / 2.0 * self.ysamp + xvals = np.linspace(startx, stopx, num = config.nx) + yvals = np.linspace(starty, stopy, num = config.ny) + + ones = np.ones((config.ny, config.nx)) + x = ones * xvals + y = np.flipud(ones * yvals.reshape(int(config.ny), 1)) + + self.nx = config.nx + self.ny = config.ny + self.x = x + self.y = y + self.row = xvals + # flip Y axis because we use Y increasing from bottom to top + self.col = yvals[::-1] + + def shift_rotate(self, yoff, xoff, rot): + """ + Return shifted/rotated (y, x) given offsets (yoff, xoff) and rotation, rot (degrees) + + Parameters + ---------- + yoff, xoff: float + yoff, xoff offsets in world coordinates + rot: float + rotation angle in degrees + + Returns + ------- + ysh_rot, xsh_rot: 2D numpy arrays + rotated and shifted copies of Grid.x and Grid.y + """ + pa_radians = np.pi * rot / 180.0 + xsh = self.x - xoff + ysh = self.y - yoff + xsh_rot = xsh * np.cos(pa_radians) + ysh * np.sin(pa_radians) + ysh_rot = -xsh * np.sin(pa_radians) + ysh * np.cos(pa_radians) + return ysh_rot, xsh_rot + + def sersic_map(self, mag = 12.0, r_eff = 2.0, n = 2.5, ellip = 0.5, theta = -50.0): + """ + Generate 2D map of Sersic model + + Parameters + ---------- + mag : float, optional + Integral magnitude, by default 12.0 + r_eff : float, optional + Effective radius in arcsec, by default 2.0 + n : float, optional + Sersic index, by default 2.5 + ellip : float, optional + Ellipcity, by default 0.5 + theta : float, optional + Position angle in degree, by default -50.0 + """ + self.mag = mag + self.reff = r_eff / self.xsamp + self.n = n + self.ellip = ellip + self.theta = theta + self.map = Sersic2D(self.x, self.y, mag = self.mag, + r_eff = self.reff, n = self.n, + ellip = self.ellip, theta = self.theta, + pixelscale = self.xsamp * self.ysamp) + + def tanh_map(self, vmax = 200.0, rt = 2.0, ellip = 0.5, theta = -50.0): + """ + Generate 2D velocity map of rotating disk according to tanh rotation curve + + Parameters + ---------- + vmax : float, optional + Maximum rotational velocity, by default 200.0km/s + rt : float, optional + Turn-over radius of rotation curve, by default 2.0 arcsec + ellip : float, optional + Apparent ellipcity of rotating disk, by default 0.5 + theta : float, optional + Position angle of rotating disk, by default -50.0 + """ + self.vmax = vmax + self.rt = rt / self.xsamp + self.ellip = ellip + self.theta = theta + self.map = VelMap2D(self.x, self.y, vmax = self.vmax, rt = self.rt, + ellip = self.ellip, theta = self.theta) + + def gred_map(self, a0 = 10, r_eff = 1, gred = -1, ellip = 0.5, theta = 0): + """ + Generate 2D maps according to the radial gradient form + + Parameters + ---------- + a0 : float, optional + Amplitude at the central pixel, by default 10 + r_eff : float, optional + Effective radius, by default 1 + gred : float, optional + Gradient of radial profile, by default -1 + ellip : float, optional + Ellipcity, by default 0.5 + theta : int, optional + Position angle, by default 0 + """ + self.a0 = a0 + self.reff = r_eff / self.xsamp + self.gred = gred + self.ellip = ellip + self.theta = theta + self.map = GradMap2D(self.x, self.y, a0 = self.a0, r_eff = self.reff, + gred = self.gred, ellip = self.ellip, theta = self.theta) + + def load_map(self, image): + """ + Generate 2D map according to input image + + Parameters + ---------- + image : 2d numpy array + The 2d array to be loaded. + """ + if np.ndim(image) == 2: + self.map = resize(image, (self.nx, self.ny)) + +class StellarPopulationMap(): + """ + Class of 2D maps for the parameters of stellar population, such as + surface brightness, median age and metallicity of stellar population, + velocity and velocity dispersion maps, and dust extinction. + + Parameters + ---------- + config : class + Class of configuration + sbright : class, optional + Class of the map of surface brightness of stellar population, by default None + logage : class, optional + Class of the map of stellar age, by default None + feh : class, optional + Class of the map of stellar metellicity, by default None + vel : class, optional + Class of the map of stellar velocity, by default None + vdisp : class, optional + Class of the map of stellar velocity dispersion, by default None + ebv : class, optional + Class of the map of dust extinction, by default None + """ + def __init__(self, config, sbright = None, logage = None, + feh = None, vel = None, vdisp = None, ebv = None): + + self.nx = config.nx + self.ny = config.ny + self.dpix = config.dpix + self.fov_x = config.fov_x + self.fov_y = config.fov_y + + self.sbright = sbright.map + self.logage = logage.map + self.feh = feh.map + self.vel = vel.map + self.vdisp = vdisp.map + self.ebv = ebv.map + + self.mag = self.sbright - 2.5 * np.log10(self.dpix * self.dpix) + self.age = 10 ** self.logage / 1e9 + + self.vdisp[self.vdisp < 10] = 10 + self.ebv[self.ebv < 0] = 0 + +class IonizedGasMap(): + """ + Class of 2D maps for the parameters of ionized gas, such as + Halpha flux map, gas-phase metallicity map, + velocity and velocity dispersion maps, and dust extinction. + + Parameters + ---------- + config : class + Class of configuration + halpha : class, optional + Class of the map of Halpha flux, by default None + zh : class, optional + Class of the map of gas-phase metallicity, by default None + vel : class, optional + Class of the map of gas velocity, by default None + vdisp : class, optional + Class of the map of gas velocity dispersion, by default None + ebv : class, optional + Class of the map of dust extinction, by default None + """ + def __init__(self, config, halpha = None, zh = None, vel = None, vdisp = None, ebv = None): + + self.nx = config.nx + self.ny = config.ny + self.dpix = config.dpix + self.fov_x = config.fov_x + self.fov_y = config.fov_y + + self.halpha = halpha.map + self.zh = zh.map + self.vel = vel.map + self.vdisp = vdisp.map + self.ebv = ebv.map + + self.vdisp[self.vdisp < 10] = 10 + self.ebv[self.ebv < 0] = 0 \ No newline at end of file diff --git a/src/gehong/spec1d.py b/src/gehong/spec1d.py new file mode 100644 index 0000000..9607e1a --- /dev/null +++ b/src/gehong/spec1d.py @@ -0,0 +1,1101 @@ +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.getenv('GEHONG_DATA_PATH') + +def readcol(filename, **kwargs): + """ + readcol, taken from ppxf. + + Parameters + ---------- + filename : string + The name of input ascii file + + Returns + ------- + float + The value of each columns + """ + 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. + This function is taken from ppxf. + + Parameters + ---------- + lamRange : array + Two elements vector containing the central wavelength + of the first and last pixels in the spectrum + spec : array + Input spectrum + oversample : bool, optional + Oversampling can be done, not to loose spectral resolution, + especally for extended wavelength ranges and to avoid aliasing, by default False + velscale : float, optional + velocity scale in km/s per pixels, by default None + flux : bool, optional + True to preserve total flux, by default False + + Returns + ------- + specNew : array + Output spectrum + logLam : array + Wavelength array in logarithm + velscale : array + velocity scale in km/s per pixels + """ + 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): + """ + One-dimensional Gaussian convolution + + Parameters + ---------- + spec : float array + vector with the spectrum to convolve + sig : float + vector of sigma values (in pixels) for every pixel + Returns + ------- + float array + Spectrum after convolution + """ + 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'): + """ + Flux calibration of spectrum + + Parameters + ---------- + wave : float array + Wavelength of input spectrum + flux : float array + Flux of input spectrum + mag : float + Magnitude used for flux calibration + filtername : str, optional + Filter band name, by default 'SLOAN_SDSS.r' + + Returns + ------- + float array + Spectrum after flux calibration + """ + # 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): + """ + Dust Extinction Curve of Calzetti et al. (2000) + + Parameters + ---------- + wave : float, or float array + 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 = (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 = (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 array + Wavelength of input spectra + flux : float array + 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 array + 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 array + Wavelength of spectrum + line_wave : float + Wavelength of emission line at the line center + FWHM_inst : float + Intrinsic broadening of emission line, units: A + + Returns + ------- + float array + 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 + + Parameters + ---------- + config : class + The class of configuration. + lam_range : list, optional + Wavelength range, by default [500, 15000] + dlam : float, optional + Wavelength width per pixel, by default 0.1A + model : str, optional + Emission line model, including 'hii' for HII region and 'nlr' for narrow line region of AGN, + by default 'hii' + """ + + def __init__(self, config, lam_range = [500, 15000], dlam = 0.1, model = 'hii'): + + self.lam_range = lam_range + self.wave = np.arange(lam_range[0], lam_range[1], 0.1) + self.FWHM_inst = config.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.logz_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.logz_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 + ---------- + config : class + Class of configuration + temp : class + Class of emission line template + halpha : float, optional + Integral flux of Halpha emission line, by default 100 * 1e-17 erg/s/cm^2 + logz : float, optional + Gas-phase metallicity, by default 0.0 + vel : float, optional + Line of sight velocity, by default 100.0km/s + vdisp : float, optional + Velocity dispersion, by default 120.0km/s + ebv : float, optional + Dust extinction, by default 0.1 + + Raises + ------ + ValueError + The value of logZ should be between -2 and 0.5. + """ + + def __init__(self, config, temp, halpha = 100.0, logz = 0.0, + vel = 100.0, vdisp = 120.0, ebv = 0.1): + + if (logz > -2) & (logz < 0.5): + indz = np.argmin(np.abs(logz - temp.logz_grid)) + flux_ratio = temp.flux_ratio[indz, :] + else: + raise ValueError('The value of logZ is not in the range of [-2, 0.5]!') + + # 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) + else: + flux_broad = flux_logwave + + # Redshift + redshift = vel / 3e5 + wave_r = np.exp(logLam) * (1 + redshift) + flux_red = np.interp(config.wave, wave_r, flux_broad) + + self.wave = config.wave + self.flux = flux_red + +############# +# AGN Spectra +############# + +class AGN_NLR(): + + """ + Class for narrow line region of AGN + + Parameters + ---------- + config : class + Class of configuration + temp : class + Class of emission line template + halpha : float, optional + Integral flux of Halpha emission line, by default 100 * 1e-17 erg/s/cm^2 + logz : float, optional + Gas-phase metallicity, by default 0.0 + vel : float, optional + Line of sight velocity, by default 100.0km/s + vdisp : float, optional + Velocity dispersion, by default 120.0km/s + ebv : float, optional + Dust extinction, by default 0.1 + + Raises + ------ + ValueError + The value of logZ should be between -2 and 0.5. + """ + + def __init__(self, config, temp, halpha = 100.0, logz = 0.0, + vel = 100.0, vdisp = 120.0, ebv = 0.1): + + if (logz > -2) & (logz < 0.5): + indz = np.argmin(np.abs(logz - temp.logz_grid)) + flux_ratio = temp.flux_ratio[indz, :] + else: + raise ValueError('The value of logZ is not in the range of [-2, 0.5]!') + + # Make emission line spectra through adding emission lines + emlines = temp.emission_lines * (flux_ratio / flux_ratio[6] ) + 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(config.wave, wave_r, flux_broad) + + self.wave = config.wave + self.flux = flux_red + +class AGN_BLR(): + + """ + Class for the broad line region of AGN + + Parameters + ---------- + config : class + Class of configuration + hbeta_flux : float, optional + Integral flux of Hbeta broad line, by default 100 * 1e-17 erg/s/cm^2 + hbeta_fwhm : float, optional + FWHM of Hbeta broad line, by default 2000.0km/s + vel : float, optional + Line of sight velocity, by default 100.0km/s + ebv : float, optional + Dust extinction, by default 0.1 + lam_range : list, optional + Wavelength range, by default [500, 15000] + """ + + def __init__(self, config, hbeta_flux = 100.0, hbeta_fwhm = 2000.0, ebv = 0.1, + 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(config.wave, wave_r, flux_dust) + + self.wave = config.wave + self.flux = flux_red + +class AGN_FeII(): + """ + Class for FeII emission lines of AGN + + Parameters + ---------- + config : class + Class of configuration + hbeta_broad : float, optional + Integral flux of Hbeta broad line, by default 100 * 1e-17 erg/s/cm^2 + r4570 : float, optional + Flux ratio between Fe4570 flux and Hbeta broad line flux, by default 0.4 + vel : float, optional + Line of sight velocity, by default 100.0km/s + ebv : float, optional + Dust extinction, by default 0.1 + """ + + def __init__(self, config, hbeta_broad = 100.0, r4570 = 0.4, ebv = 0.1, vel = 100.0): + + 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(config.wave, wave_r, flux_dust) + + self.wave = config.wave + self.flux = flux_red + +class AGN_Powerlaw(): + + """ + The class of power-law spectrum of AGN + + Parameters + ---------- + config : class + Class of configuration + M5100 : float, optional + Magnitude of power law spectrum between 5050A and 5150A, by default 1000.0 * 1e-17 erg/s/cm^2 + alpha : float, optional + Index of power law, by default -1.5 + vel : float, optional + Line of sight velocity, by default 100.0km/s + Ebv : float, optional + Dust extinction, by default 0.1 + """ + + def __init__(self, config, m5100 = 1000.0, alpha = -1.5, vel = 100.0, ebv = 0.1): + + 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(config.wave, wave_r, flux_dust) + + self.wave = config.wave + self.flux = flux_red + +class AGN(): + + """ + Class of singal spectra of AGN + + Parameters + ---------- + config : class + Class of configuration + nlr_template : class + Class of emission line template + bhmass : float, optional + Black hole mass used for calculating the luminosity of power law spectrum at 5100A, + by default 1e6 solar mass + edd_ratio : float, optional + Eddinton ratio used for calculating the luminosity of power law spectrum at 5100A, by default 0.05 + halpha_broad : float, optional + Integral flux of Halpha broad line, by default 100.0 * 1e-17 erg/s/cm^2 + halpha_narrow : float, optional + Integral flux of Halpha narrow line, by default 100.0 * 1e-17 erg/s/cm^2 + vdisp_broad : float, optional + Velocity dispersion of Halpha broad line, by default 5000.0km/s + vdisp_narrow : float, optional + Velocity dispersion of Halpha narrow line, by default 500.0km/s + vel : float, optional + Line of sight velocity, by default 1000.0km/s + logz : float, optional + Gas-phase metallicity of narrow line region, by default 0.0 + ebv : float, optional + Dust extinction, by default 0.1 + dist : float, optional + Luminosity distance of AGN, by default 20.0Mpc + """ + def __init__(self, config, nlr_template, bhmass = 1e6, edd_ratio = 0.05, + halpha_broad = 100.0, halpha_narrow = 100.0, vdisp_broad = 5000.0, vdisp_narrow = 500.0, + vel = 1000.0, logz = 0.0, ebv = 0.1, dist = 20.0): + + NLR = AGN_NLR(config, nlr_template, halpha = halpha_narrow, logz = logz, + vel = vel, vdisp = vdisp_narrow, ebv = ebv) + if halpha_broad > 0: + BLR = AGN_BLR(config, 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(config, m5100 = m5100, ebv = ebv, vel = vel) + Fe = AGN_FeII(config, hbeta_broad = halpha_broad / 2.579, ebv = ebv, vel = vel) + + self.wave = config.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.0): + """ + Caculate magnitude 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 : float, optional + Distance to the black hole, by default 21.0Mpc + + Returns + ------- + float + Magnitude 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 +################# + +""" + 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] + + """ + +def age_metal(filename): + """ + Extract the age and metallicity from the name of a file of + the MILES library of Single Stellar Population models. + + Parameters + ---------- + filename : string + Full path of template files + + Returns + ------- + age : float + Age of SSP (Gyr) + FeH : float + Metallicity of SSP + + Raises + ------ + ValueError + This is not a standard MILES filename + """ + 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): + """ + Class of single stellar population template. + + Parameters + ---------- + config : class + Class of configuration + velscale : array + velocity scale in km/s per pixels, by default 50.0km/s + pathname : string, optional + path with wildcards returning the list files to use, + by default data_path+'/data/EMILES/Ech*_baseFe.fits' + normalize : bool, optional + Set to True to normalize each template to mean=1, by default False + """ + def __init__(self, config, velscale = 50, + pathname = data_path + '/data/EMILES/Ech*_baseFe.fits', + normalize = False): + + FWHM_inst = config.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(): + """ + The class of stellar continuum + + Parameters + ---------- + config : class + Class of configuration + template : class + Class of single stellar population template + mag : float, optional + Magnitude in SDSS r-band, by default 15.0 + age : float, optional + Median age of stellar continuum, by default 1.0Gyr + feh : float, optional + Metallicity of stellar continuum, by default 0.0 + vel : float, optional + Line of sight velocity, by default 100.0km/s + vdisp : float, optional + Velocity dispersion, by default 120.0km/s + ebv : float, optional + Dust extinction, by default 0.1 + """ + def __init__(self, config, template, mag = 15.0, age = 1.0, feh = 0.0, + vel = 100.0, vdisp = 100.0, ebv = 0.1): + + # ----------------- + # 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(config.wave, wave_r, flux0) + + # Calibration + if np.isscalar(mag): + flux = calibrate(config.wave, flux, mag, filtername='SLOAN_SDSS.r') + + # Convert to input wavelength + self.wave = config.wave + self.flux = flux + +##################### +# Single Star Spectra +##################### + +class SingleStarTemplate(): + """ + Class of single stellar template + + Parameters + ---------- + config : class + Class of configuration + velscale : float, option + velocity scale in km/s per pixels, by default 20.0km/s + """ + def __init__(self, config, velscale = 20): + + FWHM_inst = config.inst_fwhm + filename = data_path + '/data/Starlib.XSL.fits' + + 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(): + + """ + Class of single stelar spectrum + + Parameters + ---------- + config : class + Class of configuration + template : class + Class of single stellar population template + mag : float, optional + Magnitude in SDSS r-band, by default 15.0 + Teff : float, optional + Effective tempreture, by default 10000.0K + FeH : float, optional + Metallicity of stellar, by default 0.0 + vel : float, optional + Line of sight velocity, by default 100.0km/s + Ebv : float, optional + Dust extinction, by default 0.1 + """ + + def __init__(self, config, template, mag = 15.0, teff = 10000.0, feh = 0.0, vel = 100.0, ebv = 0.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(config.wave, wave_r, starspec) + + # Calibration + if np.isscalar(mag): + flux = calibrate(config.wave, flux, mag, filtername='SLOAN_SDSS.r') + + # Convert to input wavelength + self.wave = config.wave + self.flux = flux \ No newline at end of file -- GitLab