Skip to content
GitLab
Commit 2c2dac03 authored by Fang Yuedong's avatar Fang Yuedong
Browse files

remove files

parent 4afd1181
Hide whitespace changes
Inline Side-by-side
ObservationSim/MockObject/CatalogBase.py deleted 100644 → 0
View file @ 4afd1181
import numpy as np
import galsim
import copy
import cmath
from astropy.table import Table
from abc import abstractmethod, ABCMeta
from ObservationSim.MockObject._util import integrate_sed_bandpass, getNormFactorForSpecWithABMAG, getABMAG
class CatalogBase(metaclass=ABCMeta):
def __init__(self):
pass
@abstractmethod
def load_sed(self, obj, **kward):
pass
@abstractmethod
def _load(self, **kward):
pass
@abstractmethod
def load_norm_filt(self, obj):
return None
@staticmethod
def initialize_param():
param = {
"star": -1,
"id": -1,
"ra": 0,
"dec": 0.,
"ra_orig": 0.,
"dec_orig": 0.,
"z": 0.,
"sed_type": -1,
"model_tag": "unknown",
"mag_use_normal": 100.,
"theta": 0.,
"kappa": 0.,
"g1": 0.,
"g2": 0.,
"bfrac": 0.,
"av": 0.,
"redden": 0.,
"hlr_bulge": 0.,
"hlr_disk": 0.,
"ell_bulge": 0.,
"ell_disk": 0.,
"ell_tot": 0.,
"e1_disk": 0.,
"e2_disk": 0.,
"e1_bulge": 0.,
"e2_bulge": 0.,
"teff": 0.,
"logg": 0.,
"feh": 0.,
"DM": 0.,
"stellarMass": 1.,
# C6 galaxies parameters
"e1": 0.,
"e2": 0.,
"bulgemass": 0.,
"diskmass": 0.,
"size": 0.,
"detA": 1.,
"type": 0,
"veldisp": 0.,
"coeff": np.zeros(20),
# Astrometry related
"pmra": 0.,
"pmdec": 0.,
"rv": 0.,
"parallax": 1e-9
}
return param
@staticmethod
def rotate_ellipticity(e1, e2, rotation=0., unit='radians'):
if unit == 'degree':
rotation = np.radians(rotation)
e_total = np.sqrt(e1**2 + e2**2)
phi = cmath.phase(complex(e1, e2))
e1 = e_total * np.cos(phi + 2*rotation)
e2 = e_total * np.sin(phi + 2*rotation)
return e1, e2, e_total
@staticmethod
def convert_sed(mag, sed, target_filt, norm_filt=None, mu=1.):
bandpass = target_filt.bandpass_full
if norm_filt is not None:
norm_thr_rang_ids = norm_filt['SENSITIVITY'] > 0.001
else:
norm_filt = Table(
np.array(np.array([bandpass.wave_list*10.0, bandpass.func(
bandpass.wave_list)])).T, names=(['WAVELENGTH', 'SENSITIVITY'])
)
norm_thr_rang_ids = norm_filt['SENSITIVITY'] > 0.001
sedNormFactor = getNormFactorForSpecWithABMAG(ABMag=mag,
spectrum=sed,
norm_thr=norm_filt,
sWave=np.floor(
norm_filt[norm_thr_rang_ids][0][0]),
eWave=np.ceil(norm_filt[norm_thr_rang_ids][-1][0]))
sed_photon = copy.copy(sed)
sed_photon = np.array(
[sed_photon['WAVELENGTH'], sed_photon['FLUX']*sedNormFactor]).T
sed_photon = galsim.LookupTable(x=np.array(sed_photon[:, 0]), f=np.array(
sed_photon[:, 1] * mu), interpolant='nearest')
# Get magnitude
sed_photon = galsim.SED(sed_photon, wave_type='A',
flux_type='1', fast=False)
interFlux = integrate_sed_bandpass(sed=sed_photon, bandpass=bandpass)
mag_csst = getABMAG(
interFlux=interFlux,
bandpass=bandpass
)
if target_filt.survey_type == "photometric":
return sed_photon, mag_csst, interFlux
elif target_filt.survey_type == "spectroscopic":
del sed_photon
return sed, mag_csst, interFlux
ObservationSim/MockObject/CosmicRay.py deleted 100755 → 0
View file @ 4afd1181
from ObservationSim.MockObject.MockObject import MockObject
class CosmicRay(MockObject):
pass
\ No newline at end of file
ObservationSim/MockObject/FlatLED.py deleted 100755 → 0
View file @ 4afd1181
import galsim
import os, sys
import numpy as np
import time
import math
import astropy.constants as cons
from astropy.io import fits
from scipy.interpolate import griddata
from astropy.table import Table
from ObservationSim.MockObject.SpecDisperser import SpecDisperser
from scipy import interpolate
import gc
from ObservationSim.MockObject.MockObject import MockObject
# from ObservationSim.Straylight import calculateSkyMap_split_g
try:
import importlib.resources as pkg_resources
except ImportError:
# Try backported to PY<37 'importlib_resources'
import importlib_resources as pkg_resources
# flatDir = '/Volumes/EAGET/LED_FLAT/'
LED_name = ['LED1', 'LED2', 'LED3', 'LED4', 'LED5', 'LED6', 'LED7', 'LED8', 'LED9', 'LED10', 'LED11', 'LED12', 'LED13',
'LED14']
cwaves_name = {'LED1': '275', 'LED2': '310', 'LED3': '430', 'LED4': '505', 'LED5': '545', 'LED6': '590', 'LED7': '670',
'LED8': '760', 'LED9': '940', 'LED10': '940', 'LED11': '1050', 'LED12': '1050',
'LED13': '340', 'LED14': '365'}
cwaves = {'LED1': 2750, 'LED2': 3100, 'LED3': 4300, 'LED4': 5050, 'LED5': 5250, 'LED6': 5900, 'LED7': 6700,
'LED8': 7600, 'LED9': 8800, 'LED10': 9400, 'LED11': 10500, 'LED12': 15500, 'LED13': 3400, 'LED14': 3650}
cwaves_fwhm = {'LED1': 110, 'LED2': 120, 'LED3': 200, 'LED4': 300, 'LED5': 300, 'LED6': 130, 'LED7': 210,
'LED8': 260, 'LED9': 400, 'LED10': 370, 'LED11': 500, 'LED12': 1400, 'LED13': 90, 'LED14': 100}
# LED_QE = {'LED1': 0.3, 'LED2': 0.4, 'LED13': 0.5, 'LED14': 0.5, 'LED10': 0.4}
# e-/ms
# fluxLED = {'LED1': 0.16478729, 'LED2': 0.084220931, 'LED3': 2.263360617, 'LED4': 2.190623489, 'LED5': 0.703504768,
# 'LED6': 0.446117963, 'LED7': 0.647122098, 'LED8': 0.922313442,
# 'LED9': 0.987278143, 'LED10': 2.043989167, 'LED11': 0.612571429, 'LED12': 1.228915663, 'LED13': 0.17029384,
# 'LED14': 0.27842925}
# e-/ms
fluxLED = {'LED1': 15, 'LED2': 15, 'LED3': 12.5, 'LED4': 9, 'LED5': 9,
'LED6': 9, 'LED7': 9, 'LED8': 9, 'LED9': 9, 'LED10': 12.5, 'LED11': 15, 'LED12':15, 'LED13': 12.5,
'LED14': 12.5}
# fluxLEDL = {'LED1': 10, 'LED2': 10, 'LED3': 10, 'LED4': 10, 'LED5': 10,
# 'LED6': 10, 'LED7': 10, 'LED8': 10, 'LED9': 10, 'LED10': 10, 'LED11': 10, 'LED12':10, 'LED13': 10,
# 'LED14': 10}
mirro_eff = {'GU':0.61, 'GV':0.8, 'GI':0.8}
# mirro_eff = {'GU':1, 'GV':1, 'GI':1}
class FlatLED(MockObject):
def __init__(self, chip,filt, flatDir = None, logger=None):
# self.led_type_list = led_type_list
self.filt = filt
self.chip = chip
self.logger = logger
if flatDir is not None:
self.flatDir = flatDir
else:
try:
with pkg_resources.files('ObservationSim.MockObject.data.led').joinpath("") as ledDir:
self.flatDir = ledDir.as_posix()
except AttributeError:
with pkg_resources.path('ObservationSim.MockObject.data.led', "") as ledDir:
self.flatDir = ledDir.as_posix()
###
### return LED flat, e/s
###
def getLEDImage(self, led_type='LED1'):
# cwave = cwaves[led_type]
flat = fits.open(os.path.join(self.flatDir, 'model_' + cwaves_name[led_type] + 'nm.fits'))
xlen = flat[0].header['NAXIS1']
ylen = 601
x = np.linspace(0, self.chip.npix_x * 6, xlen)
y = np.linspace(0, self.chip.npix_y * 5, ylen)
xx, yy = np.meshgrid(x, y)
a1 = flat[0].data[0:ylen, 0:xlen]
# z = np.sin((xx+yy+xx**2+yy**2))
# fInterp = interp2d(xx, yy, z, kind='linear')
X_ = np.hstack((xx.flatten()[:, None], yy.flatten()[:, None]))
Z_ = a1.flatten()
n_x = np.arange(0, self.chip.npix_x * 6, 1)
n_y = np.arange(0, self.chip.npix_y * 5, 1)
M, N = np.meshgrid(n_x, n_y)
i = self.chip.rowID - 1
j = self.chip.colID - 1
U = griddata(X_, Z_, (
M[self.chip.npix_y * i:self.chip.npix_y * (i + 1), self.chip.npix_x * j:self.chip.npix_x * (j + 1)],
N[self.chip.npix_y * i:self.chip.npix_y * (i + 1), self.chip.npix_x * j:self.chip.npix_x * (j + 1)]),
method='linear')
U = U/np.mean(U)
flatImage = U*fluxLED[led_type]*1000
gc.collect()
return flatImage
###
### return LED flat, e/s
###
def getLEDImage1(self, led_type='LED1'):
# cwave = cwaves[led_type]
flat = fits.open(os.path.join(self.flatDir, 'model_' + cwaves_name[led_type] + 'nm.fits'))
xlen = flat[0].header['NAXIS1']
ylen = 601
i = self.chip.rowID - 1
j = self.chip.colID - 1
x = np.linspace(0, self.chip.npix_x, int(xlen/6.))
y = np.linspace(0, self.chip.npix_y, int(ylen/5.))
xx, yy = np.meshgrid(x, y)
a1 = flat[0].data[int(ylen*i/5.):int(ylen*i/5.)+int(ylen/5.), int(xlen*j/6.):int(xlen*j/6.)+int(xlen/6.)]
# z = np.sin((xx+yy+xx**2+yy**2))
# fInterp = interp2d(xx, yy, z, kind='linear')
X_ = np.hstack((xx.flatten()[:, None], yy.flatten()[:, None]))
Z_ = a1.flatten()
n_x = np.arange(0, self.chip.npix_x , 1)
n_y = np.arange(0, self.chip.npix_y, 1)
M, N = np.meshgrid(n_x, n_y)
U = griddata(X_, Z_, (
M[0:self.chip.npix_y, 0:self.chip.npix_x],
N[0:self.chip.npix_y, 0:self.chip.npix_x ]),
method='linear')
U = U/np.mean(U)
flatImage = U*fluxLED[led_type]*1000
gc.collect()
return flatImage
def drawObj_LEDFlat_img(self, led_type_list=['LED1'], exp_t_list=[0.1]):
if len(led_type_list) > len(exp_t_list):
return np.ones([self.chip.npix_y,self.chip.npix_x])
ledFlat = np.zeros([self.chip.npix_y,self.chip.npix_x])
ledStat = '00000000000000'
ledTimes = [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0]
nledStat = '2'
for i in np.arange(len(led_type_list)):
led_type = led_type_list[i]
exp_t = exp_t_list[i]
# unitFlatImg = self.getLEDImage(led_type=led_type)
unitFlatImg = self.getLEDImage1(led_type=led_type)
# print("---------------TEST mem:",np.mean(unitFlatImg))
led_wave = cwaves[led_type]
led_fwhm = cwaves_fwhm[led_type]
led_spec = self.gaussian1d_profile_led(led_wave, led_fwhm)
speci = interpolate.interp1d(led_spec['WAVELENGTH'], led_spec['FLUX'])
w_list = np.arange(self.filt.blue_limit, self.filt.red_limit, 0.5) #A
f_spec = speci(w_list)
ccd_bp = self.chip._getChipEffCurve(self.chip.filter_type)
ccd_eff = ccd_bp.__call__(w_list / 10.)
filt_bp = self.filt.filter_bandpass
fil_eff = filt_bp.__call__(w_list / 10.)
t_spec = np.trapz(f_spec*ccd_eff*fil_eff, w_list)
# print(i, np.mean(unitFlatImg), t_spec, exp_t)
unitFlatImg = unitFlatImg * t_spec
# print("DEBUG1:---------------",np.mean(unitFlatImg))
ledFlat = ledFlat+unitFlatImg*exp_t
ledStat = ledStat[0:int(led_type[3:])-1]+nledStat+ledStat[int(led_type[3:]):]
ledTimes[int(led_type[3:])-1] = exp_t * 1000
gc.collect()
return ledFlat, ledStat, ledTimes
def drawObj_LEDFlat_slitless(self, led_type_list=['LED1'], exp_t_list=[0.1]):
if len(led_type_list) != len(exp_t_list):
return np.ones([self.chip.npix_y,self.chip.npix_x])
ledFlat = np.zeros([self.chip.npix_y,self.chip.npix_x])
ledStat = '00000000000000'
ledTimes = [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0]
nledStat = '2'
for i in np.arange(len(led_type_list)):
led_type = led_type_list[i]
exp_t = exp_t_list[i]
# unitFlatImg = self.getLEDImage(led_type=led_type)
unitFlatImg = self.getLEDImage1(led_type=led_type)
# print("---------------TEST mem:",np.mean(unitFlatImg))
ledFlat_ = unitFlatImg*exp_t
ledFlat_ = ledFlat_ / mirro_eff[self.filt.filter_type]
ledFlat_.astype(np.float32)
led_wave = cwaves[led_type]
led_fwhm = cwaves_fwhm[led_type]
led_spec = self.gaussian1d_profile_led(led_wave, led_fwhm)
# print("DEBUG1:---------------",np.mean(ledFlat_))
ledspec_map = self.calculateLEDSpec(
skyMap=ledFlat_,
blueLimit=self.filt.blue_limit,
redLimit=self.filt.red_limit,
conf=self.chip.sls_conf,
pixelSize=self.chip.pix_scale,
isAlongY=0,
flat_cube=self.chip.flat_cube, led_spec=led_spec)
ledFlat = ledFlat + ledspec_map
ledStat = ledStat[0:int(led_type[3:])-1]+nledStat+ledStat[int(led_type[3:]):]
ledTimes[int(led_type[3:])-1] = exp_t * 1000
return ledFlat, ledStat, ledTimes
def drawObj_LEDFlat(self, led_type_list=['LED1'], exp_t_list=[0.1]):
if self.chip.survey_type == "photometric":
return self.drawObj_LEDFlat_img(led_type_list=led_type_list, exp_t_list=exp_t_list)
elif self.chip.survey_type == "spectroscopic":
return self.drawObj_LEDFlat_slitless(led_type_list=led_type_list, exp_t_list=exp_t_list)
def gaussian1d_profile_led(self, xc=5050, fwhm=300):
sigma = fwhm/2.355
x_radii = int(5*sigma + 1)
xlist = np.arange(xc-x_radii, xc+x_radii, 0.5)
xlist_ = np.zeros(len(xlist) + 2)
xlist_[1:-1] = xlist
xlist_[0] = 2000
xlist_[-1] = 18000
ids1 = xlist>xc-fwhm
ids2 = xlist[ids1]<xc+fwhm
data = np.exp((-(xlist-xc)*(xlist-xc))/(2*sigma*sigma))/(np.sqrt(2*math.pi)*sigma)
scale = 1/np.trapz(data[ids1][ids2], xlist[ids1][ids2])
data_ = np.zeros(len(xlist) + 2)
data_[1:-1] = data*scale
# print("DEBUG:-------------------------------",np.sum(data_), scale)
return Table(np.array([xlist_.astype(np.float32), data_.astype(np.float32)]).T, names=('WAVELENGTH', 'FLUX'))
def calculateLEDSpec(self, skyMap=None, blueLimit=4200, redLimit=6500,
conf=[''], pixelSize=0.074, isAlongY=0,
split_pos=3685, flat_cube=None, led_spec=None):
conf1 = conf[0]
conf2 = conf[0]
if np.size(conf) == 2:
conf2 = conf[1]
skyImg = galsim.Image(skyMap, xmin=0, ymin=0)
tbstart = blueLimit
tbend = redLimit
fimg = np.zeros_like(skyMap)
fImg = galsim.Image(fimg)
spec = led_spec
if isAlongY == 0:
directParm = 0
if isAlongY == 1:
directParm = 1
if split_pos >= skyImg.array.shape[directParm]:
skyImg1 = galsim.Image(skyImg.array)
origin1 = [0, 0]
# sdp = specDisperser.specDisperser(orig_img=skyImg1, xcenter=skyImg1.center.x, ycenter=skyImg1.center.y,
# full_img=fimg, tar_spec=spec, band_start=tbstart, band_end=tbend,
# origin=origin1,
# conf=conf1)
# sdp.compute_spec_orders()
y_len = skyMap.shape[0]
x_len = skyMap.shape[1]
delt_x = 100
delt_y = 100
sub_y_start_arr = np.arange(0, y_len, delt_y)
sub_y_end_arr = sub_y_start_arr + delt_y
sub_y_end_arr[-1] = min(sub_y_end_arr[-1], y_len)
sub_x_start_arr = np.arange(0, x_len, delt_x)
sub_x_end_arr = sub_x_start_arr + delt_x
sub_x_end_arr[-1] = min(sub_x_end_arr[-1], x_len)
for i, k1 in enumerate(sub_y_start_arr):
sub_y_s = k1
sub_y_e = sub_y_end_arr[i]
sub_y_center = (sub_y_s + sub_y_e) / 2.
for j, k2 in enumerate(sub_x_start_arr):
sub_x_s = k2
sub_x_e = sub_x_end_arr[j]
skyImg_sub = galsim.Image(skyImg.array[sub_y_s:sub_y_e, sub_x_s:sub_x_e])
origin_sub = [sub_y_s, sub_x_s]
sub_x_center = (sub_x_s + sub_x_e) / 2.
sdp = SpecDisperser(orig_img=skyImg_sub, xcenter=sub_x_center, ycenter=sub_y_center,
origin=origin_sub,
tar_spec=spec,
band_start=tbstart, band_end=tbend,
conf=conf2,
flat_cube=flat_cube)
spec_orders = sdp.compute_spec_orders()
for k, v in spec_orders.items():
img_s = v[0]
origin_order_x = v[1]
origin_order_y = v[2]
ssImg = galsim.ImageF(img_s)
ssImg.setOrigin(origin_order_x, origin_order_y)
bounds = ssImg.bounds & fImg.bounds
if bounds.area() == 0:
continue
fImg[bounds] = fImg[bounds] + ssImg[bounds]
else:
# sdp.compute_spec_orders()
y_len = skyMap.shape[0]
x_len = skyMap.shape[1]
delt_x = 500
delt_y = y_len
sub_y_start_arr = np.arange(0, y_len, delt_y)
sub_y_end_arr = sub_y_start_arr + delt_y
sub_y_end_arr[-1] = min(sub_y_end_arr[-1], y_len)
delt_x = split_pos - 0
sub_x_start_arr = np.arange(0, split_pos, delt_x)
sub_x_end_arr = sub_x_start_arr + delt_x
sub_x_end_arr[-1] = min(sub_x_end_arr[-1], split_pos)
for i, k1 in enumerate(sub_y_start_arr):
sub_y_s = k1
sub_y_e = sub_y_end_arr[i]
sub_y_center = (sub_y_s + sub_y_e) / 2.
for j, k2 in enumerate(sub_x_start_arr):
sub_x_s = k2
sub_x_e = sub_x_end_arr[j]
# print(i,j,sub_y_s, sub_y_e,sub_x_s,sub_x_e)
T1 = time.time()
skyImg_sub = galsim.Image(skyImg.array[sub_y_s:sub_y_e, sub_x_s:sub_x_e])
origin_sub = [sub_y_s, sub_x_s]
sub_x_center = (sub_x_s + sub_x_e) / 2.
sdp = SpecDisperser(orig_img=skyImg_sub, xcenter=sub_x_center, ycenter=sub_y_center,
origin=origin_sub,
tar_spec=spec,
band_start=tbstart, band_end=tbend,
conf=conf1,
flat_cube=flat_cube)
spec_orders = sdp.compute_spec_orders()
for k, v in spec_orders.items():
img_s = v[0]
origin_order_x = v[1]
origin_order_y = v[2]
ssImg = galsim.ImageF(img_s)
ssImg.setOrigin(origin_order_x, origin_order_y)
bounds = ssImg.bounds & fImg.bounds
if bounds.area() == 0:
continue
fImg[bounds] = fImg[bounds] + ssImg[bounds]
T2 = time.time()
print('time: %s ms' % ((T2 - T1) * 1000))
delt_x = x_len - split_pos
sub_x_start_arr = np.arange(split_pos, x_len, delt_x)
sub_x_end_arr = sub_x_start_arr + delt_x
sub_x_end_arr[-1] = min(sub_x_end_arr[-1], x_len)
for i, k1 in enumerate(sub_y_start_arr):
sub_y_s = k1
sub_y_e = sub_y_end_arr[i]
sub_y_center = (sub_y_s + sub_y_e) / 2.
for j, k2 in enumerate(sub_x_start_arr):
sub_x_s = k2
sub_x_e = sub_x_end_arr[j]
# print(i,j,sub_y_s, sub_y_e,sub_x_s,sub_x_e)
T1 = time.time()
skyImg_sub = galsim.Image(skyImg.array[sub_y_s:sub_y_e, sub_x_s:sub_x_e])
origin_sub = [sub_y_s, sub_x_s]
sub_x_center = (sub_x_s + sub_x_e) / 2.
sdp = SpecDisperser(orig_img=skyImg_sub, xcenter=sub_x_center, ycenter=sub_y_center,
origin=origin_sub,
tar_spec=spec,
band_start=tbstart, band_end=tbend,
conf=conf2,
flat_cube=flat_cube)
spec_orders = sdp.compute_spec_orders()
for k, v in spec_orders.items():
img_s = v[0]
origin_order_x = v[1]
origin_order_y = v[2]
ssImg = galsim.ImageF(img_s)
ssImg.setOrigin(origin_order_x, origin_order_y)
bounds = ssImg.bounds & fImg.bounds
if bounds.area() == 0:
continue
fImg[bounds] = fImg[bounds] + ssImg[bounds]
T2 = time.time()
print('time: %s ms' % ((T2 - T1) * 1000))
if isAlongY == 1:
fimg, tmx, tmy = SpecDisperser.rotate90(array_orig=fImg.array, xc=0, yc=0, isClockwise=0)
else:
fimg = fImg.array
# fimg = fimg * pixelSize * pixelSize
return fimg
ObservationSim/MockObject/Galaxy.py deleted 100755 → 0
View file @ 4afd1181
import numpy as np
import galsim
from astropy.table import Table
from ObservationSim.MockObject._util import eObs, integrate_sed_bandpass, getNormFactorForSpecWithABMAG
from ObservationSim.MockObject.SpecDisperser import SpecDisperser
from ObservationSim.MockObject.MockObject import MockObject
# import tracemalloc
class Galaxy(MockObject):
def __init__(self, param, logger=None):
super().__init__(param, logger=logger)
if not hasattr(self, "disk_sersic_idx"):
self.disk_sersic_idx = 1.
if not hasattr(self, "bulge_sersic_idx"):
self.bulge_sersic_idx = 4.
if not hasattr(self, "mu"):
if hasattr(self, "detA"):
self.mu = 1./self.detA
else:
self.mu = 1.
def unload_SED(self):
"""(Test) free up SED memory
"""
del self.sed
def getGSObj_multiband(self, tel, psf_list, bandpass_list, filt, nphotons_tot=None, g1=0, g2=0, exptime=150., fd_shear=None):
if len(psf_list) != len(bandpass_list):
raise ValueError(
"!!!The number of PSF profiles and the number of bandpasses must be equal.")
objs = []
if nphotons_tot == None:
nphotons_tot = self.getElectronFluxFilt(filt, tel, exptime)
# print("nphotons_tot = ", nphotons_tot)
try:
full = integrate_sed_bandpass(
sed=self.sed, bandpass=filt.bandpass_full)
except Exception as e:
print(e)
if self.logger:
self.logger.error(e)
return -1
for i in range(len(bandpass_list)):
bandpass = bandpass_list[i]
try:
sub = integrate_sed_bandpass(sed=self.sed, bandpass=bandpass)
except Exception as e:
print(e)
if self.logger:
self.logger.error(e)
return -1
ratio = sub/full
if not (ratio == -1 or (ratio != ratio)):
nphotons = ratio * nphotons_tot
else:
return -1
psf = psf_list[i]
disk = galsim.Sersic(n=self.disk_sersic_idx,
half_light_radius=self.hlr_disk, flux=1.0)
disk_shape = galsim.Shear(g1=self.e1_disk, g2=self.e2_disk)
disk = disk.shear(disk_shape)
bulge = galsim.Sersic(n=self.bulge_sersic_idx,
half_light_radius=self.hlr_bulge, flux=1.0)
bulge_shape = galsim.Shear(g1=self.e1_bulge, g2=self.e2_bulge)
bulge = bulge.shear(bulge_shape)
if self.bfrac == 0:
gal = disk
elif self.bfrac == 1:
gal = bulge
else:
gal = self.bfrac * bulge + (1.0 - self.bfrac) * disk
if fd_shear is not None:
g1 += fd_shear.g1
g2 += fd_shear.g2
gal_shear = galsim.Shear(g1=g1, g2=g2)
gal = gal.shear(gal_shear)
# Magnification
gal = gal.magnify(self.mu)
gal = galsim.Convolve(psf, gal)
gal = gal.withFlux(nphotons)
objs.append(gal)
final = galsim.Sum(objs)
return final
def drawObj_multiband(self, tel, pos_img, psf_model, bandpass_list, filt, chip, nphotons_tot=None, g1=0, g2=0, exptime=150., fd_shear=None):
if nphotons_tot == None:
nphotons_tot = self.getElectronFluxFilt(filt, tel, exptime)
# print("nphotons_tot = ", nphotons_tot)
try:
full = integrate_sed_bandpass(
sed=self.sed, bandpass=filt.bandpass_full)
except Exception as e:
print(e)
if self.logger:
self.logger.error(e)
return 2, None
# # [C6 TEST]
# print('hlr_disk = %.4f, hlr_bulge = %.4f'%(self.hlr_disk, self.hlr_bulge))
# tracemalloc.start()
big_galaxy = False
if self.hlr_disk > 3.0 or self.hlr_bulge > 3.0: # Very big galaxy
big_galaxy = True
# Set Galsim Parameters
if self.getMagFilter(filt) <= 15 and (not big_galaxy):
folding_threshold = 5.e-4
else:
folding_threshold = 5.e-3
gsp = galsim.GSParams(folding_threshold=folding_threshold)
self.real_pos = self.getRealPos(chip.img, global_x=self.posImg.x, global_y=self.posImg.y,
img_real_wcs=self.chip_wcs)
x, y = self.real_pos.x + 0.5, self.real_pos.y + 0.5
x_nominal = int(np.floor(x + 0.5))
y_nominal = int(np.floor(y + 0.5))
dx = x - x_nominal
dy = y - y_nominal
offset = galsim.PositionD(dx, dy)
# Get real local wcs of object (deal with chip rotation w.r.t its center)
chip_wcs_local = self.chip_wcs.local(self.real_pos)
is_updated = 0
# Model the galaxy as disk + bulge
disk = galsim.Sersic(
n=self.disk_sersic_idx, half_light_radius=self.hlr_disk, flux=1.0, gsparams=gsp)
disk_shape = galsim.Shear(g1=self.e1_disk, g2=self.e2_disk)
disk = disk.shear(disk_shape)
bulge = galsim.Sersic(
n=self.bulge_sersic_idx, half_light_radius=self.hlr_bulge, flux=1.0, gsparams=gsp)
bulge_shape = galsim.Shear(g1=self.e1_bulge, g2=self.e2_bulge)
bulge = bulge.shear(bulge_shape)
# Get shear and deal with shear induced by field distortion
if fd_shear:
g1 += fd_shear.g1
g2 += fd_shear.g2
gal_shear = galsim.Shear(g1=g1, g2=g2)
# Loop over all sub-bandpasses
for i in range(len(bandpass_list)):
bandpass = bandpass_list[i]
try:
sub = integrate_sed_bandpass(sed=self.sed, bandpass=bandpass)
except Exception as e:
print(e)
if self.logger:
self.logger.error(e)
continue
ratio = sub/full
if not (ratio == -1 or (ratio != ratio)):
nphotons = ratio * nphotons_tot
else:
continue
# nphotons_sum += nphotons
# # [C6 TEST]
# print("nphotons_sub-band_%d = %.2f"%(i, nphotons))
# Get PSF model
psf, pos_shear = psf_model.get_PSF(
chip=chip, pos_img=pos_img, bandpass=bandpass, folding_threshold=folding_threshold)
if self.bfrac == 0:
gal_temp = disk
elif self.bfrac == 1:
gal_temp = bulge
else:
gal_temp = self.bfrac * bulge + (1.0 - self.bfrac) * disk
gal_temp = gal_temp.shear(gal_shear)
# Magnification
gal_temp = gal_temp.magnify(self.mu)
if not big_galaxy: # Not apply PSF for very big galaxy
gal_temp = galsim.Convolve(psf, gal_temp)
gal_temp = gal_temp.withFlux(nphotons)
if i == 0:
gal = gal_temp
else:
gal = gal + gal_temp
# (TEST) Random knots
# knots = galsim.RandomKnots(npoints=100, profile=disk)
# kfrac = np.random.random()*(1.0 - self.bfrac)
# gal = self.bfrac * bulge + (1.0 - self.bfrac - kfrac) * disk + kfrac * knots
# stamp = gal.drawImage(wcs=chip_wcs_local, method='phot', offset=offset, save_photons=True)
stamp = gal.drawImage(wcs=chip_wcs_local, offset=offset)
if np.sum(np.isnan(stamp.array)) > 0:
# ERROR happens
return 2, pos_shear
stamp.setCenter(x_nominal, y_nominal)
bounds = stamp.bounds & galsim.BoundsI(
0, chip.npix_x - 1, 0, chip.npix_y - 1)
if bounds.area() > 0:
chip.img.setOrigin(0, 0)
chip.img[bounds] += stamp[bounds]
is_updated = 1
chip.img.setOrigin(chip.bound.xmin, chip.bound.ymin)
del stamp
if is_updated == 0:
# Return code 0: object photons missed this detector
print("obj %s missed" % (self.id))
if self.logger:
self.logger.info("obj %s missed" % (self.id))
return 0, pos_shear
# # [C6 TEST]
# print("nphotons_sum = ", nphotons_sum)
return 1, pos_shear
def drawObj_slitless(self, tel, pos_img, psf_model, bandpass_list, filt, chip, nphotons_tot=None, g1=0, g2=0,
exptime=150., normFilter=None, grating_split_pos=3685, fd_shear=None):
if normFilter is not None:
norm_thr_rang_ids = normFilter['SENSITIVITY'] > 0.001
sedNormFactor = getNormFactorForSpecWithABMAG(ABMag=self.param['mag_use_normal'], spectrum=self.sed,
norm_thr=normFilter,
sWave=np.floor(
normFilter[norm_thr_rang_ids][0][0]),
eWave=np.ceil(normFilter[norm_thr_rang_ids][-1][0]))
if sedNormFactor == 0:
return 2, None
else:
sedNormFactor = 1.
normalSED = Table(np.array([self.sed['WAVELENGTH'], self.sed['FLUX'] * sedNormFactor]).T,
names=('WAVELENGTH', 'FLUX'))
self.real_pos = self.getRealPos(chip.img, global_x=self.posImg.x, global_y=self.posImg.y,
img_real_wcs=self.chip_wcs)
x, y = self.real_pos.x + 0.5, self.real_pos.y + 0.5
x_nominal = int(np.floor(x + 0.5))
y_nominal = int(np.floor(y + 0.5))
dx = x - x_nominal
dy = y - y_nominal
offset = galsim.PositionD(dx, dy)
chip_wcs_local = self.chip_wcs.local(self.real_pos)
big_galaxy = False
if self.hlr_disk > 3.0 or self.hlr_bulge > 3.0: # Very big galaxy
big_galaxy = True
if self.getMagFilter(filt) <= 15 and (not big_galaxy):
folding_threshold = 5.e-4
else:
folding_threshold = 5.e-3
gsp = galsim.GSParams(folding_threshold=folding_threshold)
# nphotons_sum = 0
flat_cube = chip.flat_cube
xOrderSigPlus = {'A': 1.3909419820029296, 'B': 1.4760376591236062,
'C': 4.035447379743442, 'D': 5.5684364343742825, 'E': 16.260021029735388}
grating_split_pos_chip = 0 + grating_split_pos
branges = np.zeros([len(bandpass_list), 2])
# print(hasattr(psf_model, 'bandranges'))
if hasattr(psf_model, 'bandranges'):
if psf_model.bandranges is None:
return 2, None
if len(psf_model.bandranges) != len(bandpass_list):
return 2, None
branges = psf_model.bandranges
else:
for i in range(len(bandpass_list)):
branges[i, 0] = bandpass_list[i].blue_limit * 10
branges[i, 1] = bandpass_list[i].red_limit * 10
for i in range(len(bandpass_list)):
# bandpass = bandpass_list[i]
brange = branges[i]
# psf, pos_shear = psf_model.get_PSF(chip=chip, pos_img=pos_img, bandpass=bandpass, folding_threshold=folding_threshold)
disk = galsim.Sersic(
n=self.disk_sersic_idx, half_light_radius=self.hlr_disk, flux=1.0, gsparams=gsp)
disk_shape = galsim.Shear(g1=self.e1_disk, g2=self.e2_disk)
disk = disk.shear(disk_shape)
bulge = galsim.Sersic(
n=self.bulge_sersic_idx, half_light_radius=self.hlr_bulge, flux=1.0, gsparams=gsp)
bulge_shape = galsim.Shear(g1=self.e1_bulge, g2=self.e2_bulge)
bulge = bulge.shear(bulge_shape)
if self.bfrac == 0:
gal = disk
elif self.bfrac == 1:
gal = bulge
else:
gal = self.bfrac * bulge + (1.0 - self.bfrac) * disk
# (TEST) Random knots
# knots = galsim.RandomKnots(npoints=100, profile=disk)
# kfrac = np.random.random()*(1.0 - self.bfrac)
# gal = self.bfrac * bulge + (1.0 - self.bfrac - kfrac) * disk + kfrac * knots
if fd_shear:
g1 += fd_shear.g1
g2 += fd_shear.g2
gal_shear = galsim.Shear(g1=g1, g2=g2)
gal = gal.shear(gal_shear)
gal = gal.magnify(self.mu)
gal = gal.withFlux(tel.pupil_area * exptime)
# gal = galsim.Convolve(psf, gal)
# if not big_galaxy: # Not apply PSF for very big galaxy
# gal = galsim.Convolve(psf, gal)
# # if fd_shear is not None:
# # gal = gal.shear(fd_shear)
starImg = gal.drawImage(
wcs=chip_wcs_local, offset=offset, method='real_space')
origin_star = [y_nominal - (starImg.center.y - starImg.ymin),
x_nominal - (starImg.center.x - starImg.xmin)]
starImg.setOrigin(0, 0)
gal_origin = [origin_star[0], origin_star[1]]
gal_end = [origin_star[0] + starImg.array.shape[0] -
1, origin_star[1] + starImg.array.shape[1] - 1]
if gal_origin[1] < grating_split_pos_chip < gal_end[1]:
subSlitPos = int(grating_split_pos_chip - gal_origin[1] + 1)
# part img disperse
subImg_p1 = starImg.array[:, 0:subSlitPos]
star_p1 = galsim.Image(subImg_p1)
star_p1.setOrigin(0, 0)
origin_p1 = origin_star
xcenter_p1 = min(x_nominal, grating_split_pos_chip-1) - 0
ycenter_p1 = y_nominal-0
sdp_p1 = SpecDisperser(orig_img=star_p1, xcenter=xcenter_p1,
ycenter=ycenter_p1, origin=origin_p1,
tar_spec=normalSED,
band_start=brange[0], band_end=brange[1],
conf=chip.sls_conf[0],
isAlongY=0,
flat_cube=flat_cube)
# self.addSLStoChipImage(sdp=sdp_p1, chip=chip, xOrderSigPlus = xOrderSigPlus, local_wcs=chip_wcs_local)
pos_shear = self.addSLStoChipImageWithPSF(sdp=sdp_p1, chip=chip, pos_img_local=[xcenter_p1, ycenter_p1],
psf_model=psf_model, bandNo=i + 1,
grating_split_pos=grating_split_pos,
local_wcs=chip_wcs_local, pos_img=pos_img)
subImg_p2 = starImg.array[:,
subSlitPos+1:starImg.array.shape[1]]
star_p2 = galsim.Image(subImg_p2)
star_p2.setOrigin(0, 0)
origin_p2 = [origin_star[0], grating_split_pos_chip]
xcenter_p2 = max(x_nominal, grating_split_pos_chip - 1) - 0
ycenter_p2 = y_nominal - 0
sdp_p2 = SpecDisperser(orig_img=star_p2, xcenter=xcenter_p2,
ycenter=ycenter_p2, origin=origin_p2,
tar_spec=normalSED,
band_start=brange[0], band_end=brange[1],
conf=chip.sls_conf[1],
isAlongY=0,
flat_cube=flat_cube)
# self.addSLStoChipImage(sdp=sdp_p2, chip=chip, xOrderSigPlus = xOrderSigPlus, local_wcs=chip_wcs_local)
pos_shear = self.addSLStoChipImageWithPSF(sdp=sdp_p2, chip=chip, pos_img_local=[xcenter_p2, ycenter_p2],
psf_model=psf_model, bandNo=i + 1,
grating_split_pos=grating_split_pos,
local_wcs=chip_wcs_local, pos_img=pos_img)
del sdp_p1
del sdp_p2
elif grating_split_pos_chip <= gal_origin[1]:
sdp = SpecDisperser(orig_img=starImg, xcenter=x_nominal - 0,
ycenter=y_nominal - 0, origin=origin_star,
tar_spec=normalSED,
band_start=brange[0], band_end=brange[1],
conf=chip.sls_conf[1],
isAlongY=0,
flat_cube=flat_cube)
# self.addSLStoChipImage(sdp=sdp, chip=chip, xOrderSigPlus = xOrderSigPlus, local_wcs=chip_wcs_local)
pos_shear = self.addSLStoChipImageWithPSF(sdp=sdp, chip=chip, pos_img_local=[x_nominal, y_nominal],
psf_model=psf_model, bandNo=i + 1,
grating_split_pos=grating_split_pos,
local_wcs=chip_wcs_local, pos_img=pos_img)
del sdp
elif grating_split_pos_chip >= gal_end[1]:
sdp = SpecDisperser(orig_img=starImg, xcenter=x_nominal - 0,
ycenter=y_nominal - 0, origin=origin_star,
tar_spec=normalSED,
band_start=brange[0], band_end=brange[1],
conf=chip.sls_conf[0],
isAlongY=0,
flat_cube=flat_cube)
# self.addSLStoChipImage(sdp=sdp, chip=chip, xOrderSigPlus = xOrderSigPlus, local_wcs=chip_wcs_local)
pos_shear = self.addSLStoChipImageWithPSF(sdp=sdp, chip=chip, pos_img_local=[x_nominal, y_nominal],
psf_model=psf_model, bandNo=i + 1,
grating_split_pos=grating_split_pos,
local_wcs=chip_wcs_local, pos_img=pos_img)
del sdp
# print(self.y_nominal, starImg.center.y, starImg.ymin)
# del psf
return 1, pos_shear
def getGSObj(self, psf, g1=0, g2=0, flux=None, filt=None, tel=None, exptime=150.):
if flux == None:
flux = self.getElectronFluxFilt(filt, tel, exptime)
disk = galsim.Sersic(n=self.disk_sersic_idx,
half_light_radius=self.hlr_disk, flux=1.0)
disk_shape = galsim.Shear(g1=self.e1_disk, g2=self.e2_disk)
disk = disk.shear(disk_shape)
bulge = galsim.Sersic(n=self.bulge_sersic_idx,
half_light_radius=self.hlr_bulge, flux=1.0)
bulge_shape = galsim.Shear(g1=self.e1_bulge, g2=self.e2_bulge)
bulge = bulge.shear(bulge_shape)
gal = self.bfrac * bulge + (1.0 - self.bfrac) * disk
gal = gal.withFlux(flux)
gal_shear = galsim.Shear(g1=g1, g2=g2)
gal = gal.shear(gal_shear)
final = galsim.Convolve(psf, gal)
return final
def getObservedEll(self, g1=0, g2=0):
e1_obs, e2_obs, e_obs, theta = eObs(
self.e1_total, self.e2_total, g1, g2)
return self.e1_total, self.e2_total, g1, g2, e1_obs, e2_obs
ObservationSim/MockObject/MockObject.py deleted 100755 → 0
View file @ 4afd1181
import os
import galsim
import numpy as np
import astropy.constants as cons
from astropy import wcs
from astropy.table import Table
import astropy.io.fits as fitsio
from ObservationSim.MockObject._util import magToFlux, VC_A, convolveGaussXorders, convolveImg
from ObservationSim.MockObject._util import integrate_sed_bandpass, getNormFactorForSpecWithABMAG, getObservedSED, \
getABMAG
from ObservationSim.MockObject.SpecDisperser import SpecDisperser
class MockObject(object):
def __init__(self, param, logger=None):
self.param = param
for key in self.param:
setattr(self, key, self.param[key])
if self.param["star"] == 0:
self.type = "galaxy"
elif self.param["star"] == 1:
self.type = "star"
elif self.param["star"] == 2:
self.type = "quasar"
###mock_stamp_START
elif self.param["star"] == 3:
self.type = "stamp"
###mock_stamp_END
###for calibration
elif self.param["star"] == 4:
self.type = "calib"
###END
self.sed = None
self.fd_shear = None
# Place holder for outputs
self.additional_output_str = ""
self.logger = logger
def getMagFilter(self, filt):
if filt.filter_type in ["GI", "GV", "GU"]:
return self.param["mag_use_normal"]
return self.param["mag_%s" % filt.filter_type.lower()]
def getFluxFilter(self, filt):
return self.param["flux_%s" % filt.filter_type.lower()]
def getNumPhotons(self, flux, tel, exptime=150.):
pupil_area = tel.pupil_area * (100.) ** 2 # m^2 to cm^2
return flux * pupil_area * exptime
def getElectronFluxFilt(self, filt, tel, exptime=150.):
# photonEnergy = filt.getPhotonE()
# flux = magToFlux(self.getMagFilter(filt))
# factor = 1.0e4 * flux/photonEnergy * VC_A * (1.0/filt.blue_limit - 1.0/filt.red_limit)
# return factor * filt.efficiency * tel.pupil_area * exptime
flux = self.getFluxFilter(filt)
return flux * tel.pupil_area * exptime
def getPosWorld(self):
ra = self.param["ra"]
dec = self.param["dec"]
return galsim.CelestialCoord(ra=ra * galsim.degrees, dec=dec * galsim.degrees)
def getPosImg_Offset_WCS(self, img, fdmodel=None, chip=None, verbose=True, chip_wcs=None, img_header=None):
self.posImg = img.wcs.toImage(self.getPosWorld())
self.localWCS = img.wcs.local(self.posImg)
# Apply field distortion model
if (fdmodel is not None) and (chip is not None):
if verbose:
print("\n")
print("Before field distortion:\n")
print("x = %.2f, y = %.2f\n" % (self.posImg.x, self.posImg.y), flush=True)
self.posImg, self.fd_shear = fdmodel.get_distorted(chip=chip, pos_img=self.posImg)
if verbose:
print("After field distortion:\n")
print("x = %.2f, y = %.2f\n" % (self.posImg.x, self.posImg.y), flush=True)
x, y = self.posImg.x + 0.5, self.posImg.y + 0.5
self.x_nominal = int(np.floor(x + 0.5))
self.y_nominal = int(np.floor(y + 0.5))
dx = x - self.x_nominal
dy = y - self.y_nominal
self.offset = galsim.PositionD(dx, dy)
# Deal with chip rotation
if chip_wcs is not None:
self.chip_wcs = chip_wcs
elif img_header is not None:
self.chip_wcs = galsim.FitsWCS(header=img_header)
else:
self.chip_wcs = None
return self.posImg, self.offset, self.localWCS, self.chip_wcs, self.fd_shear
def getRealPos(self, img, global_x=0., global_y=0., img_real_wcs=None):
img_global_pos = galsim.PositionD(global_x, global_y)
cel_pos = img.wcs.toWorld(img_global_pos)
realPos = img_real_wcs.toImage(cel_pos)
return realPos
def drawObj_multiband(self, tel, pos_img, psf_model, bandpass_list, filt, chip, nphotons_tot=None, g1=0, g2=0,
exptime=150., fd_shear=None):
if nphotons_tot == None:
nphotons_tot = self.getElectronFluxFilt(filt, tel, exptime)
# print("nphotons_tot = ", nphotons_tot)
try:
full = integrate_sed_bandpass(sed=self.sed, bandpass=filt.bandpass_full)
except Exception as e:
print(e)
if self.logger:
self.logger.error(e)
return 2, None
# Set Galsim Parameters
if self.getMagFilter(filt) <= 15:
folding_threshold = 5.e-4
else:
folding_threshold = 5.e-3
gsp = galsim.GSParams(folding_threshold=folding_threshold)
# Get real image position of object (deal with chip rotation w.r.t its center)
self.real_pos = self.getRealPos(chip.img, global_x=self.posImg.x, global_y=self.posImg.y,
img_real_wcs=self.chip_wcs)
x, y = self.real_pos.x + 0.5, self.real_pos.y + 0.5
x_nominal = int(np.floor(x + 0.5))
y_nominal = int(np.floor(y + 0.5))
dx = x - x_nominal
dy = y - y_nominal
offset = galsim.PositionD(dx, dy)
# Get real local wcs of object (deal with chip rotation w.r.t its center)
chip_wcs_local = self.chip_wcs.local(self.real_pos)
is_updated = 0
# Loop over all sub-bandpasses
for i in range(len(bandpass_list)):
bandpass = bandpass_list[i]
try:
sub = integrate_sed_bandpass(sed=self.sed, bandpass=bandpass)
except Exception as e:
print(e)
if self.logger:
self.logger.error(e)
continue
ratio = sub / full
if not (ratio == -1 or (ratio != ratio)):
nphotons = ratio * nphotons_tot
else:
continue
# nphotons_sum += nphotons
# print("nphotons_sub-band_%d = %.2f"%(i, nphotons))
# Get PSF model
psf, pos_shear = psf_model.get_PSF(chip=chip, pos_img=pos_img, bandpass=bandpass,
folding_threshold=folding_threshold)
# star = galsim.DeltaFunction(gsparams=gsp)
# star = star.withFlux(nphotons)
# star = galsim.Convolve(psf, star)
star = psf.withFlux(nphotons)
stamp = star.drawImage(wcs=chip_wcs_local, offset=offset)
if np.sum(np.isnan(stamp.array)) > 0:
continue
stamp.setCenter(x_nominal, y_nominal)
bounds = stamp.bounds & galsim.BoundsI(0, chip.npix_x - 1, 0, chip.npix_y - 1)
if bounds.area() > 0:
chip.img.setOrigin(0, 0)
chip.img[bounds] += stamp[bounds]
is_updated = 1
chip.img.setOrigin(chip.bound.xmin, chip.bound.ymin)
del stamp
if is_updated == 0:
# Return code 0: object has missed this detector
print("obj %s missed"%(self.id))
if self.logger:
self.logger.info("obj %s missed"%(self.id))
return 0, pos_shear
return 1, pos_shear # Return code 1: draw sucesss
def addSLStoChipImage(self, sdp=None, chip=None, xOrderSigPlus=None, local_wcs=None):
spec_orders = sdp.compute_spec_orders()
for k, v in spec_orders.items():
img_s = v[0]
#########################################################
# DEBUG
#########################################################
# print("before convolveGaussXorders, img_s:", img_s)
nan_ids = np.isnan(img_s)
if img_s[nan_ids].shape[0] > 0:
# img_s[nan_ids] = 0
print("DEBUG: before convolveGaussXorders specImg nan num is", img_s[nan_ids].shape[0])
#########################################################
img_s, orig_off = convolveGaussXorders(img_s, xOrderSigPlus[k])
origin_order_x = v[1] - orig_off
origin_order_y = v[2] - orig_off
#########################################################
# DEBUG
#########################################################
# print("DEBUG: orig_off is", orig_off)
nan_ids = np.isnan(img_s)
if img_s[nan_ids].shape[0] > 0:
img_s[nan_ids] = 0
print("DEBUG: specImg nan num is", img_s[nan_ids].shape[0])
#########################################################
specImg = galsim.ImageF(img_s)
photons = galsim.PhotonArray.makeFromImage(specImg)
photons.x += origin_order_x
photons.y += origin_order_y
xlen_imf = int(specImg.xmax - specImg.xmin + 1)
ylen_imf = int(specImg.ymax - specImg.ymin + 1)
stamp = galsim.ImageF(xlen_imf, ylen_imf)
stamp.wcs = local_wcs
stamp.setOrigin(origin_order_x, origin_order_y)
bounds = stamp.bounds & galsim.BoundsI(0, chip.npix_x - 1, 0, chip.npix_y - 1)
if bounds.area() == 0:
continue
chip.img.setOrigin(0, 0)
stamp[bounds] = chip.img[bounds]
chip.sensor.accumulate(photons, stamp)
chip.img[bounds] = stamp[bounds]
chip.img.setOrigin(chip.bound.xmin, chip.bound.ymin)
del stamp
del spec_orders
def addSLStoChipImageWithPSF(self, sdp=None, chip=None, pos_img_local = [1,1], psf_model=None, bandNo = 1, grating_split_pos=3685, local_wcs=None, pos_img=None):
spec_orders = sdp.compute_spec_orders()
for k, v in spec_orders.items():
img_s = v[0]
# print(bandNo,k)
try:
psf, pos_shear = psf_model.get_PSF(chip, pos_img_local = pos_img_local, bandNo = bandNo, galsimGSObject=True, g_order = k, grating_split_pos=grating_split_pos)
except:
psf, pos_shear = psf_model.get_PSF(chip=chip, pos_img=pos_img)
psf_img = psf.drawImage(nx=100, ny=100, wcs = local_wcs)
psf_img_m = psf_img.array
#########################################################
# DEBUG
#########################################################
# ids_p = psf_img_m < 0
# psf_img_m[ids_p] = 0
# from astropy.io import fits
# fits.writeto(str(bandNo) + '_' + str(k) + '_psf.fits', psf_img_m)
# print("DEBUG: orig_off is", orig_off)
nan_ids = np.isnan(img_s)
if img_s[nan_ids].shape[0] > 0:
img_s[nan_ids] = 0
print("DEBUG: specImg nan num is", img_s[nan_ids].shape[0])
#########################################################
img_s, orig_off = convolveImg(img_s, psf_img_m)
origin_order_x = v[1] - orig_off[0]
origin_order_y = v[2] - orig_off[1]
specImg = galsim.ImageF(img_s)
# photons = galsim.PhotonArray.makeFromImage(specImg)
# photons.x += origin_order_x
# photons.y += origin_order_y
# xlen_imf = int(specImg.xmax - specImg.xmin + 1)
# ylen_imf = int(specImg.ymax - specImg.ymin + 1)
# stamp = galsim.ImageF(xlen_imf, ylen_imf)
# stamp.wcs = local_wcs
# stamp.setOrigin(origin_order_x, origin_order_y)
specImg.wcs = local_wcs
specImg.setOrigin(origin_order_x, origin_order_y)
bounds = specImg.bounds & galsim.BoundsI(0, chip.npix_x - 1, 0, chip.npix_y - 1)
if bounds.area() == 0:
continue
chip.img.setOrigin(0, 0)
chip.img[bounds] = chip.img[bounds] + specImg[bounds]
# stamp[bounds] = chip.img[bounds]
# # chip.sensor.accumulate(photons, stamp)
# chip.img[bounds] = stamp[bounds]
chip.img.setOrigin(chip.bound.xmin, chip.bound.ymin)
# del stamp
del spec_orders
return pos_shear
def drawObj_slitless(self, tel, pos_img, psf_model, bandpass_list, filt, chip, nphotons_tot=None, g1=0, g2=0,
exptime=150., normFilter=None, grating_split_pos=3685, fd_shear=None):
if normFilter is not None:
norm_thr_rang_ids = normFilter['SENSITIVITY'] > 0.001
sedNormFactor = getNormFactorForSpecWithABMAG(ABMag=self.param['mag_use_normal'], spectrum=self.sed,
norm_thr=normFilter,
sWave=np.floor(normFilter[norm_thr_rang_ids][0][0]),
eWave=np.ceil(normFilter[norm_thr_rang_ids][-1][0]))
if sedNormFactor == 0:
return 2, None
else:
sedNormFactor = 1.
if self.getMagFilter(filt) <= 15:
folding_threshold = 5.e-4
else:
folding_threshold = 5.e-3
gsp = galsim.GSParams(folding_threshold=folding_threshold)
normalSED = Table(np.array([self.sed['WAVELENGTH'], self.sed['FLUX'] * sedNormFactor]).T,
names=('WAVELENGTH', 'FLUX'))
self.real_pos = self.getRealPos(chip.img, global_x=self.posImg.x, global_y=self.posImg.y,
img_real_wcs=self.chip_wcs)
x, y = self.real_pos.x + 0.5, self.real_pos.y + 0.5
x_nominal = int(np.floor(x + 0.5))
y_nominal = int(np.floor(y + 0.5))
dx = x - x_nominal
dy = y - y_nominal
offset = galsim.PositionD(dx, dy)
chip_wcs_local = self.chip_wcs.local(self.real_pos)
flat_cube = chip.flat_cube
xOrderSigPlus = {'A': 1.3909419820029296, 'B': 1.4760376591236062, 'C': 4.035447379743442,
'D': 5.5684364343742825, 'E': 16.260021029735388}
grating_split_pos_chip = 0 + grating_split_pos
branges = np.zeros([len(bandpass_list),2])
if hasattr(psf_model,'bandranges'):
if psf_model.bandranges is None:
return 2, None
if len(psf_model.bandranges) != len(bandpass_list):
return 2, None
branges = psf_model.bandranges
else:
for i in range(len(bandpass_list)):
branges[i, 0] = bandpass_list[i].blue_limit * 10
branges[i, 1] = bandpass_list[i].red_limit * 10
for i in range(len(bandpass_list)):
# bandpass = bandpass_list[i]
brange = branges[i]
# psf, pos_shear = psf_model.get_PSF(chip=chip, pos_img=pos_img, bandpass=bandpass,
# folding_threshold=folding_threshold)
star = galsim.DeltaFunction(gsparams=gsp)
star = star.withFlux(tel.pupil_area * exptime)
psf_tmp = galsim.Gaussian(sigma=0.002)
star = galsim.Convolve(psf_tmp, star)
starImg = star.drawImage(nx=60, ny=60, wcs=chip_wcs_local, offset=offset)
origin_star = [y_nominal - (starImg.center.y - starImg.ymin),
x_nominal - (starImg.center.x - starImg.xmin)]
starImg.setOrigin(0,0)
gal_origin = [origin_star[0], origin_star[1]]
gal_end = [origin_star[0] + starImg.array.shape[0] - 1, origin_star[1] + starImg.array.shape[1] - 1]
if gal_origin[1] < grating_split_pos_chip < gal_end[1]:
subSlitPos = int(grating_split_pos_chip - gal_origin[1] + 1)
## part img disperse
subImg_p1 = starImg.array[:, 0:subSlitPos]
star_p1 = galsim.Image(subImg_p1)
origin_p1 = origin_star
star_p1.setOrigin(0,0)
xcenter_p1 = min(x_nominal, grating_split_pos_chip - 1) - 0
ycenter_p1 = y_nominal - 0
sdp_p1 = SpecDisperser(orig_img=star_p1, xcenter=xcenter_p1,
ycenter=ycenter_p1, origin=origin_p1,
tar_spec=normalSED,
band_start=brange[0], band_end=brange[1],
conf=chip.sls_conf[0],
isAlongY=0,
flat_cube=flat_cube)
# self.addSLStoChipImage(sdp=sdp_p1, chip=chip, xOrderSigPlus=xOrderSigPlus, local_wcs=chip_wcs_local)
pos_shear=self.addSLStoChipImageWithPSF(sdp=sdp_p1, chip=chip, pos_img_local = [xcenter_p1,ycenter_p1],
psf_model=psf_model, bandNo = i+1, grating_split_pos=grating_split_pos,
local_wcs=chip_wcs_local, pos_img = pos_img)
subImg_p2 = starImg.array[:, subSlitPos + 1:starImg.array.shape[1]]
star_p2 = galsim.Image(subImg_p2)
star_p2.setOrigin(0, 0)
origin_p2 = [origin_star[0], grating_split_pos_chip]
xcenter_p2 = max(x_nominal, grating_split_pos_chip - 1) - 0
ycenter_p2 = y_nominal - 0
sdp_p2 = SpecDisperser(orig_img=star_p2, xcenter=xcenter_p2,
ycenter=ycenter_p2, origin=origin_p2,
tar_spec=normalSED,
band_start=brange[0], band_end=brange[1],
conf=chip.sls_conf[1],
isAlongY=0,
flat_cube=flat_cube)
# self.addSLStoChipImage(sdp=sdp_p2, chip=chip, xOrderSigPlus=xOrderSigPlus, local_wcs=chip_wcs_local)
pos_shear=self.addSLStoChipImageWithPSF(sdp=sdp_p2, chip=chip, pos_img_local=[xcenter_p2, ycenter_p2],
psf_model=psf_model, bandNo=i + 1, grating_split_pos=grating_split_pos,
local_wcs=chip_wcs_local, pos_img = pos_img)
del sdp_p1
del sdp_p2
elif grating_split_pos_chip <= gal_origin[1]:
sdp = SpecDisperser(orig_img=starImg, xcenter=x_nominal - 0,
ycenter=y_nominal - 0, origin=origin_star,
tar_spec=normalSED,
band_start=brange[0], band_end=brange[1],
conf=chip.sls_conf[1],
isAlongY=0,
flat_cube=flat_cube)
# self.addSLStoChipImage(sdp=sdp, chip=chip, xOrderSigPlus=xOrderSigPlus, local_wcs=chip_wcs_local)
pos_shear=self.addSLStoChipImageWithPSF(sdp=sdp, chip=chip, pos_img_local=[x_nominal, y_nominal],
psf_model=psf_model, bandNo=i + 1, grating_split_pos=grating_split_pos,
local_wcs=chip_wcs_local, pos_img = pos_img)
del sdp
elif grating_split_pos_chip >= gal_end[1]:
sdp = SpecDisperser(orig_img=starImg, xcenter=x_nominal - 0,
ycenter=y_nominal - 0, origin=origin_star,
tar_spec=normalSED,
band_start=brange[0], band_end=brange[1],
conf=chip.sls_conf[0],
isAlongY=0,
flat_cube=flat_cube)
# self.addSLStoChipImage(sdp=sdp, chip=chip, xOrderSigPlus=xOrderSigPlus, local_wcs=chip_wcs_local)
pos_shear=self.addSLStoChipImageWithPSF(sdp=sdp, chip=chip, pos_img_local=[x_nominal, y_nominal],
psf_model=psf_model, bandNo=i + 1, grating_split_pos=grating_split_pos,
local_wcs=chip_wcs_local, pos_img = pos_img)
del sdp
# del psf
return 1, pos_shear
def SNRestimate(self, img_obj, flux, noise_level=0.0, seed=31415):
img_flux = img_obj.added_flux
stamp = img_obj.copy() * 0.0
rng = galsim.BaseDeviate(seed)
gaussianNoise = galsim.GaussianNoise(rng, sigma=noise_level)
stamp.addNoise(gaussianNoise)
sig_obj = np.std(stamp.array)
snr_obj = img_flux / sig_obj
return snr_obj
def drawObj_PSF(self, tel, pos_img, psf_model, bandpass_list, filt, chip, nphotons_tot=None, g1=0, g2=0,
exptime=150., fd_shear=None, chip_output=None):
if nphotons_tot == None:
nphotons_tot = self.getElectronFluxFilt(filt, tel, exptime)
# print("nphotons_tot = ", nphotons_tot)
try:
full = integrate_sed_bandpass(sed=self.sed, bandpass=filt.bandpass_full)
except Exception as e:
print(e)
if self.logger:
self.logger.error(e)
return 2, None
# Set Galsim Parameters
if self.getMagFilter(filt) <= 15:
folding_threshold = 5.e-4
else:
folding_threshold = 5.e-3
gsp = galsim.GSParams(folding_threshold=folding_threshold)
# Get real image position of object (deal with chip rotation w.r.t its center)
self.real_pos = self.getRealPos(chip.img, global_x=self.posImg.x, global_y=self.posImg.y,
img_real_wcs=self.chip_wcs)
x, y = self.real_pos.x + 0.5, self.real_pos.y + 0.5
x_nominal = int(np.floor(x + 0.5))
y_nominal = int(np.floor(y + 0.5))
dx = x - x_nominal
dy = y - y_nominal
offset = galsim.PositionD(dx, dy)
# Get real local wcs of object (deal with chip rotation w.r.t its center)
chip_wcs_local = self.chip_wcs.local(self.real_pos)
is_updated = 0
# Loop over all sub-bandpasses
for i in range(len(bandpass_list)):
bandpass = bandpass_list[i]
try:
sub = integrate_sed_bandpass(sed=self.sed, bandpass=bandpass)
except Exception as e:
print(e)
if self.logger:
self.logger.error(e)
continue
ratio = sub / full
if not (ratio == -1 or (ratio != ratio)):
nphotons = ratio * nphotons_tot
else:
continue
# Get PSF model
psf, pos_shear = psf_model.get_PSF(chip=chip, pos_img=pos_img, bandpass=bandpass,
folding_threshold=folding_threshold)
star_temp = psf.withFlux(nphotons)
if i==0:
star = star_temp
else:
star = star+star_temp
pixelScale = 0.074
stamp = star.drawImage(wcs=chip_wcs_local, offset=offset)
#stamp = star.drawImage(nx=256, ny=256, scale=pixelScale)
if np.sum(np.isnan(stamp.array)) > 0:
return None
fn = chip_output.subdir + "/psfIDW"
os.makedirs(fn, exist_ok=True)
fn = fn + "/ccd_{:}".format(chip.chipID)+"_psf_"+str(self.param['id'])+".fits"
if fn != None:
if os.path.exists(fn):
os.remove(fn)
hdu = fitsio.PrimaryHDU()
hdu.data = stamp.array
hdu.header.set('name', self.type)
hdu.header.set('pixScale', pixelScale)
hdu.header.set('objID', self.param['id'])
hdu.writeto(fn)
del stamp
return None
ObservationSim/MockObject/Quasar.py deleted 100755 → 0
View file @ 4afd1181
import galsim
import os
import sys
import numpy as np
import astropy.constants as cons
from astropy.table import Table
from scipy import interpolate
from ObservationSim.MockObject.MockObject import MockObject
from ObservationSim.MockObject._util import integrate_sed_bandpass, getNormFactorForSpecWithABMAG, getObservedSED, getABMAG
class Quasar(MockObject):
def __init__(self, param, logger=None):
super().__init__(param, logger=logger)
if not hasattr(self, "mu"):
if hasattr(self, "detA"):
self.mu = 1./self.detA
else:
self.mu = 1.
def load_SED(self, survey_type, sed_path=None, cosids=None, objtypes=None, sed_templates=None, normFilter=None, target_filt=None):
'''
--------------------------------------------------------------------
(Deprecated) Left over codes from cycle-1 to cycle-3 implementations
--------------------------------------------------------------------
'''
if survey_type == "photometric":
norm_thr_rang_ids = normFilter['SENSITIVITY'] > 0.001
if sed_templates is None:
# Read SED data directly
itype = objtypes[cosids == self.sed_type][0]
sed_file = os.path.join(
sed_path, itype + "_ID%s.sed" % (self.sed_type))
if not os.path.exists(sed_file):
raise ValueError("!!! No SED found.")
sed_data = Table.read(sed_file, format="ascii")
wave, flux = sed_data["observedLambda"].data, sed_data["observedFlux"].data
else:
# Load SED from templates
sed_data = sed_templates[self.sed_type]
# redshift, intrinsic extinction
sed_data = getObservedSED(
sedCat=sed_data,
redshift=self.z,
av=self.param['av'],
redden=self.param['redden'])
wave, flux = sed_data[0], sed_data[1]
flux_photon = flux * (wave / (cons.h.value * cons.c.value)) * 1e-13
sed_photon = Table(
np.array([wave, flux_photon]).T, names=('WAVELENGTH', 'FLUX'))
# Get scaling factor for SED
sedNormFactor = getNormFactorForSpecWithABMAG(ABMag=self.param['mag_use_normal'],
spectrum=sed_photon,
norm_thr=normFilter,
sWave=np.floor(
normFilter[norm_thr_rang_ids][0][0]),
eWave=np.ceil(normFilter[norm_thr_rang_ids][-1][0]))
sed_photon = np.array(
[sed_photon['WAVELENGTH'], sed_photon['FLUX']*sedNormFactor]).T
# Convert to galsim.SED object
spec = galsim.LookupTable(x=np.array(sed_photon[:, 0]), f=np.array(
sed_photon[:, 1]), interpolant='nearest')
self.sed = galsim.SED(spec, wave_type='A',
flux_type='1', fast=False)
# Get magnitude
interFlux = integrate_sed_bandpass(
sed=self.sed, bandpass=target_filt.bandpass_full)
self.param['mag_%s' % target_filt.filter_type] = getABMAG(
interFlux=interFlux,
bandpass=target_filt.bandpass_full)
# print('mag_use_normal = ', self.param['mag_use_normal'])
# print('mag_%s = '%target_filt.filter_type, self.param['mag_%s'%target_filt.filter_type])
elif survey_type == "spectroscopic":
if sed_templates is None:
self.sedPhotons(sed_path=sed_path,
cosids=cosids, objtypes=objtypes)
else:
sed_data = sed_templates[self.sed_type]
sed_data = getObservedSED(
sedCat=sed_data,
redshift=self.z,
av=self.param['av'],
redden=self.param['redden'])
speci = interpolate.interp1d(sed_data[0], sed_data[1])
lamb = np.arange(2500, 10001 + 0.5, 0.5)
y = speci(lamb)
# erg/s/cm2/A --> photo/s/m2/A
all_sed = y * lamb / (cons.h.value * cons.c.value) * 1e-13
self.sed = Table(
np.array([lamb, all_sed]).T, names=('WAVELENGTH', 'FLUX'))
def unload_SED(self):
"""(Test) free up SED memory
"""
del self.sed
def getGSObj(self, psf, g1=0, g2=0, flux=None, filt=None, tel=None, exptime=150.):
if flux == None:
flux = self.getElectronFluxFilt(filt, tel, exptime)
qso = galsim.Gaussian(sigma=1.e-8, flux=1.)
qso = qso.withFlux(flux)
final = galsim.Convolve(psf, qso)
return final
ObservationSim/MockObject/SpecDisperser/SpecDisperser.py deleted 100644 → 0
View file @ 4afd1181
from astropy.table import Table
import astropy.constants as cons
import collections
from collections import OrderedDict
from scipy import interpolate
from astropy import units as u
# from numpy import *
import numpy as np
from pylab import *
import galsim
import sys
import os
def rotate90(array_orig=None, xc=0, yc=0, isClockwise=0):
if array_orig is None:
return
l1 = array_orig.shape[0]
l2 = array_orig.shape[1]
if xc < 0 or xc > l2 - 1 or yc < 0 or yc > l1 - 1:
return
n_xc = xc
n_yc = yc
array_final = np.zeros_like(array_orig.T)
if isClockwise == 1:
for i in np.arange(l2):
array_final[i] = array_orig[:, l2 - i - 1]
n_xc = yc
n_yc = l2 - 1 - xc
else:
for i in np.arange(l2):
array_final[i] = array_orig[::-1, i]
n_xc = l1 - 1 - yc
n_yc = xc
return array_final, n_xc, n_yc
class SpecDisperser(object):
def __init__(self, orig_img=None, xcenter=0, ycenter=0, origin=[100, 100], tar_spec=None, band_start=6200,
band_end=10000, isAlongY=0, conf='../param/CONF/csst.conf', gid=0, flat_cube=None, ignoreBeam=[]):
"""
orig_img: normal image,galsim image
xcenter, ycenter: the center of galaxy in orig_img
origin : [int, int]
`origin` defines the lower left pixel index (y,x) of the `direct`
cutout from a larger detector-frame image
tar_spec: galsim.SED
"""
# self.img_x = orig_img.shape[1]
# self.img_y = orig_img.shape[0]
self.thumb_img = np.abs(orig_img.array)
self.thumb_x = orig_img.center.x
self.thumb_y = orig_img.center.y
self.img_sh = orig_img.array.shape
self.id = gid
self.xcenter = xcenter
self.ycenter = ycenter
self.isAlongY = isAlongY
self.flat_cube = flat_cube
if self.isAlongY == 1:
self.thumb_img, self.thumb_x, self.thumb_y = rotate90(array_orig=self.thumb_img, xc=orig_img.center.x,
yc=orig_img.center.y, isClockwise=1)
self.img_sh = orig_img.array.T.shape
self.xcenter = ycenter
self.ycenter = xcenter
self.origin = origin
self.band_start = band_start
self.band_end = band_end
self.spec = tar_spec
self.beam_flux = OrderedDict()
self.grating_conf = aXeConf(conf)
self.grating_conf.get_beams()
self.grating_conf_file = conf
self.ignoreBeam = ignoreBeam
def compute_spec_orders(self):
all_orders = OrderedDict()
beam_names = self.grating_conf.beams
for beam in beam_names:
if beam in self.ignoreBeam:
continue
all_orders[beam] = self.compute_spec(beam)
return all_orders
def compute_spec(self, beam):
from .disperse_c import interp
from .disperse_c import disperse
# from MockObject.disperse_c import disperse
dx = self.grating_conf.dxlam[beam]
xoff = 0
ytrace_beam, lam_beam = self.grating_conf.get_beam_trace(x=self.xcenter, y=self.ycenter, dx=(dx + xoff),
beam=beam)
### Account for pixel centering of the trace
yfrac_beam = ytrace_beam - floor(ytrace_beam+0.5)
ysens = lam_beam * 0
lam_index = argsort(lam_beam)
conf_sens = self.grating_conf.sens[beam]
lam_intep = np.linspace(self.band_start, self.band_end, int((self.band_end - self.band_start) / 0.1))
thri = interpolate.interp1d(conf_sens['WAVELENGTH'], conf_sens['SENSITIVITY'])
spci = interpolate.interp1d(self.spec['WAVELENGTH'], self.spec['FLUX'])
beam_thr = thri(lam_intep)
spec_sample = spci(lam_intep)
bean_thr_spec = beam_thr * spec_sample
###generate sensitivity file for aXe
# ysensitivity = lam_beam * 0
#
# ysensitivity[lam_index] = interp.interp_conserve_c(lam_beam[lam_index], lam_intep,
# beam_thr * math.pi * 100 * 100 * 1e-7 / (
# cons.h.value * cons.c.value / (
# lam_intep * 1e-10)), integrate=0, left=0,
# right=0)
#
# self.writerSensitivityFile(conffile = self.grating_conf_file, beam = beam, w = lam_beam[lam_index], sens = ysensitivity[lam_index])
ysens[lam_index] = interp.interp_conserve_c(lam_beam[lam_index], lam_intep, bean_thr_spec, integrate=1, left=0,
right=0)
sensitivity_beam = ysens
len_spec_x = len(dx)
len_spec_y = int(abs(ceil(ytrace_beam[-1]) - floor(ytrace_beam[0])) + 1)
beam_sh = (self.img_sh[0] + len_spec_y, self.img_sh[1] + len_spec_x)
modelf = zeros(product(beam_sh), dtype=float)
model = modelf.reshape(beam_sh)
idx = np.arange(modelf.size, dtype=int64).reshape(beam_sh)
x0 = array((self.thumb_y, self.thumb_x), dtype=int64)
dxpix = dx - dx[0] + x0[1]
dyc = cast[int](np.floor(ytrace_beam+0.5))
dypix = cast[int](np.floor(ytrace_beam - dyc[0] + x0[0] + 0.5))
frac_ids = yfrac_beam<0
dypix[frac_ids] = dypix[frac_ids] - 1
yfrac_beam[frac_ids] = 1+yfrac_beam[frac_ids]
flat_index = idx[dypix, dxpix]
nonz = sensitivity_beam != 0
origin_in = zeros_like(self.origin)
dx0_in = dx[0]
dy0_in = dyc[0]
if self.isAlongY == 1:
origin_in[0] = self.origin[0]
origin_in[1] = self.origin[1] - len_spec_y
dx0_in = -dyc[0]
dy0_in = dx[0]
else:
origin_in[0] = self.origin[0]
origin_in[1] = self.origin[1]
dx0_in = dx[0]
dy0_in = dyc[0]
originOut_x = origin_in[1] + dx0_in
originOut_y = origin_in[0] + dy0_in
if self.flat_cube is None:
beam_flat = None
else:
beam_flat = zeros([len(modelf), len(self.flat_cube)])
sub_flat_cube = zeros([len(self.flat_cube),beam_sh[0], beam_sh[1]])
sub_flat_cube[0] = sub_flat_cube[0] + 1.
overlap_flag = 1
sub_y_s = originOut_y
sub_y_e = originOut_y + beam_sh[0] - 1
sub_x_s = originOut_x
sub_x_e = originOut_x + beam_sh[1] - 1
beam_x_s = max(sub_x_s, 0)
if beam_x_s > self.flat_cube[0].shape[1] - 1: overlap_flag = 0
if overlap_flag == 1:
beam_x_e = min(sub_x_e, self.flat_cube[0].shape[1] - 1)
if beam_x_e < 0: overlap_flag = 0
if overlap_flag == 1:
beam_y_s = max(sub_y_s, 0)
if beam_y_s > self.flat_cube[0].shape[0] - 1: overlap_flag = 0
if overlap_flag == 1:
beam_y_e = min(sub_y_e, self.flat_cube[0].shape[0] - 1)
if beam_y_e < 0: overlap_flag = 0
if overlap_flag == 1:
sub_flat_cube[:,beam_y_s-originOut_y:beam_y_e-originOut_y+1,beam_x_s-originOut_x:beam_x_e-originOut_x+1] = self.flat_cube[:,beam_y_s:beam_y_e+1,beam_x_s:beam_x_e+1]
for i in arange(0, len(self.flat_cube), 1):
beam_flat[:,i] = sub_flat_cube[i].flatten()
# beam_flat = zeros([len(modelf), len(self.flat_cube)])
# flat_sh = self.flat_cube[0].shape
# for i in arange(0, beam_sh[0], 1):
# for j in arange(0, beam_sh[1], 1):
# k = i * beam_sh[1] + j
# if originOut_y + i >= flat_sh[0] or originOut_y + i < 0 or originOut_x + j>= flat_sh[1] or originOut_x+j < 0:
# temp_bf = np.zeros_like(self.flat_cube[:, 0, 0])
# temp_bf[0] = 1.0
# beam_flat[k] = temp_bf
# else:
# beam_flat[k] = self.flat_cube[:, originOut_y + i, originOut_x + j]
status = disperse.disperse_grism_object(self.thumb_img.astype(np.float32),
flat_index[nonz], yfrac_beam[nonz],
sensitivity_beam[nonz],
modelf, x0,
array(self.img_sh, dtype=int64),
array(beam_sh, dtype=int64),
beam_flat,
lam_beam[lam_index][nonz])
model = modelf.reshape(beam_sh)
self.beam_flux[beam] = sum(modelf)
if self.isAlongY == 1:
model, _, _ = rotate90(array_orig=model, isClockwise=0)
return model, originOut_x, originOut_y, dxpix, dypix, lam_beam,ysens
def writerSensitivityFile(self, conffile='', beam='', w=None, sens=None):
orders = {'A': '1st', 'B': '0st', 'C': '2st', 'D': '-1st', 'E': '-2st'}
sens_file_name = conffile[0:-5] + '_sensitivity_' + orders[beam] + '.fits'
if not os.path.exists(sens_file_name) == True:
senstivity_out = Table(array([w, sens]).T, names=('WAVELENGTH', 'SENSITIVITY'))
senstivity_out.write(sens_file_name, format='fits')
"""
Demonstrate aXe trace polynomials.
"""
class aXeConf():
def __init__(self, conf_file='WFC3.IR.G141.V2.5.conf'):
"""Read an aXe-compatible configuration file
Parameters
----------
conf_file: str
Filename of the configuration file to read
"""
if conf_file is not None:
self.conf = self.read_conf_file(conf_file)
self.conf_file = conf_file
self.count_beam_orders()
## Global XOFF/YOFF offsets
if 'XOFF' in self.conf.keys():
self.xoff = np.float(self.conf['XOFF'])
else:
self.xoff = 0.
if 'YOFF' in self.conf.keys():
self.yoff = np.float(self.conf['YOFF'])
else:
self.yoff = 0.
def read_conf_file(self, conf_file='WFC3.IR.G141.V2.5.conf'):
"""Read an aXe config file, convert floats and arrays
Parameters
----------
conf_file: str
Filename of the configuration file to read.
Parameters are stored in an OrderedDict in `self.conf`.
"""
from collections import OrderedDict
conf = OrderedDict()
lines = open(conf_file).readlines()
for line in lines:
## empty / commented lines
if (line.startswith('#')) | (line.strip() == '') | ('"' in line):
continue
## split the line, taking out ; and # comments
spl = line.split(';')[0].split('#')[0].split()
param = spl[0]
if len(spl) > 2:
value = np.cast[float](spl[1:])
else:
try:
value = float(spl[1])
except:
value = spl[1]
conf[param] = value
return conf
def count_beam_orders(self):
"""Get the maximum polynomial order in DYDX or DLDP for each beam
"""
self.orders = {}
for beam in ['A', 'B', 'C', 'D', 'E', 'F', 'G', 'H', 'I', 'J']:
order = 0
while 'DYDX_{0:s}_{1:d}'.format(beam, order) in self.conf.keys():
order += 1
while 'DLDP_{0:s}_{1:d}'.format(beam, order) in self.conf.keys():
order += 1
self.orders[beam] = order - 1
def get_beams(self):
"""Get beam parameters and read sensitivity curves
"""
import os
from collections import OrderedDict
from astropy.table import Table, Column
self.dxlam = OrderedDict()
self.nx = OrderedDict()
self.sens = OrderedDict()
self.beams = []
for beam in self.orders:
if self.orders[beam] > 0:
self.beams.append(beam)
self.dxlam[beam] = np.arange(self.conf['BEAM{0}'.format(beam)].min(),
self.conf['BEAM{0}'.format(beam)].max(), dtype=int)
self.nx[beam] = int(self.dxlam[beam].max() - self.dxlam[beam].min()) + 1
self.sens[beam] = Table.read(
'{0}/{1}'.format(os.path.dirname(self.conf_file), self.conf['SENSITIVITY_{0}'.format(beam)]))
# self.sens[beam].wave = np.cast[np.double](self.sens[beam]['WAVELENGTH'])
# self.sens[beam].sens = np.cast[np.double](self.sens[beam]['SENSITIVITY'])
### Need doubles for interpolating functions
for col in self.sens[beam].colnames:
data = np.cast[np.double](self.sens[beam][col])
self.sens[beam].remove_column(col)
self.sens[beam].add_column(Column(data=data, name=col))
self.beams.sort()
def field_dependent(self, xi, yi, coeffs):
"""aXe field-dependent coefficients
See the `aXe manual <http://axe.stsci.edu/axe/manual/html/node7.html#SECTION00721200000000000000>`_ for a description of how the field-dependent coefficients are specified.
Parameters
----------
xi, yi : float or array-like
Coordinate to evaluate the field dependent coefficients, where
`xi = x-REFX` and `yi = y-REFY`.
coeffs : array-like
Field-dependency coefficients
Returns
-------
a : float or array-like
Evaluated field-dependent coefficients
"""
## number of coefficients for a given polynomial order
## 1:1, 2:3, 3:6, 4:10, order:order*(order+1)/2
if isinstance(coeffs, float):
order = 1
else:
order = int(-1 + np.sqrt(1 + 8 * len(coeffs))) // 2
## Build polynomial terms array
## $a = a_0+a_1x_i+a_2y_i+a_3x_i^2+a_4x_iy_i+a_5yi^2+$ ...
xy = []
for p in range(order):
for px in range(p + 1):
# print 'x**%d y**%d' %(p-px, px)
xy.append(xi ** (p - px) * yi ** (px))
## Evaluate the polynomial, allowing for N-dimensional inputs
a = np.sum((np.array(xy).T * coeffs).T, axis=0)
return a
def evaluate_dp(self, dx, dydx):
"""Evalate arc length along the trace given trace polynomial coefficients
Parameters
----------
dx : array-like
x pixel to evaluate
dydx : array-like
Coefficients of the trace polynomial
Returns
-------
dp : array-like
Arc length along the trace at position `dx`.
For `dydx` polynomial orders 0, 1 or 2, integrate analytically.
Higher orders must be integrated numerically.
**Constant:**
.. math:: dp = dx
**Linear:**
.. math:: dp = \sqrt{1+\mathrm{DYDX}[1]}\cdot dx
**Quadratic:**
.. math:: u = \mathrm{DYDX}[1] + 2\ \mathrm{DYDX}[2]\cdot dx
.. math:: dp = (u \sqrt{1+u^2} + \mathrm{arcsinh}\ u) / (4\cdot \mathrm{DYDX}[2])
"""
## dp is the arc length along the trace
## $\lambda = dldp_0 + dldp_1 dp + dldp_2 dp^2$ ...
poly_order = len(dydx) - 1
if (poly_order == 2):
if np.abs(np.unique(dydx[2])).max() == 0:
poly_order = 1
if poly_order == 0: ## dy=0
dp = dx
elif poly_order == 1: ## constant dy/dx
dp = np.sqrt(1 + dydx[1] ** 2) * (dx)
elif poly_order == 2: ## quadratic trace
u0 = dydx[1] + 2 * dydx[2] * (0)
dp0 = (u0 * np.sqrt(1 + u0 ** 2) + np.arcsinh(u0)) / (4 * dydx[2])
u = dydx[1] + 2 * dydx[2] * (dx)
dp = (u * np.sqrt(1 + u ** 2) + np.arcsinh(u)) / (4 * dydx[2]) - dp0
else:
## high order shape, numerical integration along trace
## (this can be slow)
xmin = np.minimum((dx).min(), 0)
xmax = np.maximum((dx).max(), 0)
xfull = np.arange(xmin, xmax)
dyfull = 0
for i in range(1, poly_order):
dyfull += i * dydx[i] * (xfull - 0.5) ** (i - 1)
## Integrate from 0 to dx / -dx
dpfull = xfull * 0.
lt0 = xfull < 0
if lt0.sum() > 1:
dpfull[lt0] = np.cumsum(np.sqrt(1 + dyfull[lt0][::-1] ** 2))[::-1]
dpfull[lt0] *= -1
#
gt0 = xfull > 0
if gt0.sum() > 0:
dpfull[gt0] = np.cumsum(np.sqrt(1 + dyfull[gt0] ** 2))
dp = np.interp(dx, xfull, dpfull)
if dp[-1] == dp[-2]:
dp[-1] = dp[-2] + np.diff(dp)[-2]
return dp
def get_beam_trace(self, x=507, y=507, dx=0., beam='A'):
"""Get an aXe beam trace for an input reference pixel and list of output x pixels `dx`
Parameters
----------
x, y : float or array-like
Evaluate trace definition at detector coordinates `x` and `y`.
dx : float or array-like
Offset in x pixels from `(x,y)` where to compute trace offset and
effective wavelength
beam : str
Beam name (i.e., spectral order) to compute. By aXe convention,
`beam='A'` is the first order, 'B' is the zeroth order and
additional beams are the higher positive and negative orders.
Returns
-------
dy : float or array-like
Center of the trace in y pixels offset from `(x,y)` evaluated at
`dx`.
lam : float or array-like
Effective wavelength along the trace evaluated at `dx`.
"""
NORDER = self.orders[beam] + 1
xi, yi = x - self.xoff, y - self.yoff
xoff_beam = self.field_dependent(xi, yi, self.conf['XOFF_{0}'.format(beam)])
yoff_beam = self.field_dependent(xi, yi, self.conf['YOFF_{0}'.format(beam)])
## y offset of trace (DYDX)
dydx = np.zeros(NORDER) # 0 #+1.e-80
dydx = [0] * NORDER
for i in range(NORDER):
if 'DYDX_{0:s}_{1:d}'.format(beam, i) in self.conf.keys():
coeffs = self.conf['DYDX_{0:s}_{1:d}'.format(beam, i)]
dydx[i] = self.field_dependent(xi, yi, coeffs)
# $dy = dydx_0+dydx_1 dx+dydx_2 dx^2+$ ...
dy = yoff_beam
for i in range(NORDER):
dy += dydx[i] * (dx - xoff_beam) ** i
## wavelength solution
dldp = np.zeros(NORDER)
dldp = [0] * NORDER
for i in range(NORDER):
if 'DLDP_{0:s}_{1:d}'.format(beam, i) in self.conf.keys():
coeffs = self.conf['DLDP_{0:s}_{1:d}'.format(beam, i)]
dldp[i] = self.field_dependent(xi, yi, coeffs)
self.eval_input = {'x': x, 'y': y, 'beam': beam, 'dx': dx}
self.eval_output = {'xi': xi, 'yi': yi, 'dldp': dldp, 'dydx': dydx,
'xoff_beam': xoff_beam, 'yoff_beam': yoff_beam,
'dy': dy}
dp = self.evaluate_dp(dx - xoff_beam, dydx)
# ## dp is the arc length along the trace
# ## $\lambda = dldp_0 + dldp_1 dp + dldp_2 dp^2$ ...
# if self.conf['DYDX_ORDER_%s' %(beam)] == 0: ## dy=0
# dp = dx-xoff_beam
# elif self.conf['DYDX_ORDER_%s' %(beam)] == 1: ## constant dy/dx
# dp = np.sqrt(1+dydx[1]**2)*(dx-xoff_beam)
# elif self.conf['DYDX_ORDER_%s' %(beam)] == 2: ## quadratic trace
# u0 = dydx[1]+2*dydx[2]*(0)
# dp0 = (u0*np.sqrt(1+u0**2)+np.arcsinh(u0))/(4*dydx[2])
# u = dydx[1]+2*dydx[2]*(dx-xoff_beam)
# dp = (u*np.sqrt(1+u**2)+np.arcsinh(u))/(4*dydx[2])-dp0
# else:
# ## high order shape, numerical integration along trace
# ## (this can be slow)
# xmin = np.minimum((dx-xoff_beam).min(), 0)
# xmax = np.maximum((dx-xoff_beam).max(), 0)
# xfull = np.arange(xmin, xmax)
# dyfull = 0
# for i in range(1, NORDER):
# dyfull += i*dydx[i]*(xfull-0.5)**(i-1)
#
# ## Integrate from 0 to dx / -dx
# dpfull = xfull*0.
# lt0 = xfull <= 0
# if lt0.sum() > 1:
# dpfull[lt0] = np.cumsum(np.sqrt(1+dyfull[lt0][::-1]**2))[::-1]
# dpfull[lt0] *= -1
# #
# gt0 = xfull >= 0
# if gt0.sum() > 0:
# dpfull[gt0] = np.cumsum(np.sqrt(1+dyfull[gt0]**2))
#
# dp = np.interp(dx-xoff_beam, xfull, dpfull)
## Evaluate dldp
lam = dp * 0.
for i in range(NORDER):
lam += dldp[i] * dp ** i
return dy, lam
def show_beams(self, beams=['E', 'D', 'C', 'B', 'A']):
"""
Make a demo plot of the beams of a given configuration file
"""
import matplotlib.pyplot as plt
x0, x1 = 507, 507
dx = np.arange(-800, 1200)
if 'WFC3.UV' in self.conf_file:
x0, x1 = 2073, 250
dx = np.arange(-1200, 1200)
if 'G800L' in self.conf_file:
x0, x1 = 2124, 1024
dx = np.arange(-1200, 1200)
s = 200 # marker size
fig = plt.figure(figsize=[10, 3])
plt.scatter(0, 0, marker='s', s=s, color='black', edgecolor='0.8',
label='Direct')
for beam in beams:
if 'XOFF_{0}'.format(beam) not in self.conf.keys():
continue
xoff = self.field_dependent(x0, x1, self.conf['XOFF_{0}'.format(beam)])
dy, lam = self.get_beam_trace(x0, x1, dx=dx, beam=beam)
xlim = self.conf['BEAM{0}'.format(beam)]
ok = (dx >= xlim[0]) & (dx <= xlim[1])
plt.scatter(dx[ok] + xoff, dy[ok], c=lam[ok] / 1.e4, marker='s', s=s,
alpha=0.5, edgecolor='None')
plt.text(np.median(dx[ok]), np.median(dy[ok]) + 1, beam,
ha='center', va='center', fontsize=14)
print('Beam {0}, lambda=({1:.1f} - {2:.1f})'.format(beam, lam[ok].min(), lam[ok].max()))
plt.grid()
plt.xlabel(r'$\Delta x$')
plt.ylabel(r'$\Delta y$')
cb = plt.colorbar(pad=0.01, fraction=0.05)
cb.set_label(r'$\lambda\,(\mu\mathrm{m})$')
plt.title(self.conf_file)
plt.tight_layout()
plt.savefig('{0}.pdf'.format(self.conf_file))
# def load_grism_config(conf_file):
# """Load parameters from an aXe configuration file
# Parameters
# ----------
# conf_file : str
# Filename of the configuration file
# Returns
# -------
# conf : `~grizli.grismconf.aXeConf`
# Configuration file object. Runs `conf.get_beams()` to read the
# sensitivity curves.
# """
# conf = aXeConf(conf_file)
# conf.get_beams()
# return conf
ObservationSim/MockObject/SpecDisperser/__init__.py deleted 100644 → 0
View file @ 4afd1181
from .SpecDisperser import *
from .disperse_c import disperse, interp
\ No newline at end of file
ObservationSim/MockObject/SpecDisperser/disperse_c/__init__.py deleted 100644 → 0
View file @ 4afd1181
#!/usr/bin/env python
# encoding: utf-8
"""
Cython
"""
# from . import disperse
# from . import interp
from .disperse import *
from .interp import *
ObservationSim/MockObject/SpecDisperser/disperse_c/disperse.pyx deleted 100644 → 0
View file @ 4afd1181
from __future__ import division
import numpy as np
cimport numpy as np
DTYPE = np.double
ITYPE = np.int64
ctypedef np.double_t DTYPE_t
ctypedef np.uint_t UINT_t
ctypedef np.int_t INT_t
ctypedef np.int64_t LINT_t
ctypedef np.int32_t FINT_t
ctypedef np.float32_t FTYPE_t
import cython
cdef extern from "math.h":
double sqrt(double x)
double exp(double x)
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.embedsignature(True)
def disperse_grism_object(np.ndarray[FTYPE_t, ndim=2] flam,
np.ndarray[LINT_t, ndim=1] idxl,
np.ndarray[DTYPE_t, ndim=1] yfrac,
np.ndarray[DTYPE_t, ndim=1] ysens,
np.ndarray[DTYPE_t, ndim=1] full,
np.ndarray[LINT_t, ndim=1] x0,
np.ndarray[LINT_t, ndim=1] shd,
np.ndarray[LINT_t, ndim=1] shg,
np.ndarray[DTYPE_t, ndim=2] flat,
np.ndarray[DTYPE_t, ndim=1] wlambda):
"""Compute a dispersed 2D spectrum
Parameters
----------
flam : direct image matrix, 2 dim [y,x]
idxl: grating disperse light to pixel, pixel index, 1 dim, length = ysens, yfrac
yfrac:
ysens: sensitivity use pixel describe
full： output result ，1 dim, y_beam * x_beam
x0: the center of gal in image thumbnail (y, x)
shd: shape of direct image
shg: shape of grating image
"""
cdef int i,j,k1,k2
cdef unsigned int nk,nl,k,shx,shy
cdef double fl_ij
nk = len(idxl)
nl = len(full)
if (flat is not None):
nlamb = len(wlambda)
nflat = len(flat)
flat_lamb_min = 2500
flat_lamb_max = 10000
lambda_co = np.zeros([len(flat[0]),nlamb])
lambda_co[0] = lambda_co[0] + ysens
lambda_co[1] = (wlambda-flat_lamb_min)/(flat_lamb_max-flat_lamb_min)*ysens
lambda_co[2] = lambda_co[1]*lambda_co[1]*ysens
lambda_co[3] = lambda_co[1]*lambda_co[2]*ysens
flat_eff_all = np.zeros(nflat*nlamb)
for i in range(0, nflat):
flat_eff_all[i*nlamb:(i+1)*nlamb]=np.dot(flat[i], lambda_co)
for i in range(0-x0[1], x0[1]):
if (x0[1]+i < 0) | (x0[1]+i >= shd[1]):
continue
for j in range(0-x0[0], x0[0]):
if (x0[0]+j < 0) | (x0[0]+j >= shd[0]):
continue
fl_ij = flam[x0[0]+j, x0[1]+i] #/1.e-17
if (fl_ij == 0):
continue
for k in range(nk):
k1 = idxl[k]+j*shg[1]+i
if (k1 >= 0) & (k1 < nl):
flat_ids = k1*nlamb+k
full[k1] += fl_ij*(1-yfrac[k])*flat_eff_all[flat_ids]
k2 = idxl[k]+(j+1)*shg[1]+i
if (k2 >= 0) & (k2 < nl):
flat_ids = k2*nlamb+k
full[k2] += fl_ij*yfrac[k]*flat_eff_all[flat_ids]
else:
for i in range(0-x0[1], x0[1]):
if (x0[1]+i < 0) | (x0[1]+i >= shd[1]):
continue
for j in range(0-x0[0], x0[0]):
if (x0[0]+j < 0) | (x0[0]+j >= shd[0]):
continue
fl_ij = flam[x0[0]+j, x0[1]+i] #/1.e-17
if (fl_ij == 0):
continue
for k in range(nk):
k1 = idxl[k]+j*shg[1]+i
if (k1 >= 0) & (k1 < nl):
full[k1] += ysens[k]*fl_ij*(1-yfrac[k])
k2 = idxl[k]+(j+1)*shg[1]+i
if (k2 >= 0) & (k2 < nl):
full[k2] += ysens[k]*fl_ij*yfrac[k]
return True
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.embedsignature(True)
def compute_segmentation_limits(np.ndarray[FTYPE_t, ndim=2] segm, int seg_id, np.ndarray[FTYPE_t, ndim=2] flam, np.ndarray[LINT_t, ndim=1] shd):
"""Find pixel limits of a segmentation region
Parameters
----------
segm: ndarray (np.float32)
segmentation array
seg_id: int
ID to test
flam: ndarray (float)
Flux array to compute weighted centroid within segmentation region
shd: [int, int]
Shape of segm
"""
cdef int i, j, imin, imax, jmin, jmax, area
cdef double inumer, jnumer, denom, wht_ij
area = 0
imin = shd[0]
imax = 0
jmin = shd[1]
jmax = 0
inumer = 0.
jnumer = 0.
denom = 0.
for i in range(shd[0]):
for j in range(shd[1]):
if segm[i,j] != seg_id:
continue
area += 1
wht_ij = flam[i,j]
inumer += i*wht_ij
jnumer += j*wht_ij
denom += wht_ij
if i < imin:
imin = i
if i > imax:
imax = i
if j < jmin:
jmin = j
if j > jmax:
jmax = j
### No matched pixels
if denom == 0:
denom = -99
return imin, imax, inumer/denom, jmin, jmax, jnumer/denom, area, denom
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.embedsignature(True)
def seg_flux(np.ndarray[FTYPE_t, ndim=2] flam, np.ndarray[LINT_t, ndim=1] idxl, np.ndarray[DTYPE_t, ndim=1] yfrac, np.ndarray[DTYPE_t, ndim=1] ysens, np.ndarray[DTYPE_t, ndim=1] full, np.ndarray[LINT_t, ndim=1] x0, np.ndarray[LINT_t, ndim=1] shd, np.ndarray[LINT_t, ndim=1] shg):
pass
ObservationSim/MockObject/SpecDisperser/disperse_c/interp.pyx deleted 100644 → 0
View file @ 4afd1181
"""
Pythonic utilities ported to C [Cython] for speedup.
"""
import numpy as np
cimport numpy as np
DTYPE = np.double
ctypedef np.double_t DTYPE_t
ctypedef np.int_t ITYPE_t
ctypedef np.uint_t UINT_t
import cython
cdef extern from "math.h":
double fabs(double)
cdef extern from"stdio.h":
extern int printf(const char *format, ...)
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
@cython.embedsignature(True)
def interp_c(np.ndarray[DTYPE_t, ndim=1] x, np.ndarray[DTYPE_t, ndim=1] xp, np.ndarray[DTYPE_t, ndim=1] fp, double extrapolate=0., short int assume_sorted=1):
"""
interp_c(x, xp, fp, extrapolate=0., assume_sorted=0)
Fast interpolation: [`xp`, `fp`] interpolated at `x`.
Extrapolated values are set to `extrapolate`.
The default `assume_sorted`=1 assumes that the `x` array is sorted and single-
valued, providing a significant gain in speed. (xp is always assumed to be sorted)
"""
cdef unsigned long i, j, N, Np
cdef DTYPE_t x1,x2,y1,y2,out
cdef DTYPE_t fout, xval, xmin
N, Np = len(x), len(xp)
cdef np.ndarray[DTYPE_t, ndim=1] f = np.zeros(N)
i=0
j=0
### Handle left extrapolation
xmin = xp[0]
if assume_sorted == 1:
while x[j] < xmin:
f[j] = extrapolate
j+=1
if j>=N:
break
while j < N:
xval = x[j]
if assume_sorted == 0:
if x[j] < xmin:
f[j] = extrapolate
j+=1
continue
else:
i=0
while (xp[i] <= xval) & (i < Np-1): i+=1;
if i == (Np-1):
if x[j] != xp[i]:
f[j] = extrapolate
else:
f[j] = fp[i]
j+=1
continue
#### x[i] is now greater than xval because the
#### expression (x[i]<xval) is false, assuming
#### that xval < max(x).
# x1 = xp[i];
# x2 = xp[i+1];
# y1 = fp[i];
# y2 = fp[i+1];
x1 = xp[i-1];
x2 = xp[i];
y1 = fp[i-1];
y2 = fp[i];
out = ((y2-y1)/(x2-x1))*(xval-x1)+y1;
f[j] = out
j+=1
return f
@cython.boundscheck(False)
def interp_conserve_c(np.ndarray[DTYPE_t, ndim=1] x, np.ndarray[DTYPE_t, ndim=1] tlam, np.ndarray[DTYPE_t, ndim=1] tf, double left=0, double right=0, double integrate=0):
"""
interp_conserve_c(x, xp, fp, left=0, right=0, integrate=0)
Interpolate `xp`,`yp` array to the output x array, conserving flux.
`xp` can be irregularly spaced.
"""
cdef np.ndarray[DTYPE_t, ndim=1] templmid
cdef np.ndarray[DTYPE_t, ndim=1] tempfmid
cdef np.ndarray[DTYPE_t, ndim=1] outy
cdef unsigned long i,k,istart,ntlam,NTEMPL
cdef DTYPE_t h, numsum
# templmid = (x[1:]+x[:-1])/2. #2.+x[:-1]
# templmid = np.append(templmid, np.array([x[0], x[-1]]))
# templmid = templmid[np.argsort(templmid)]
NTEMPL = len(x)
ntlam = len(tlam)
templmid = midpoint_c(x, NTEMPL)
#tempfmid = np.interp(templmid, tlam, tf, left=left, right=right)
tempfmid = interp_c(templmid, tlam, tf, extrapolate=0.)
outy = np.zeros(NTEMPL, dtype=DTYPE)
###### Rebin template grid to master wavelength grid, conserving template flux
i=0
k=0
while templmid[k] < tlam[0]:
outy[k] = left
k+=1
if k >NTEMPL-1:
break
if(k>0) & (templmid[k-1] < tlam[0]) & (templmid[k] > tlam[0]):
m = 1;
numsum=0.;
while (tlam[m] < templmid[k]):
h = tlam[m]-tlam[m-1];
numsum+=h*(tf[m]+tf[m-1]);
m+=1;
if m >= ntlam:
break;
#print 'test #%d, %d' %(m, ntlam)
if m == 1:
h = templmid[k]-tlam[0];
numsum+=h*(tempfmid[k]+tf[0]);
else:
##### Last point
if m < ntlam:
if (templmid[k] == tlam[m]):
h = tlam[m]-tlam[m-1];
numsum+=h*(tf[m]+tf[m-1]);
else:
m-=1;
h = templmid[k]-tlam[m];
numsum+=h*(tempfmid[k]+tf[m]);
outy[k-1] = numsum*0.5;#/(templmid[k+1]-templmid[k]);
if integrate == 0.:
outy[k-1] /= (templmid[k]-templmid[k-1]);
for k in range(k, NTEMPL):
if templmid[k] > tlam[ntlam-1]:
break
numsum=0.;
#### Go to where tlam is greater than the first midpoint
while (tlam[i] < templmid[k]) & (i < ntlam): i+=1;
istart=i;
####### First point
if tlam[i] < templmid[k+1]:
h = tlam[i]-templmid[k];
numsum+=h*(tf[i]+tempfmid[k]);
i+=1;
if i==0: i+=1;
####### Template points between master grid points
while (tlam[i] < templmid[k+1]) & (i < ntlam):
h = tlam[i]-tlam[i-1];
numsum+=h*(tf[i]+tf[i-1]);
i+=1;
#### If no template points between master grid points, then just use interpolated midpoints
if i == istart:
h = templmid[k+1]-templmid[k];
numsum=h*(tempfmid[k+1]+tempfmid[k]);
else:
##### Last point
if (templmid[k+1] == tlam[i]) & (i < ntlam):
h = tlam[i]-tlam[i-1];
numsum+=h*(tf[i]+tf[i-1]);
elif (i < ntlam):
i-=1;
h = templmid[k+1]-tlam[i];
numsum+=h*(tempfmid[k+1]+tf[i]);
outy[k] = numsum*0.5;#/(templmid[k+1]-templmid[k]);
if integrate == 0.:
outy[k] /= (templmid[k+1]-templmid[k]);
return outy
def midpoint(x):
mp = (x[1:]+x[:-1])/2.
mp = np.append(mp, np.array([x[0],x[-1]]))
mp = mp[np.argsort(mp)]
return mp
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
@cython.embedsignature(True)
def midpoint_c(np.ndarray[DTYPE_t, ndim=1] x, long N):
cdef long i
cdef DTYPE_t xi,xi1
# N = len(x)
cdef np.ndarray[DTYPE_t, ndim=1] midpoint = np.zeros(N+1, dtype=DTYPE)
midpoint[0] = x[0]
midpoint[N] = x[N-1]
xi1 = x[0]
for i in range(1, N):
xi = x[i]
midpoint[i] = 0.5*xi+0.5*xi1
xi1 = xi
midpoint[0] = 2*x[0]-midpoint[1]
midpoint[N] = 2*x[N-1]-midpoint[N-1]
return midpoint
ObservationSim/MockObject/SpecDisperser/setup_c.py deleted 100644 → 0
View file @ 4afd1181
from setuptools import setup
from setuptools.extension import Extension
from setuptools.config import read_configuration
from Cython.Build import cythonize
import numpy
extensions = [
Extension("disperse_c.interp", ["disperse_c/interp.pyx"],
include_dirs = [numpy.get_include()],
libraries=["m"]),
Extension("disperse_c.disperse", ["disperse_c/disperse.pyx"],
include_dirs = [numpy.get_include()],
libraries=["m"]),
]
setup(
name = "slssim_disperse",
ext_modules = cythonize(extensions),
)
ObservationSim/MockObject/Stamp.py deleted 100644 → 0
View file @ 4afd1181
import os, sys
import numpy as np
import galsim
import astropy.constants as cons
from astropy.table import Table
from scipy import interpolate
from ObservationSim.MockObject.MockObject import MockObject
from ObservationSim.MockObject.SpecDisperser import SpecDisperser
from ObservationSim.MockObject._util import eObs, integrate_sed_bandpass, getNormFactorForSpecWithABMAG, getObservedSED, getABMAG,convolveGaussXorders
class Stamp(MockObject):
def __init__(self, param, logger=None):
super().__init__(param, logger=logger)
def unload_SED(self):
"""(Test) free up SED memory
"""
del self.sed
def drawObj_multiband(self, tel, pos_img, psf_model, bandpass_list, filt, chip, nphotons_tot=None, g1=0, g2=0, exptime=150., fd_shear=None):
if nphotons_tot == None:
nphotons_tot = self.getElectronFluxFilt(filt, tel, exptime)
try:
full = integrate_sed_bandpass(sed=self.sed, bandpass=filt.bandpass_full)
except Exception as e:
print(e)
self.logger.error(e)
return False
#nphotons_sum = 0
#photons_list = []
#xmax, ymax = 0, 0
if self.getMagFilter(filt) <= 15:
folding_threshold = 5.e-4
else:
folding_threshold = 5.e-3
gsp = galsim.GSParams(folding_threshold=folding_threshold)
self.real_pos = self.getRealPos(chip.img, global_x=self.posImg.x, global_y=self.posImg.y,
img_real_wcs=self.chip_wcs)
x, y = self.real_pos.x + 0.5, self.real_pos.y + 0.5
x_nominal = int(np.floor(x + 0.5))
y_nominal = int(np.floor(y + 0.5))
dx = x - x_nominal
dy = y - y_nominal
offset = galsim.PositionD(dx, dy)
chip_wcs_local = self.chip_wcs.local(self.real_pos)
is_updated = 0
if fd_shear:
g1 += fd_shear.g1
g2 += fd_shear.g2
gal_shear = galsim.Shear(g1=g1, g2=g2)
for i in range(len(bandpass_list)):
bandpass = bandpass_list[i]
try:
sub = integrate_sed_bandpass(sed=self.sed, bandpass=bandpass)
except Exception as e:
print(e)
if self.logger:
self.logger.error(e)
continue
ratio = sub/full
if not (ratio == -1 or (ratio != ratio)):
nphotons = ratio * nphotons_tot
else:
continue
#nphotons_sum += nphotons
psf, pos_shear = psf_model.get_PSF(chip=chip, pos_img=pos_img, bandpass=bandpass, folding_threshold=folding_threshold)
_gal = self.param['image']
galImg= galsim.ImageF(_gal, scale=self.param['pixScale'])
gal_temp= galsim.InterpolatedImage(galImg)
gal_temp= gal_temp.shear(gal_shear)
gal_temp= gal_temp.withFlux(nphotons)
gal_temp= galsim.Convolve(psf, gal_temp)
if i == 0:
gal = gal_temp
else:
gal = gal + gal_temp
stamp = gal.drawImage(wcs=chip_wcs_local, offset=offset)
if np.sum(np.isnan(stamp.array)) > 0:
# ERROR happens
return 2, pos_shear
stamp.setCenter(x_nominal, y_nominal)
bounds = stamp.bounds & galsim.BoundsI(0, chip.npix_x - 1, 0, chip.npix_y - 1)
if bounds.area() > 0:
chip.img.setOrigin(0, 0)
chip.img[bounds] += stamp[bounds]
is_updated = 1
chip.img.setOrigin(chip.bound.xmin, chip.bound.ymin)
del stamp
if is_updated == 0:
print("fits obj %s missed"%(self.id))
if self.logger:
self.logger.info("fits obj %s missed"%(self.id))
return 0, pos_shear
return 1, pos_shear
def drawObj_slitless(self, tel, pos_img, psf_model, bandpass_list, filt, chip, nphotons_tot=None, g1=0, g2=0,
exptime=150., normFilter=None, grating_split_pos=3685, fd_shear=None):
if normFilter is not None:
norm_thr_rang_ids = normFilter['SENSITIVITY'] > 0.001
sedNormFactor = getNormFactorForSpecWithABMAG(ABMag=self.param['mag_use_normal'], spectrum=self.sed,
norm_thr=normFilter,
sWave=np.floor(normFilter[norm_thr_rang_ids][0][0]),
eWave=np.ceil(normFilter[norm_thr_rang_ids][-1][0]))
if sedNormFactor == 0:
return 2, None
else:
sedNormFactor = 1.
normalSED = Table(np.array([self.sed['WAVELENGTH'], self.sed['FLUX'] * sedNormFactor]).T,
names=('WAVELENGTH', 'FLUX'))
self.real_pos = self.getRealPos(chip.img, global_x=self.posImg.x, global_y=self.posImg.y,
img_real_wcs=self.chip_wcs)
x, y = self.real_pos.x + 0.5, self.real_pos.y + 0.5
x_nominal = int(np.floor(x + 0.5))
y_nominal = int(np.floor(y + 0.5))
dx = x - x_nominal
dy = y - y_nominal
offset = galsim.PositionD(dx, dy)
chip_wcs_local = self.chip_wcs.local(self.real_pos)
if self.getMagFilter(filt) <= 15:
folding_threshold = 5.e-4
else:
folding_threshold = 5.e-3
gsp = galsim.GSParams(folding_threshold=folding_threshold)
# nphotons_sum = 0
flat_cube = chip.flat_cube
xOrderSigPlus = {'A':1.3909419820029296,'B':1.4760376591236062,'C':4.035447379743442,'D':5.5684364343742825,'E':16.260021029735388}
grating_split_pos_chip = 0 + grating_split_pos
branges = np.zeros([len(bandpass_list), 2])
# print(hasattr(psf_model, 'bandranges'))
if hasattr(psf_model, 'bandranges'):
if psf_model.bandranges is None:
return 2, None
if len(psf_model.bandranges) != len(bandpass_list):
return 2, None
branges = psf_model.bandranges
else:
for i in range(len(bandpass_list)):
branges[i, 0] = bandpass_list[i].blue_limit * 10
branges[i, 1] = bandpass_list[i].red_limit * 10
for i in range(len(bandpass_list)):
# bandpass = bandpass_list[i]
brange = branges[i]
# psf, pos_shear = psf_model.get_PSF(chip=chip, pos_img=pos_img, bandpass=bandpass, folding_threshold=folding_threshold)
_gal = self.param['image']
galImg= galsim.ImageF(_gal, scale=self.param['pixScale'])
gal = galsim.InterpolatedImage(galImg)
# (TEST) Random knots
# knots = galsim.RandomKnots(npoints=100, profile=disk)
# kfrac = np.random.random()*(1.0 - self.bfrac)
# gal = self.bfrac * bulge + (1.0 - self.bfrac - kfrac) * disk + kfrac * knots
gal = gal.withFlux(tel.pupil_area * exptime)
if fd_shear:
g1 += fd_shear.g1
g2 += fd_shear.g2
gal_shear = galsim.Shear(g1=g1, g2=g2)
gal = gal.shear(gal_shear)
# gal = galsim.Convolve(psf, gal)
# if not big_galaxy: # Not apply PSF for very big galaxy
# gal = galsim.Convolve(psf, gal)
# # if fd_shear is not None:
# # gal = gal.shear(fd_shear)
starImg = gal.drawImage(wcs=chip_wcs_local, offset=offset,method = 'real_space')
origin_star = [y_nominal - (starImg.center.y - starImg.ymin),
x_nominal - (starImg.center.x - starImg.xmin)]
starImg.setOrigin(0, 0)
gal_origin = [origin_star[0], origin_star[1]]
gal_end = [origin_star[0] + starImg.array.shape[0] - 1, origin_star[1] + starImg.array.shape[1] - 1]
if gal_origin[1] < grating_split_pos_chip < gal_end[1]:
subSlitPos = int(grating_split_pos_chip - gal_origin[1] + 1)
## part img disperse
subImg_p1 = starImg.array[:, 0:subSlitPos]
star_p1 = galsim.Image(subImg_p1)
star_p1.setOrigin(0, 0)
origin_p1 = origin_star
xcenter_p1 = min(x_nominal,grating_split_pos_chip-1) - 0
ycenter_p1 = y_nominal-0
sdp_p1 = SpecDisperser(orig_img=star_p1, xcenter=xcenter_p1,
ycenter=ycenter_p1, origin=origin_p1,
tar_spec=normalSED,
band_start=brange[0], band_end=brange[1],
conf=chip.sls_conf[0],
isAlongY=0,
flat_cube=flat_cube)
# self.addSLStoChipImage(sdp=sdp_p1, chip=chip, xOrderSigPlus = xOrderSigPlus, local_wcs=chip_wcs_local)
pos_shear = self.addSLStoChipImageWithPSF(sdp=sdp_p1, chip=chip, pos_img_local=[xcenter_p1, ycenter_p1],
psf_model=psf_model, bandNo=i + 1,
grating_split_pos=grating_split_pos,
local_wcs=chip_wcs_local, pos_img = pos_img)
subImg_p2 = starImg.array[:, subSlitPos+1:starImg.array.shape[1]]
star_p2 = galsim.Image(subImg_p2)
star_p2.setOrigin(0, 0)
origin_p2 = [origin_star[0], grating_split_pos_chip]
xcenter_p2 = max(x_nominal, grating_split_pos_chip - 1) - 0
ycenter_p2 = y_nominal - 0
sdp_p2 = SpecDisperser(orig_img=star_p2, xcenter=xcenter_p2,
ycenter=ycenter_p2, origin=origin_p2,
tar_spec=normalSED,
band_start=brange[0], band_end=brange[1],
conf=chip.sls_conf[1],
isAlongY=0,
flat_cube=flat_cube)
# self.addSLStoChipImage(sdp=sdp_p2, chip=chip, xOrderSigPlus = xOrderSigPlus, local_wcs=chip_wcs_local)
pos_shear = self.addSLStoChipImageWithPSF(sdp=sdp_p2, chip=chip, pos_img_local=[xcenter_p2, ycenter_p2],
psf_model=psf_model, bandNo=i + 1,
grating_split_pos=grating_split_pos,
local_wcs=chip_wcs_local, pos_img = pos_img)
del sdp_p1
del sdp_p2
elif grating_split_pos_chip<=gal_origin[1]:
sdp = SpecDisperser(orig_img=starImg, xcenter=x_nominal - 0,
ycenter=y_nominal - 0, origin=origin_star,
tar_spec=normalSED,
band_start=brange[0], band_end=brange[1],
conf=chip.sls_conf[1],
isAlongY=0,
flat_cube=flat_cube)
# self.addSLStoChipImage(sdp=sdp, chip=chip, xOrderSigPlus = xOrderSigPlus, local_wcs=chip_wcs_local)
pos_shear = self.addSLStoChipImageWithPSF(sdp=sdp, chip=chip, pos_img_local=[x_nominal, y_nominal],
psf_model=psf_model, bandNo=i + 1,
grating_split_pos=grating_split_pos,
local_wcs=chip_wcs_local, pos_img = pos_img)
del sdp
elif grating_split_pos_chip>=gal_end[1]:
sdp = SpecDisperser(orig_img=starImg, xcenter=x_nominal - 0,
ycenter=y_nominal - 0, origin=origin_star,
tar_spec=normalSED,
band_start=brange[0], band_end=brange[1],
conf=chip.sls_conf[0],
isAlongY=0,
flat_cube=flat_cube)
# self.addSLStoChipImage(sdp=sdp, chip=chip, xOrderSigPlus = xOrderSigPlus, local_wcs=chip_wcs_local)
pos_shear = self.addSLStoChipImageWithPSF(sdp=sdp, chip=chip, pos_img_local=[x_nominal, y_nominal],
psf_model=psf_model, bandNo=i + 1,
grating_split_pos=grating_split_pos,
local_wcs=chip_wcs_local, pos_img = pos_img)
del sdp
# print(self.y_nominal, starImg.center.y, starImg.ymin)
# del psf
return 1, pos_shear
ObservationSim/MockObject/Star.py deleted 100755 → 0
View file @ 4afd1181
import galsim
import os
import sys
import numpy as np
import astropy.constants as cons
from astropy.table import Table
from scipy import interpolate
from ObservationSim.MockObject._util import integrate_sed_bandpass, getNormFactorForSpecWithABMAG, getObservedSED, getABMAG, tag_sed
from ObservationSim.MockObject.MockObject import MockObject
class Star(MockObject):
def __init__(self, param, logger=None):
super().__init__(param, logger=logger)
if not hasattr(self, "mu"):
self.mu = 1.
def unload_SED(self):
"""(Test) free up SED memory
"""
del self.sed
def getGSObj(self, psf, g1=0, g2=0, flux=None, filt=None, tel=None, exptime=150.):
if flux == None:
flux = self.getElectronFluxFilt(filt, tel, exptime)
# star = galsim.Gaussian(sigma=1.e-8, flux=1.)
star = galsim.DeltaFunction()
star = star.withFlux(flux)
final = galsim.Convolve(psf, star)
return final
def getGSObj_multiband(self, tel, psf_list, bandpass_list, filt, nphotons_tot=None, g1=0, g2=0, exptime=150.):
if len(psf_list) != len(bandpass_list):
raise ValueError(
"!!!The number of PSF profiles and the number of bandpasses must be equal.")
objs = []
if nphotons_tot == None:
nphotons_tot = self.getElectronFluxFilt(filt, tel, exptime)
try:
full = integrate_sed_bandpass(
sed=self.sed, bandpass=filt.bandpass_full)
except Exception as e:
print(e)
self.logger.error(e)
return -1
for i in range(len(bandpass_list)):
bandpass = bandpass_list[i]
psf = psf_list[i]
try:
sub = integrate_sed_bandpass(sed=self.sed, bandpass=bandpass)
except Exception as e:
print(e)
self.logger.error(e)
return -1
ratio = sub/full
if not (ratio == -1 or (ratio != ratio)):
nphotons = ratio * nphotons_tot
else:
return -1
star = galsim.DeltaFunction()
star = star.withFlux(nphotons)
star = galsim.Convolve(psf, star)
objs.append(star)
final = galsim.Sum(objs)
return final
ObservationSim/MockObject/__init__.py deleted 100755 → 0
View file @ 4afd1181
from .MockObject import MockObject
from .Galaxy import Galaxy
from .CatalogBase import CatalogBase
from .Quasar import Quasar
from .Star import Star
from .Stamp import Stamp
from .FlatLED import FlatLED
# from .SkybackgroundMap import *
# from .CosmicRay import CosmicRay
ObservationSim/MockObject/_util.py deleted 100755 → 0
View file @ 4afd1181
import numpy as np
import os
from scipy.interpolate import InterpolatedUnivariateSpline, interp1d
from scipy import interpolate, integrate
from astropy.table import Table
import galsim
VC_A = 2.99792458e+18 # speed of light: A/s
VC_M = 2.99792458e+8 # speed of light: m/s
H_PLANK = 6.626196e-27 # Plank constant: erg s
def comoving_dist(z, om_m=0.3111, om_L=0.6889, h=0.6766):
# Return comving distance in pc
H0 = h*100. # km / (s Mpc)
def dist_int(z):
return 1./np.sqrt(om_m*(1.+z)**3 + om_L)
res, err = integrate.quad(dist_int, 0., z)
return [res * (VC_M/1e3/H0) * 1e6, err * (VC_M/1e3/H0) * 1e6]
def magToFlux(mag):
"""
flux of a given AB magnitude
Parameters:
mag: magnitude in unit of AB
Return:
flux: flux in unit of erg/s/cm^2/Hz
"""
flux = 10**(-0.4*(mag+48.6))
return flux
def extAv(nav, seed=1212123):
"""
Generate random intrinsic extinction Av
following the distribution from Holwerda et al, 2015
"""
np.random.seed(seed)
tau = 0.4
peak, a = 0.1, 0.5
b = a*(tau-peak)
pav = lambda av: (a*av+b)*np.exp(-av/tau)
avmin, avmax = 0., 3.
avs = np.linspace(avmin, avmax, int((avmax-avmin)/0.001)+1)
norm = np.trapz(pav(avs), avs)
pav_base = pav(avs)/np.sum(pav(avs))
rav = np.random.choice(avs, nav, p=pav_base)
return rav
###########################################
def seds(sedlistn, seddir="./", unit="A"):
"""
read SEDs and save into Python dictionary
Parameters:
sedlistn: filename of the sed template list and corresponding intrinsic
extinction model, see tmag_dz.py for detailes
listdir: directory of the list
unit: wavelength unit of the input templates
Return:
SED dictionary, the output wavelength unit is 'A'
"""
seds = {}
reds = {}
sedlist = seddir + sedlistn
sedn = open(sedlist).read().splitlines()
sedtype = range(1,len(sedn)+1)
for i in range(len(sedn)):
xxx = sedn[i].split()
isedn = seddir+xxx[0]
itype = sedtype[i]
ised = np.loadtxt(isedn)
if unit=="nm": ised[:,0] *= 10.0
seds[itype] = ised
reds[itype] = int(xxx[1])
return seds, reds
def sed_assign(phz, btt, rng):
"""
assign SED template to a galaxy.
"""
sedid = list(range(1, 34))
lzid = sedid[0:13] + sedid[23:28]
hzid = sedid[13:23]
if btt == 1.0:
sedtype = rng.sample(sedid[0:5] + sedid[28:33], 1)[0]
if phz > 2.0 and sedtype in sedid[0:5]:
sedtype = rng.sample(sedid[28:33], 1)[0]
elif btt > 0.3 and btt < 1.0:
sedtype = rng.sample(sedid, 1)[0]
if phz > 2.0 and sedtype in lzid:
sedtype = rng.sample(hzid, 1)[0]
elif btt >= 0.1 and btt < 0.3:
sedtype = rng.sample(sedid[5:28], 1)[0]
if phz > 1.5 and sedtype in lzid:
sedtype = rng.sample(hzid, 1)[0]
elif btt >= 0.0 and btt < 0.1:
sedtype = rng.sample(sedid[5:23], 1)[0]
if phz > 1.5 and sedtype in lzid:
sedtype = rng.sample(hzid, 1)[0]
else:
sedtype = 0
return sedtype
###########################################
def tflux(filt, sed, redshift=0.0, av=0.0, redden=0):
"""
calculate the theoretical SED for given filter set and template
Only AB magnitude is support!!!
Parameters:
filt: 2d array
fliter transmissions: lambda(A), T
sed: 2d array
sed templateL lambda (A), flux (erg s-1 cm-2 A-1)
redshift: float
redshift of the corresponding source, default is zero
av: float
extinction value for intrincal extinction
redden: int
reddening model, see Function 'reddening' for details
return:
tflux: float
theoretical flux
sedObs: array
SED in observed frame
"""
z = redshift + 1.0
sw, sf = sed[:,0], sed[:,1]
# reddening
sf = reddening(sw, sf, av=av, model=redden)
sw, sf = sw*z, sf*(z**3)
# lyman forest correction
sf = lyman_forest(sw, sf, redshift)
sedxx = (sw.copy(), sf.copy())
sw = VC_A/sw
sf = sf*(VC_A/sw**2) # convert flux unit to erg/s/cm^s/Hz
sw, sf = sw[::-1], sf[::-1]
sfun = interp1d(sw, sf, kind='linear')
fwave, fresp = filt[:,0], filt[:,1]
fwave = VC_A/fwave
fwave, fresp = fwave[::-1], fresp[::-1]
tflux = sfun(fwave)
zpflux = 3.631*1.0e-20
tflux = np.trapz(tflux*fresp/fwave,fwave)/np.trapz(zpflux*fresp/fwave,fwave)
#tflux = np.trapz(tflux*fresp,fwave)/np.trapz(zpflux*fresp,fwave)
return tflux, sedxx
###########################################
def lyman_forest(wavelen, flux, z):
"""
Compute the Lyman forest mean absorption of an input spectrum,
according to D_A and D_B evolution from Madau (1995).
The waveeln and flux are in observed frame
"""
if z<=0:
flux0 = flux
else:
nw = 200
istep = np.linspace(0,nw-1,nw)
w1a, w2a = 1050.0*(1.0+z), 1170.0*(1.0+z)
w1b, w2b = 920.0*(1.0+z), 1015.0*(1.0+z)
wstepa = (w2a-w1a)/float(nw)
wstepb = (w2b-w1b)/float(nw)
wtempa = w1a + istep*wstepa
ptaua = np.exp(-3.6e-03*(wtempa/1216.0)**3.46)
wtempb = w1b + istep*wstepb
ptaub = np.exp(-1.7e-3*(wtempb/1026.0)**3.46\
-1.2e-3*(wtempb/972.50)**3.46\
-9.3e-4*(wtempb/950.00)**3.46)
da = (1.0/(120.0*(1.0+z)))*np.trapz(ptaua, wtempa)
db = (1.0/(95.0*(1.0+z)))*np.trapz(ptaub, wtempb)
if da>1.0: da=1.0
if db>1.0: db=1.0
if da<0.0: da=0.0
if db<0.0: db=0.0
flux0 = flux.copy()
id0 = wavelen<=1026.0*(1.0+z)
id1 = np.logical_and(wavelen<1216.0*(1.0+z),wavelen>=1026.0*(1.0+z))
flux0[id0] = db*flux[id0]
flux0[id1] = da*flux[id1]
return flux0
###########################################
def reddening(sw, sf, av=0.0, model=0):
"""
calculate the intrinsic extinction of a given template
Parameters:
sw: array
wavelength
sf: array
flux
av: float or array
model: int
Five models will be used:
1: Allen (1976) for the Milky Way
2: Seaton (1979) fit by Fitzpatrick (1986) for the Milky Way
3: Fitzpatrick (1986) for the Large Magellanic Cloud (LMC)
4: Prevot et al (1984) and Bouchet (1985) for the Small Magellanic Cloud (SMC)
5: Calzetti et al (2000) for starburst galaxies
6: Reddy et al (2015) for star forming galaxies
Return:
reddening-corrected flux or observed flux
"""
if model==0 or av==0.0:
flux=sf
elif model==1: # Allen (1976) for the Milky Way
lambda0 = np.array([1000, 1110, 1250, 1430, 1670, \
2000, 2220, 2500, 2850, 3330, \
3650, 4000, 4400, 5000, 5530, \
6700, 9000, 10000, 20000, 100000], dtype=float)
kR = np.array([4.20, 3.70, 3.30, 3.00, 2.70, \
2.80, 2.90, 2.30, 1.97, 1.69, \
1.58, 1.45, 1.32, 1.13, 1.00, \
0.74, 0.46, 0.38, 0.11, 0.00],dtype=float)
ext0 = InterpolatedUnivariateSpline(lambda0, kR, k=1)
A_lambda = av*ext0(sw)
A_lambda[A_lambda<0.0] = 0.0
flux = sf*10**(-0.4*A_lambda)
elif model==2: # Seaton (1979) fit by Fitzpatrick (1986) for the Milky Way
Rv=3.1
al0, ga, c1, c2, c3, c4 = 4.595, 1.051, -0.38, 0.74, 3.96, 0.26
ff11 = __red(1100.0,al0,ga,c1,c2,c3,c4)
ff12 = __red(1200.0,al0,ga,c1,c2,c3,c4)
slope=(ff12-ff11)/100.0
lambda0 = np.array([3650, 4000, 4400, 5000, 5530, \
6700, 9000, 10000, 20000, 100000], dtype=float)
kR = np.array([1.58, 1.45, 1.32, 1.13, 1.00, \
0.74, 0.46, 0.38, 0.11, 0.00],dtype=float)
fun = interp1d(lambda0, kR, kind='linear')
sw0 = sw[sw<1200.0]
A_lambda0 = (ff11+(sw0-1100.0)*slope)/Rv+1.0
sw1 = sw[np.logical_and(sw>=1200.0, sw<=3650.0)]
ff = __red(sw1,al0,ga,c1,c2,c3,c4)
A_lambda1 = ff/Rv+1.0
sw2 = sw[np.logical_and(sw>3650.0, sw<=100000.0)]
A_lambda2 = fun(sw2)
A_lambda3 = sw[sw>100000.0]*0.0
A_lambda = av*np.hstack([A_lambda0,A_lambda1,A_lambda2,A_lambda3])
A_lambda[A_lambda<0.0] = 0.0
flux = sf*10**(-0.4*A_lambda)
elif model==3: # Fitzpatrick (1986) for the Large Magellanic Cloud (LMC)
Rv=3.1
al0, ga, c1, c2, c3, c4 = 4.608, 0.994, -0.69, 0.89, 2.55, 0.50
ff11 = __red(1100.0,al0,ga,c1,c2,c3,c4)
ff12 = __red(1200.0,al0,ga,c1,c2,c3,c4)
slope=(ff12-ff11)/100.0
lambda0 = np.array([3330, 3650, 4000, 4400, 5000, 5530, \
6700, 9000, 10000, 20000, 100000], dtype=float)
kR = np.array([1.682, 1.58, 1.45, 1.32, 1.13, 1.00, \
0.74, 0.46, 0.38, 0.11, 0.00],dtype=float)
fun = interp1d(lambda0, kR, kind='linear')
sw0 = sw[sw<1200.0]
A_lambda0 = (ff11+(sw0-1100.0)*slope)/Rv+1.0
sw1 = sw[np.logical_and(sw>=1200.0, sw<=3330.0)]
ff = __red(sw1,al0,ga,c1,c2,c3,c4)
A_lambda1 = ff/Rv+1.0
sw2 = sw[np.logical_and(sw>3330.0, sw<=100000.0)]
A_lambda2 = fun(sw2)
A_lambda3 = sw[sw>100000.0]*0.0
A_lambda = av*np.hstack([A_lambda0,A_lambda1,A_lambda2,A_lambda3])
A_lambda[A_lambda<0.0] = 0.0
flux = sf*10**(-0.4*A_lambda)
elif model==4: # Prevot et al (1984) and Bouchet (1985) for the Small Magellanic Cloud (SMC)
Rv = 2.72
lambda0 = np.array([1275, 1330, 1385, 1435, 1490, 1545, \
1595, 1647, 1700, 1755, 1810, 1860, \
1910, 2000, 2115, 2220, 2335, 2445, \
2550, 2665, 2778, 2890, 2995, 3105, \
3704, 4255, 5291, 12500, 16500, 22000], dtype=float)
kR = np.array([13.54, 12.52, 11.51, 10.80, 9.84, 9.28, \
9.06, 8.49, 8.01, 7.71, 7.17, 6.90, 6.76, \
6.38, 5.85, 5.30, 4.53, 4.24, 3.91, 3.49, \
3.15, 3.00, 2.65, 2.29, 1.81, 1.00, 0.00, \
-2.02, -2.36, -2.47],dtype=float)
kR = kR/Rv+1.0
ext0 = InterpolatedUnivariateSpline(lambda0, kR, k=1)
A_lambda = av*ext0(sw)
A_lambda[A_lambda<0.0] = 0.0
flux = sf*10**(-0.4*A_lambda)
elif model==5: # Calzetti et al (2000) for starburst galaxies
Rv = 4.05
sw = sw*1.0e-04 #wavelength in microns
fun1 = lambda x: 2.659*(-2.156+1.509/x-0.198/x**2+0.011/x**3)+Rv
fun2 = lambda x: 2.659*(-1.857+1.040/x)+Rv
ff11, ff12 = fun1(0.11), fun1(0.12)
slope1=(ff12-ff11)/0.01
ff99, ff100 = fun2(2.19), fun2(2.2)
slope2=(ff100-ff99)/0.01
sw0 = sw[sw<0.12]
sw1 = sw[np.logical_and(sw>=0.12, sw<=0.63)]
sw2 = sw[np.logical_and(sw>0.63, sw<=2.2)]
sw3 = sw[sw>2.2]
k_lambda0 = ff11+(sw0-0.11)*slope1
k_lambda1, k_lambda2 = fun1(sw1), fun2(sw2)
k_lambda3 = ff99+(sw3-2.19)*slope2
A_lambda = av*np.hstack([k_lambda0,k_lambda1,k_lambda2,k_lambda3])/Rv
A_lambda[A_lambda<0.0] = 0.0
flux = sf*10**(-0.4*A_lambda)
elif model==6: # Reddy et al (2015) for satr forming galaxies
Rv = 2.505
sw = sw*1.0e-04
fun1 = lambda x: -5.726+4.004/x-0.525/x**2+0.029/x**3+Rv
fun2 = lambda x: -2.672-0.010/x+1.532/x**2-0.412/x**3+Rv
ff11, ff12 = fun1(0.14), fun1(0.15)
slope1=(ff12-ff11)/0.01
ff99, ff100 = fun2(2.84), fun2(2.85)
slope2=(ff100-ff99)/0.01
sw0 = sw[sw<0.15]
sw1 = sw[np.logical_and(sw>=0.15, sw<0.60)]
sw2 = sw[np.logical_and(sw>=0.60, sw<2.85)]
sw3 = sw[sw>=2.85]
k_lambda0 = ff11+(sw0-0.14)*slope1
k_lambda1, k_lambda2 = fun1(sw1), fun2(sw2)
k_lambda3 = ff99+(sw3-2.84)*slope2
A_lambda = av*np.hstack([k_lambda0,k_lambda1,k_lambda2,k_lambda3])/Rv
A_lambda[A_lambda<0.0] = 0.0
flux = sf*10**(-0.4*A_lambda)
else:
raise ValueError("!!! Please select a proper reddening model")
return flux
###########################################
def __red(alan,al0,ga,c1,c2,c3,c4):
fun1 = lambda x: c3/(((x-(al0**2/x))**2)+ga*ga)
fun2 = lambda x,cc: cc*(0.539*((x-5.9)**2)+0.0564*((x-5.9)**3))
fun = lambda x,cc: c1+c2*x+fun1(x)+fun2(x,cc)
ala = alan*1.0e-04 #wavelength in microns
p = 1.0/ala
if np.size(p)>1:
p1, p2 = p[p>=5.9], p[p<5.9]
ff = np.append(fun(p1,c4), fun(p2,0.0))
elif np.size(p)==1:
if p<5.9: c4 = 0.0
ff = fun(p, c4)
else:
return
return ff
###########################################
def sed2mag(mag_i, sedCat, filter_list, redshift=0.0, av=0.0, redden=0):
# load the filters
nfilt = len(filter_list)
aflux = np.zeros(nfilt)
nid = -1
for k in range(nfilt):
if filter_list[k].filter_type == 'i':
nid = k
bandpass = filter_list[k].bandpass_full
ktrans = np.transpose(np.array([bandpass.wave_list*10.0, bandpass.func(bandpass.wave_list)]))
aflux[k], isedObs = tflux(ktrans, sedCat, redshift=redshift, av=av, redden=redden)
# normalize to i-band
aflux = aflux / aflux[nid]
# magnitudes in all filters
amag = -2.5*np.log10(aflux) + mag_i
spec = galsim.LookupTable(x=np.array(isedObs[0]), f=np.array(isedObs[1]), interpolant='nearest')
isedObs = galsim.SED(spec, wave_type='A', flux_type='1', fast=False)
return amag, isedObs
def eObs(e1,e2,g1,g2):
"""
Calculate the sheared (observed) ellipticity using the
intrinsic ellipticity and cosmic shear components.
Parameters:
e1, e2: scalar or numpy array
g1, g2: scalar or numpy array
Return:
Sheared (observed) ellipticity components, absolute value, and orientation
in format of scalar or numpy array
** NOTE: e1, e2, g1, and g2 should have the same dimensions.
"""
if np.isscalar(e1):
e1 = np.array([e1])
e2 = np.array([e2])
g1 = np.array([g1])
g2 = np.array([g2])
else:
e1 = e1.flatten()
e2 = e2.flatten()
g1 = g1.flatten()
g2 = g2.flatten()
# calculate the sheared (observed) ellipticity using complex rule
nobj = len(e1)
e1obs, e2obs = [], []
eeobs, theta = [], []
for i in range(nobj):
e = complex(e1[i], e2[i])
g = complex(g1[i], g2[i])
e, gg = abs(e), abs(g)
if gg<=1.0:
tt = e + g
bb = 1.0 + e*g.conjugate()
eobs = tt/bb
else:
tt = 1.0 + g*e.conjugate()
bb = e.conjugate() + g.conjugate()
eobs = tt/bb
# derive the orientation
dd = 0.5*np.arctan(abs(eobs.imag/eobs.real))*180.0/np.pi
if eobs.imag>0 and eobs.real>0: dd = dd
if eobs.imag>0 and eobs.real<0: dd = 90.0 - dd
if eobs.imag<0 and eobs.real>0: dd = 0.0 - dd
if eobs.imag<0 and eobs.real<0: dd = dd - 90.0
e1obs += [eobs.real]
e2obs += [eobs.imag]
eeobs += [abs(eobs)]
theta += [dd]
e1obs,e2obs,eeobs,theta = np.array(e1obs),np.array(e2obs),np.array(eeobs),np.array(theta)
if nobj == 1: e1obs,e2obs,eeobs,theta = e1obs[0],e2obs[0],eeobs[0],theta[0]
return e1obs, e2obs, eeobs, theta
def getObservedSED(sedCat, redshift=0.0, av=0.0, redden=0):
z = redshift + 1.0
sw, sf = sedCat[:,0], sedCat[:,1]
# reddening
sf = reddening(sw, sf, av=av, model=redden)
# sw, sf = sw*z, sf*(z**3)
sw, sf = sw*z, sf/z
# sw, sf = sw*z, sf
# lyman forest correction
sf = lyman_forest(sw, sf, redshift)
isedObs = (sw.copy(), sf.copy())
return isedObs
def integrate_sed_bandpass(sed, bandpass):
wave = np.linspace(bandpass.blue_limit, bandpass.red_limit, 1000) # in nm
flux_normalized = sed(wave)*bandpass(wave)
# print('in integrate_sed_bandpass', bandpass.blue_limit, bandpass.red_limit)
int_flux = np.trapz(y=flux_normalized, x=wave) * 10. # convert to photons s-1 m-2 A-1
return int_flux
def getABMAG(interFlux, bandpass):
throughtput = Table(np.array(np.array([bandpass.wave_list*10.0, bandpass.func(bandpass.wave_list)])).T, names=(['WAVELENGTH', 'SENSITIVITY']))
sWave = bandpass.blue_limit*10.0
eWave = bandpass.red_limit*10.0
# print('in getABMAG', sWave, eWave)
ABMAG_zero = getABMagAverageVal(
ABmag=0,
norm_thr=throughtput,
sWave=sWave,
eWave=eWave)
flux_ave = interFlux / (eWave-sWave)
ABMAG_spec = -2.5 * np.log10(flux_ave/ABMAG_zero)
return ABMAG_spec
def getABMagAverageVal(ABmag=20.,norm_thr=None, sWave=6840, eWave=8250):
"""
norm_thr: astropy.table, 2 colum, 'WAVELENGTH', 'SENSITIVITY'
Return:
the integerate flux of AB magnitude in the norm_thr(the throughtput of band), photos s-1 m-2 A-1
"""
inverseLambda = norm_thr['SENSITIVITY']/norm_thr['WAVELENGTH']
norm_thr_i = interpolate.interp1d(norm_thr['WAVELENGTH'], inverseLambda)
x = np.linspace(sWave,eWave, int(eWave)-int(sWave)+1)
y = norm_thr_i(x)
AverageLamdaInverse = np.trapz(y,x)/(eWave-sWave)
norm = 54798696332.52474 * pow(10.0, -0.4 * ABmag) * AverageLamdaInverse
# print('AverageLamdaInverse = ', AverageLamdaInverse)
# print('norm = ', norm)
return norm
def getNormFactorForSpecWithABMAG(ABMag=20., spectrum=None, norm_thr=None, sWave=6840, eWave=8250):
"""
Use AB magnitude system (zero point, fv = 3631 janskys) in the normal band(norm_thr) normalize the spectrum by inpute ABMag
Parameters
----------
spectrum: astropy.table, 2 colum, 'WAVELENGTH', 'FLUX'
norm_thr: astropy.table, 2 colum, 'WAVELENGTH', 'SENSITIVITY'
sWave: the start of norm_thr
eWave: the end of norm_thr
Return:
the normalization factor flux of AB system(fix a band and magnitude ) /the flux of inpute spectrum(fix a band)
"""
spectrumi = interpolate.interp1d(spectrum['WAVELENGTH'], spectrum['FLUX'])
norm_thri = interpolate.interp1d(norm_thr['WAVELENGTH'], norm_thr['SENSITIVITY'])
x = np.linspace(sWave,eWave, int(eWave)-int(sWave)+1)
y_spec = spectrumi(x)
y_thr = norm_thri(x)
y = y_spec*y_thr
specAve = np.trapz(y,x)/(eWave-sWave)
norm = getABMagAverageVal(ABmag=ABMag, norm_thr=norm_thr, sWave=sWave, eWave=eWave)
if specAve == 0:
return 0
# print('specAve = ', specAve)
return norm / specAve
def tag_sed(h5file, model_tag, teff=5000, logg=2, feh=0):
model_tag_str = model_tag.decode("utf-8").strip()
teff_grid = np.unique(h5file["teff"][model_tag_str])
logg_grid = np.unique(h5file["logg"][model_tag_str])
feh_grid = np.unique(h5file["feh"][model_tag_str])
close_teff = teff_grid[np.argmin(np.abs(teff_grid - teff))]
close_logg = logg_grid[np.argmin(np.abs(logg_grid - logg))]
if model_tag_str == "koester" or model_tag_str == "MC":
close_feh = 99
else:
close_feh = feh_grid[np.argmin(np.abs(feh_grid - feh))]
path = model_tag_str + f"_teff_{close_teff:.1f}_logg_{close_logg:.2f}_feh_{close_feh:.1f}"
wave = np.array(h5file["wave"][model_tag_str][()]).ravel()
flux = np.array(h5file["sed"][path][()]).ravel()
return path, wave, flux
def convolveGaussXorders(img=None, sigma = 1):
from astropy.modeling.models import Gaussian2D
from scipy import signal
offset = int(np.ceil(sigma*10))
g_size = 2*offset+1
m_cen = int(g_size/2)
g_PSF_ = Gaussian2D(1, m_cen, m_cen, sigma, sigma)
yp, xp = np.mgrid[0:g_size, 0:g_size]
g_PSF = g_PSF_(xp, yp)
psf = g_PSF/g_PSF.sum()
convImg = signal.fftconvolve(img, psf, mode='full', axes=None)
return convImg, offset
def convolveImg(img=None, psf = None):
from astropy.modeling.models import Gaussian2D
from scipy import signal
convImg = signal.fftconvolve(img, psf, mode='full', axes=None)
offset_x = int(psf.shape[1]/2. + 0.5) - 1
offset_y = int(psf.shape[0]/2. + 0.5) - 1
offset = [offset_x,offset_y]
return convImg, offset
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment