Skip to content
GitLab
Commit 64eb7956 authored by Fang Yuedong's avatar Fang Yuedong
Browse files

Revert "Merge branch 'milky_way_extinction' into 'release_v3.0'" 

This reverts merge request !28
parent 914004c1
Show whitespace changes
Inline Side-by-side
catalog/C9_Catalog.py
View file @ 64eb7956
...@@ -12,7 +12,7 @@ from astropy.table import Table ...@@ -12,7 +12,7 @@ from astropy.table import Table
from scipy import interpolate from scipy import interpolate
from datetime import datetime from datetime import datetime
from observation_sim.mock_objects import CatalogBase, Star, Galaxy, Quasar, ExtinctionMW from observation_sim.mock_objects import CatalogBase, Star, Galaxy, Quasar
from observation_sim.mock_objects._util import tag_sed, getObservedSED, getABMAG, integrate_sed_bandpass, comoving_dist from observation_sim.mock_objects._util import tag_sed, getObservedSED, getABMAG, integrate_sed_bandpass, comoving_dist
from observation_sim.astrometry.Astrometry_util import on_orbit_obs_position from observation_sim.astrometry.Astrometry_util import on_orbit_obs_position
...@@ -110,16 +110,6 @@ class Catalog(CatalogBase): ...@@ -110,16 +110,6 @@ class Catalog(CatalogBase):
self.max_size = 0. self.max_size = 0.
# [TODO] Milky Way extinction
if "enable_mw_ext_gal" in config["catalog_options"] and config["catalog_options"]["enable_mw_ext_gal"]:
if "planck_ebv_map" not in config["catalog_options"]:
raise ValueError(
"Planck dust map must be given to enable Milky Way extinction calculation for galaxies.")
self.mw_ext = ExtinctionMW()
self.mw_ext.init_ext_model(model_name="odonnell")
self.mw_ext.load_Planck_ext(
file_path=config["catalog_options"]["planck_ebv_map"])
if "star_cat" in config["catalog_options"]["input_path"] and config["catalog_options"]["input_path"]["star_cat"] and not config["catalog_options"]["galaxy_only"]: if "star_cat" in config["catalog_options"]["input_path"] and config["catalog_options"]["input_path"]["star_cat"] and not config["catalog_options"]["galaxy_only"]:
# Get the cloest star catalog file # Get the cloest star catalog file
star_file_name = get_star_cat( star_file_name = get_star_cat(
...@@ -258,12 +248,6 @@ class Catalog(CatalogBase): ...@@ -258,12 +248,6 @@ class Catalog(CatalogBase):
input_time_str=time_str input_time_str=time_str
) )
# [TODO] get Milky Way extinction AVs
if "enable_mw_ext_gal" in self.config["catalog_options"] and self.config["catalog_options"]["enable_mw_ext_gal"]:
MW_Av_arr = self.mw_ext.Av_from_Planck(ra=ra_arr, dec=dec_arr)
else:
MW_Av_arr = np.zeros(len(ra_arr))
for igals in range(ngals): for igals in range(ngals):
# # (TEST) # # (TEST)
# if igals > 100: # if igals > 100:
...@@ -275,9 +259,6 @@ class Catalog(CatalogBase): ...@@ -275,9 +259,6 @@ class Catalog(CatalogBase):
param['ra_orig'] = gals['ra'][igals] param['ra_orig'] = gals['ra'][igals]
param['dec_orig'] = gals['dec'][igals] param['dec_orig'] = gals['dec'][igals]
# [TODO] get Milky Way extinction AVs
param['mw_Av'] = MW_Av_arr[igals]
if not self.chip.isContainObj(ra_obj=param['ra'], dec_obj=param['dec'], margin=200): if not self.chip.isContainObj(ra_obj=param['ra'], dec_obj=param['dec'], margin=200):
continue continue
...@@ -549,30 +530,11 @@ class Catalog(CatalogBase): ...@@ -549,30 +530,11 @@ class Catalog(CatalogBase):
speci = interpolate.interp1d(wave, flux) speci = interpolate.interp1d(wave, flux)
lamb = np.arange(2000, 11001+0.5, 0.5) lamb = np.arange(2000, 11001+0.5, 0.5)
y = speci(lamb) y = speci(lamb)
# [TODO] Apply Milky Way extinction
if obj.type != 'star' and ("enable_mw_ext_gal" in self.config["catalog_options"] and self.config["catalog_options"]["enable_mw_ext_gal"]):
y = self.mw_ext.apply_extinction(y, Av=obj.mw_Av)
# erg/s/cm2/A --> photon/s/m2/A # erg/s/cm2/A --> photon/s/m2/A
all_sed = y * lamb / (cons.h.value * cons.c.value) * 1e-13 all_sed = y * lamb / (cons.h.value * cons.c.value) * 1e-13
sed = Table(np.array([lamb, all_sed]).T, names=('WAVELENGTH', 'FLUX')) sed = Table(np.array([lamb, all_sed]).T, names=('WAVELENGTH', 'FLUX'))
# if obj.type == 'quasar': if obj.type == 'quasar':
# # integrate to get the magnitudes
# sed_photon = np.array([sed['WAVELENGTH'], sed['FLUX']]).T
# sed_photon = galsim.LookupTable(x=np.array(sed_photon[:, 0]), f=np.array(
# sed_photon[:, 1]), interpolant='nearest')
# sed_photon = galsim.SED(
# sed_photon, wave_type='A', flux_type='1', fast=False)
# interFlux = integrate_sed_bandpass(
# sed=sed_photon, bandpass=self.filt.bandpass_full)
# obj.param['mag_use_normal'] = getABMAG(
# interFlux, self.filt.bandpass_full)
# # mag = getABMAG(interFlux, self.filt.bandpass_full)
# # print("mag diff = %.3f"%(mag - obj.param['mag_use_normal']))
# integrate to get the magnitudes # integrate to get the magnitudes
if obj.type == 'quasar' or obj.type == 'galaxy':
sed_photon = np.array([sed['WAVELENGTH'], sed['FLUX']]).T sed_photon = np.array([sed['WAVELENGTH'], sed['FLUX']]).T
sed_photon = galsim.LookupTable(x=np.array(sed_photon[:, 0]), f=np.array( sed_photon = galsim.LookupTable(x=np.array(sed_photon[:, 0]), f=np.array(
sed_photon[:, 1]), interpolant='nearest') sed_photon[:, 1]), interpolant='nearest')
......
config/config_overall.yaml
View file @ 64eb7956
...@@ -11,7 +11,7 @@ ...@@ -11,7 +11,7 @@
# can add some of the command-line arguments here as well; # can add some of the command-line arguments here as well;
# ok to pass either way or both, as long as they are consistent # ok to pass either way or both, as long as they are consistent
work_dir: "/public/home/fangyuedong/project/workplace/" work_dir: "/public/home/fangyuedong/project/workplace/"
run_name: "ext_on" run_name: "no_nonlinearity_test"
# Project cycle and run counter are used to name the outputs # Project cycle and run counter are used to name the outputs
project_cycle: 9 project_cycle: 9
...@@ -49,10 +49,6 @@ catalog_options: ...@@ -49,10 +49,6 @@ catalog_options:
# rotate galaxy ellipticity # rotate galaxy ellipticity
rotateEll: 0. # [degree] rotateEll: 0. # [degree]
# Whether to apply milky way extinction to galaxies
enable_mw_ext_gal: YES
planck_ebv_map: "/public/home/fangyuedong/project/ext_maps/planck/HFI_CompMap_ThermalDustModel_2048_R1.20.fits"
############################################### ###############################################
# Observation setting # Observation setting
############################################### ###############################################
...@@ -69,10 +65,10 @@ obs_setting: ...@@ -69,10 +65,10 @@ obs_setting:
# - give a list of indexes of pointings: [ip_1, ip_2...] # - give a list of indexes of pointings: [ip_1, ip_2...]
# - run all pointings: null # - run all pointings: null
# Note: only valid when a pointing list is specified # Note: only valid when a pointing list is specified
run_pointings: [0, 1, 2, 3, 4] run_pointings: [1]
# Whether to enable astrometric modeling # Whether to enable astrometric modeling
enable_astrometric_model: YES enable_astrometric_model: True
# Cut by saturation magnitude in which band? # Cut by saturation magnitude in which band?
cut_in_band: "z" cut_in_band: "z"
......
config/obs_config_SCI.yaml
View file @ 64eb7956
...@@ -16,7 +16,7 @@ obs_id: "00000001" # this setting will only be used if pointing list file is not ...@@ -16,7 +16,7 @@ obs_id: "00000001" # this setting will only be used if pointing list file is not
# Define list of chips # Define list of chips
# run_chips: [6,7,8,9,11,12,13,14,15,16,17,18,19,20,22,23,24,25] # Photometric chips # run_chips: [6,7,8,9,11,12,13,14,15,16,17,18,19,20,22,23,24,25] # Photometric chips
#run_chips: [1,2,3,4,5,10,21,26,27,28,29,30] # Spectroscopic chips #run_chips: [1,2,3,4,5,10,21,26,27,28,29,30] # Spectroscopic chips
run_chips: [17, 22] run_chips: [1]
# Define observation sequence # Define observation sequence
call_sequence: call_sequence:
......
observation_sim/ObservationSim.py
View file @ 64eb7956
...@@ -23,7 +23,7 @@ class Observation(object): ...@@ -23,7 +23,7 @@ class Observation(object):
self.filter_param = FilterParam() self.filter_param = FilterParam()
self.Catalog = Catalog self.Catalog = Catalog
def prepare_chip_for_exposure(self, chip, ra_cen, dec_cen, pointing, wcs_fp=None, slsPSFOptim=False): def prepare_chip_for_exposure(self, chip, ra_cen, dec_cen, pointing, wcs_fp=None):
# Get WCS for the focal plane # Get WCS for the focal plane
if wcs_fp == None: if wcs_fp == None:
wcs_fp = self.focal_plane.getTanWCS( wcs_fp = self.focal_plane.getTanWCS(
...@@ -34,27 +34,6 @@ class Observation(object): ...@@ -34,27 +34,6 @@ class Observation(object):
chip.img.setOrigin(chip.bound.xmin, chip.bound.ymin) chip.img.setOrigin(chip.bound.xmin, chip.bound.ymin)
chip.img.wcs = wcs_fp chip.img.wcs = wcs_fp
chip.slsPSFOptim = slsPSFOptim
if chip.chipID in [1, 2, 3, 4, 5, 10, 21, 26, 27, 28, 29, 30] and slsPSFOptim:
chip.img_stack = {}
for id1 in np.arange(2):
gn = chip_utils.getChipSLSGratingID(chip.chipID)[id1]
orders = {}
# for id2 in ['-2','-1','0','1','2']:
for id2 in ['0', '1']:
o_n = "order"+id2
allbands = {}
for id3 in ['1', '2', '3', '4']:
w_n = "w"+id3
allbands[w_n] = galsim.ImageF(chip.npix_x, chip.npix_y)
allbands[w_n].setOrigin(
chip.bound.xmin, chip.bound.ymin)
allbands[w_n].wcs = wcs_fp
orders[o_n] = allbands
chip.img_stack[gn] = orders
else:
chip.img_stack = {}
# Get random generators for this chip # Get random generators for this chip
chip.rng_poisson, chip.poisson_noise = chip_utils.get_poisson( chip.rng_poisson, chip.poisson_noise = chip_utils.get_poisson(
seed=int(self.config["random_seeds"]["seed_poisson"]) + pointing.id*30 + chip.chipID, sky_level=0.) seed=int(self.config["random_seeds"]["seed_poisson"]) + pointing.id*30 + chip.chipID, sky_level=0.)
...@@ -66,7 +45,7 @@ class Observation(object): ...@@ -66,7 +45,7 @@ class Observation(object):
chip.shutter_img = np.ones_like(chip.img.array) chip.shutter_img = np.ones_like(chip.img.array)
else: else:
chip.shutter_img = effects.ShutterEffectArr( chip.shutter_img = effects.ShutterEffectArr(
chip.img, t_exp=pointing.exp_time, t_shutter=1.3, dist_bearing=735, dt=1E-3) chip.img, t_shutter=1.3, dist_bearing=735, dt=1E-3)
chip.prnu_img = effects.PRNU_Img(xsize=chip.npix_x, ysize=chip.npix_y, sigma=0.01, chip.prnu_img = effects.PRNU_Img(xsize=chip.npix_x, ysize=chip.npix_y, sigma=0.01,
seed=int(self.config["random_seeds"]["seed_prnu"]+chip.chipID)) seed=int(self.config["random_seeds"]["seed_prnu"]+chip.chipID))
...@@ -112,24 +91,17 @@ class Observation(object): ...@@ -112,24 +91,17 @@ class Observation(object):
input_date_str=date_str, input_date_str=date_str,
input_time_str=time_str input_time_str=time_str
) )
ra_offset, dec_offset = pointing.ra - \ ra_cen, dec_cen = ra_cen[0], dec_cen[0]
ra_cen[0], pointing.dec - dec_cen[0]
else: else:
ra_offset, dec_offset = 0., 0.
ra_cen = pointing.ra ra_cen = pointing.ra
dec_cen = pointing.dec dec_cen = pointing.dec
slsPSFOpt = False
# Prepare necessary chip properties for simulation # Prepare necessary chip properties for simulation
chip = self.prepare_chip_for_exposure( chip = self.prepare_chip_for_exposure(chip, ra_cen, dec_cen, pointing)
chip, ra_cen, dec_cen, pointing, slsPSFOptim=slsPSFOpt)
# Initialize SimSteps # Initialize SimSteps
sim_steps = SimSteps(overall_config=self.config, sim_steps = SimSteps(overall_config=self.config,
chip_output=chip_output, chip_output=chip_output, all_filters=self.all_filters)
all_filters=self.all_filters,
ra_offset=ra_offset,
dec_offset=dec_offset)
for step in pointing.obs_param["call_sequence"]: for step in pointing.obs_param["call_sequence"]:
if self.config["run_option"]["out_cat_only"]: if self.config["run_option"]["out_cat_only"]:
...@@ -158,9 +130,6 @@ class Observation(object): ...@@ -158,9 +130,6 @@ class Observation(object):
chip_output.Log_info("check running:1: pointing-%d chip-%d pid-%d memory-%6.2fGB" % (pointing.id, chip_output.Log_info("check running:1: pointing-%d chip-%d pid-%d memory-%6.2fGB" % (pointing.id,
chip.chipID, os.getpid(), (psutil.Process(os.getpid()).memory_info().rss / 1024 / 1024 / 1024))) chip.chipID, os.getpid(), (psutil.Process(os.getpid()).memory_info().rss / 1024 / 1024 / 1024)))
del chip.img del chip.img
del chip.flat_img
del chip.prnu_img
del chip.shutter_img
def runExposure_MPI_PointingList(self, pointing_list, chips=None): def runExposure_MPI_PointingList(self, pointing_list, chips=None):
comm = MPI.COMM_WORLD comm = MPI.COMM_WORLD
......
observation_sim/config/ChipOutput.py
View file @ 64eb7956
...@@ -78,12 +78,12 @@ class ChipOutput(object): ...@@ -78,12 +78,12 @@ class ChipOutput(object):
self.hdr += "\n" self.hdr += "\n"
self.cat.write(self.hdr) self.cat.write(self.hdr)
def cat_add_obj(self, obj, pos_img, pos_shear, ra_offset=0., dec_offset=0.): def cat_add_obj(self, obj, pos_img, pos_shear):
ximg = obj.real_pos.x + 1.0 ximg = obj.real_pos.x + 1.0
yimg = obj.real_pos.y + 1.0 yimg = obj.real_pos.y + 1.0
line = self.fmt % ( line = self.fmt % (
obj.id, int(self.chip_label), self.filt.filter_type, ximg, yimg, obj.ra + ra_offset, obj.dec + dec_offset, obj.ra_orig, obj.dec_orig, obj.z, obj.getMagFilter( obj.id, int(self.chip_label), self.filt.filter_type, ximg, yimg, obj.ra, obj.dec, obj.ra_orig, obj.dec_orig, obj.z, obj.getMagFilter(
self.filt), obj.type, self.filt), obj.type,
obj.pmra, obj.pmdec, obj.rv, obj.parallax) obj.pmra, obj.pmdec, obj.rv, obj.parallax)
line += obj.additional_output_str line += obj.additional_output_str
......
observation_sim/config/header/ImageHeader.py
View file @ 64eb7956
...@@ -434,7 +434,7 @@ def generatePrimaryHeader(xlen=9216, ylen=9232, pointing_id='00000001', pointing ...@@ -434,7 +434,7 @@ def generatePrimaryHeader(xlen=9216, ylen=9232, pointing_id='00000001', pointing
co = coord.SkyCoord(ra, dec, unit='deg') co = coord.SkyCoord(ra, dec, unit='deg')
ra_hms = format(co.ra.hms.h, '02.0f') + format(co.ra.hms.m, ra_hms = format(co.ra.hms.h, '02.0f') + format(co.ra.hms.m,
'02.0f') + format(co.ra.hms.s, '04.1f') '02.0f') + format(co.ra.hms.s, '02.1f')
dec_hms = format(co.dec.dms.d, '02.0f') + format(abs(co.dec.dms.m), dec_hms = format(co.dec.dms.d, '02.0f') + format(abs(co.dec.dms.m),
'02.0f') + format(abs(co.dec.dms.s), '02.0f') '02.0f') + format(abs(co.dec.dms.s), '02.0f')
if dec >= 0: if dec >= 0:
......
observation_sim/instruments/chip/effects.py
View file @ 64eb7956
...@@ -793,7 +793,7 @@ def ShutterEffectArr(GSImage, t_exp=150, t_shutter=1.3, dist_bearing=735, dt=1E- ...@@ -793,7 +793,7 @@ def ShutterEffectArr(GSImage, t_exp=150, t_shutter=1.3, dist_bearing=735, dt=1E-
sizey = ymax-ymin+1 sizey = ymax-ymin+1
xnewgrid = np.mgrid[xmin:(xmin+sizex)] xnewgrid = np.mgrid[xmin:(xmin+sizex)]
expeffect = interpolate.splev(xnewgrid, intp, der=0) expeffect = interpolate.splev(xnewgrid, intp, der=0)
expeffect /= t_exp expeffect /= np.max(expeffect)
exparrnormal = np.tile(expeffect, (sizey,1)) exparrnormal = np.tile(expeffect, (sizey,1))
# GSImage *= exparrnormal # GSImage *= exparrnormal
......
observation_sim/instruments/data/filters/g_sub.list
View file @ 64eb7956
2000 3800
4217 4217
4432 4432
4631 4631
...@@ -6,4 +6,4 @@ ...@@ -6,4 +6,4 @@
5002 5002
5179 5179
5354 5354
11000 5799
\ No newline at end of file \ No newline at end of file
observation_sim/instruments/data/filters/i_sub.list
View file @ 64eb7956
2000 6600
7061 7061
7255 7255
7448 7448
...@@ -6,4 +6,4 @@ ...@@ -6,4 +6,4 @@
7833 7833
8027 8027
8226 8226
11000 8999
\ No newline at end of file \ No newline at end of file
observation_sim/instruments/data/filters/nuv_sub.list
View file @ 64eb7956
2000 2513
2621 2621
2716 2716
2805 2805
...@@ -6,4 +6,4 @@ ...@@ -6,4 +6,4 @@
2969 2969
3050 3050
3132 3132
11000 3499
\ No newline at end of file \ No newline at end of file
observation_sim/instruments/data/filters/r_sub.list
View file @ 64eb7956
2000 5100
5642 5642
5821 5821
6001 6001
...@@ -6,4 +6,4 @@ ...@@ -6,4 +6,4 @@
6363 6363
6547 6547
6735 6735
11000 7199
\ No newline at end of file \ No newline at end of file
observation_sim/instruments/data/filters/u_sub.list
View file @ 64eb7956
2000 3000
3277 3277
3380 3380
3485 3485
...@@ -6,4 +6,4 @@ ...@@ -6,4 +6,4 @@
3703 3703
3813 3813
3918 3918
11000 4499
\ No newline at end of file \ No newline at end of file
observation_sim/instruments/data/filters/y_sub.list
View file @ 64eb7956
2000 9000
9322 9322
9405 9405
9489 9489
...@@ -6,4 +6,4 @@ ...@@ -6,4 +6,4 @@
9695 9695
9832 9832
10024 10024
11000 10590
\ No newline at end of file \ No newline at end of file
observation_sim/instruments/data/filters/z_sub.list
View file @ 64eb7956
2000 7800
8494 8494
8638 8638
8790 8790
...@@ -6,4 +6,4 @@ ...@@ -6,4 +6,4 @@
9141 9141
9363 9363
9663 9663
11000 10590
\ No newline at end of file \ No newline at end of file
observation_sim/mock_objects/ExtinctionMW.py deleted 100644 → 0
View file @ 914004c1
import os
import numpy as np
import healpy as hp
from astropy.io import fits
from astropy import units as u
from astropy.coordinates import SkyCoord
class ExtinctionMW(object):
def __init__(self):
self.Rv = 3.1
self.lamb = np.arange(2000, 11001+0.5, 0.5)
@staticmethod
def radec2pix(ra, dec, NSIDE=2048):
return hp.pixelfunc.ang2pix(nside=NSIDE, theta=ra, phi=dec, nest=True, lonlat=True)
def init_ext_model(self, model_name="odonnell", Av=1.0, Rv=3.1, lamb=None):
self.model_name = model_name
self.Av = Av
self.Rv = Rv
if lamb is not None:
self.lamb = lamb
if self.model_name == "odonnell":
alpha = np.zeros(self.lamb.shape)
beta = np.zeros(self.lamb.shape)
x = 1.e4 / self.lamb
def f_alpha_1(x): return 0.574 * (x**1.61)
def f_beta_1(x): return -0.527 * (x**1.61)
def f_alpha_2(x):
y = x - 1.82
return 1 + 0.104*y - 0.609*(y**2) + 0.701*(y**3) + \
1.137*(y**4) - 1.718*(y**5) - 0.827 * \
(y**6) + 1.647*(y**7) - 0.505*(y**8)
def f_beta_2(x):
y = x - 1.82
return 1.952*y + 2.908*(y**2) - 3.989*(y**3) - 7.985 * \
(y**4) + 11.102*(y**5) + 5.491 * \
(y**6) - 10.805*(y**7) + 3.347*(y**8)
def f_alpha_3(x): return 1.752 - 0.316*x - \
0.104 / ((x - 4.67)**2 + 0.341)
def f_beta_3(x): return -3.090 + 1.825*x + \
1.206 / ((x - 4.62)**2 + 0.262)
def f_alpha_4(x): return f_alpha_3(x) + -0.04473 * \
(x - 5.9)**2 - 0.009779 * (x - 5.9)**3
def f_beta_4(x): return f_beta_3(x) + 0.2130 * \
(x - 5.9)**2 + 0.1207 * (x - 5.9)**3
def f_alpha_5(x): return -1.073 - 0.628*(x - 8) + \
0.137*(x - 8)**2 - 0.070*(x - 8)**3
def f_beta_5(x): return 13.670 + 4.257*(x - 8) - \
0.420*(x - 8)**2 + 0.374 * (x - 8)**3
alpha = np.piecewise(
x, [x <= 1.1, (x > 1.1)*(x <= 3.3), (x > 3.3)*(x <= 5.9), (x > 5.9)*(x <= 8.), (x > 8.)], [f_alpha_1, f_alpha_2, f_alpha_3, f_alpha_4, f_alpha_5])
beta = np.piecewise(
x, [x <= 1.1, (x > 1.1)*(x <= 3.3), (x > 3.3)*(x <= 5.9), (x > 5.9)*(x <= 8.), (x > 8.)], [f_beta_1, f_beta_2, f_beta_3, f_beta_4, f_beta_5])
self.ext = (alpha + beta / self.Rv) * self.Av
def load_Planck_ext(self, file_path):
hdu = fits.open(file_path)
self.ebv_planck = hdu[1].data["EBV"]
def Av_from_Planck(self, ra, dec):
if not hasattr(self, 'ebv_planck'):
raise ValueError(
"Need to load planck dust map first")
# Convert to galactic coordinates
c = SkyCoord(ra=ra * u.degree, dec=dec *
u.degree, frame='icrs').galactic
l, b = c.l.radian, c.b.radian
NSIDE = hp.pixelfunc.get_nside(self.ebv_planck)
pix = hp.pixelfunc.ang2pix(
nside=NSIDE, theta=np.pi/2. - b, phi=l, nest=True)
return self.ebv_planck[pix] * self.Rv
def apply_extinction(self, spec, Av=1.0):
if len(spec) != len(self.lamb):
raise ValueError(
"Dimension of spec do not match that of ExtinctionMW.lamb: try initilize (init_ext_model) by the actual size of spec which you want to apply the extinction to")
if not hasattr(self, 'ext'):
raise ValueError(
"Need to initialize the extinction model (init_ext_model) first")
scale = 10**(-.4 * self.ext * Av)
# print("scale = ", scale)
spec *= scale
return spec
observation_sim/mock_objects/Galaxy.py
View file @ 64eb7956
...@@ -323,73 +323,28 @@ class Galaxy(MockObject): ...@@ -323,73 +323,28 @@ class Galaxy(MockObject):
# # if fd_shear is not None: # # if fd_shear is not None:
# # gal = gal.shear(fd_shear) # # gal = gal.shear(fd_shear)
galImg_List = [] starImg = gal.drawImage(
try: wcs=chip_wcs_local, offset=offset, method='real_space')
pos_img_local = [0,0]
x_start = chip.x_cen/chip.pix_size - chip.npix_x / 2.
y_start = chip.y_cen/chip.pix_size - chip.npix_y / 2.
pos_img_local[0] = pos_img.x - x_start
pos_img_local[1] = pos_img.y - y_start
nnx = 0
nny = 0
for order in ["A","B"]:
psf, pos_shear = psf_model.get_PSF(
chip, pos_img_local=pos_img_local, bandNo=i+1, galsimGSObject=True, g_order=order, grating_split_pos=grating_split_pos)
star_p = galsim.Convolve(psf, gal)
if nnx == 0:
galImg = star_p.drawImage(wcs=chip_wcs_local, offset=offset)
nnx = galImg.xmax - galImg.xmin + 1
nny = galImg.ymax - galImg.ymin + 1
else:
galImg = star_p.drawImage(nx = nnx, ny = nny, wcs=chip_wcs_local, offset=offset)
galImg.setOrigin(0, 0)
# n1 = np.sum(np.isinf(galImg.array))
# n2 = np.sum(np.isnan(galImg.array))
# if n1>0 or n2 > 0:
# print("DEBUG: Galaxy, inf:%d, nan:%d"%(n1, n2))
if np.sum(np.isnan(galImg.array)) > 0:
# ERROR happens
return 2, pos_shear
galImg_List.append(galImg)
for order in ["C","D","E"]:
galImg_List.append(galImg)
except:
psf, pos_shear = psf_model.get_PSF(chip=chip, pos_img=pos_img)
star_p = galsim.Convolve(psf, gal)
galImg = star_p.drawImage(wcs=chip_wcs_local, offset=offset)
galImg.setOrigin(0, 0)
if np.sum(np.isnan(galImg.array)) > 0:
# ERROR happens
return 2, pos_shear
for order in ["A","B","C","D","E"]:
galImg_List.append(galImg)
# starImg = gal.drawImage(
# wcs=chip_wcs_local, offset=offset, method='real_space')
origin_star = [y_nominal - (galImg.center.y - galImg.ymin), origin_star = [y_nominal - (starImg.center.y - starImg.ymin),
x_nominal - (galImg.center.x - galImg.xmin)] x_nominal - (starImg.center.x - starImg.xmin)]
galImg.setOrigin(0, 0) starImg.setOrigin(0, 0)
gal_origin = [origin_star[0], origin_star[1]] gal_origin = [origin_star[0], origin_star[1]]
gal_end = [origin_star[0] + galImg.array.shape[0] - gal_end = [origin_star[0] + starImg.array.shape[0] -
1, origin_star[1] + galImg.array.shape[1] - 1] 1, origin_star[1] + starImg.array.shape[1] - 1]
if gal_origin[1] < grating_split_pos_chip < gal_end[1]: if gal_origin[1] < grating_split_pos_chip < gal_end[1]:
subSlitPos = int(grating_split_pos_chip - gal_origin[1] + 1) subSlitPos = int(grating_split_pos_chip - gal_origin[1] + 1)
# part img disperse # part img disperse
star_p1s=[] subImg_p1 = starImg.array[:, 0:subSlitPos]
for galImg in galImg_List:
subImg_p1 = galImg.array[:, 0:subSlitPos]
star_p1 = galsim.Image(subImg_p1) star_p1 = galsim.Image(subImg_p1)
star_p1.setOrigin(0, 0) star_p1.setOrigin(0, 0)
star_p1s.append(star_p1)
origin_p1 = origin_star origin_p1 = origin_star
xcenter_p1 = min(x_nominal, grating_split_pos_chip-1) - 0 xcenter_p1 = min(x_nominal, grating_split_pos_chip-1) - 0
ycenter_p1 = y_nominal-0 ycenter_p1 = y_nominal-0
sdp_p1 = SpecDisperser(orig_img=star_p1s, xcenter=xcenter_p1, sdp_p1 = SpecDisperser(orig_img=star_p1, xcenter=xcenter_p1,
ycenter=ycenter_p1, origin=origin_p1, ycenter=ycenter_p1, origin=origin_p1,
tar_spec=normalSED, tar_spec=normalSED,
band_start=brange[0], band_end=brange[1], band_start=brange[0], band_end=brange[1],
...@@ -397,25 +352,21 @@ class Galaxy(MockObject): ...@@ -397,25 +352,21 @@ class Galaxy(MockObject):
isAlongY=0, isAlongY=0,
flat_cube=flat_cube) flat_cube=flat_cube)
self.addSLStoChipImage(sdp=sdp_p1, chip=chip, xOrderSigPlus = xOrderSigPlus, local_wcs=chip_wcs_local) # 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], pos_shear = self.addSLStoChipImageWithPSF(sdp=sdp_p1, chip=chip, pos_img_local=[xcenter_p1, ycenter_p1],
# psf_model=psf_model, bandNo=i + 1, psf_model=psf_model, bandNo=i + 1,
# grating_split_pos=grating_split_pos, grating_split_pos=grating_split_pos,
# local_wcs=chip_wcs_local, pos_img=pos_img) local_wcs=chip_wcs_local, pos_img=pos_img)
star_p2s=[]
for galImg in galImg_List:
subImg_p2 = galImg.array[:, subImg_p2 = starImg.array[:,
subSlitPos + 1:galImg.array.shape[1]] subSlitPos+1:starImg.array.shape[1]]
star_p2 = galsim.Image(subImg_p2) star_p2 = galsim.Image(subImg_p2)
star_p2.setOrigin(0, 0) star_p2.setOrigin(0, 0)
star_p2s.append(star_p2)
origin_p2 = [origin_star[0], grating_split_pos_chip] origin_p2 = [origin_star[0], grating_split_pos_chip]
xcenter_p2 = max(x_nominal, grating_split_pos_chip - 1) - 0 xcenter_p2 = max(x_nominal, grating_split_pos_chip - 1) - 0
ycenter_p2 = y_nominal - 0 ycenter_p2 = y_nominal - 0
sdp_p2 = SpecDisperser(orig_img=star_p2s, xcenter=xcenter_p2, sdp_p2 = SpecDisperser(orig_img=star_p2, xcenter=xcenter_p2,
ycenter=ycenter_p2, origin=origin_p2, ycenter=ycenter_p2, origin=origin_p2,
tar_spec=normalSED, tar_spec=normalSED,
band_start=brange[0], band_end=brange[1], band_start=brange[0], band_end=brange[1],
...@@ -423,41 +374,41 @@ class Galaxy(MockObject): ...@@ -423,41 +374,41 @@ class Galaxy(MockObject):
isAlongY=0, isAlongY=0,
flat_cube=flat_cube) flat_cube=flat_cube)
self.addSLStoChipImage(sdp=sdp_p2, chip=chip, xOrderSigPlus = xOrderSigPlus, local_wcs=chip_wcs_local) # 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], pos_shear = self.addSLStoChipImageWithPSF(sdp=sdp_p2, chip=chip, pos_img_local=[xcenter_p2, ycenter_p2],
# psf_model=psf_model, bandNo=i + 1, psf_model=psf_model, bandNo=i + 1,
# grating_split_pos=grating_split_pos, grating_split_pos=grating_split_pos,
# local_wcs=chip_wcs_local, pos_img=pos_img) local_wcs=chip_wcs_local, pos_img=pos_img)
del sdp_p1 del sdp_p1
del sdp_p2 del sdp_p2
elif grating_split_pos_chip <= gal_origin[1]: elif grating_split_pos_chip <= gal_origin[1]:
sdp = SpecDisperser(orig_img=galImg_List, xcenter=x_nominal - 0, sdp = SpecDisperser(orig_img=starImg, xcenter=x_nominal - 0,
ycenter=y_nominal - 0, origin=origin_star, ycenter=y_nominal - 0, origin=origin_star,
tar_spec=normalSED, tar_spec=normalSED,
band_start=brange[0], band_end=brange[1], band_start=brange[0], band_end=brange[1],
conf=chip.sls_conf[1], conf=chip.sls_conf[1],
isAlongY=0, isAlongY=0,
flat_cube=flat_cube) flat_cube=flat_cube)
self.addSLStoChipImage(sdp=sdp, chip=chip, xOrderSigPlus = xOrderSigPlus, local_wcs=chip_wcs_local) # 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], pos_shear = self.addSLStoChipImageWithPSF(sdp=sdp, chip=chip, pos_img_local=[x_nominal, y_nominal],
# psf_model=psf_model, bandNo=i + 1, psf_model=psf_model, bandNo=i + 1,
# grating_split_pos=grating_split_pos, grating_split_pos=grating_split_pos,
# local_wcs=chip_wcs_local, pos_img=pos_img) local_wcs=chip_wcs_local, pos_img=pos_img)
del sdp del sdp
elif grating_split_pos_chip >= gal_end[1]: elif grating_split_pos_chip >= gal_end[1]:
sdp = SpecDisperser(orig_img=galImg_List, xcenter=x_nominal - 0, sdp = SpecDisperser(orig_img=starImg, xcenter=x_nominal - 0,
ycenter=y_nominal - 0, origin=origin_star, ycenter=y_nominal - 0, origin=origin_star,
tar_spec=normalSED, tar_spec=normalSED,
band_start=brange[0], band_end=brange[1], band_start=brange[0], band_end=brange[1],
conf=chip.sls_conf[0], conf=chip.sls_conf[0],
isAlongY=0, isAlongY=0,
flat_cube=flat_cube) flat_cube=flat_cube)
self.addSLStoChipImage(sdp=sdp, chip=chip, xOrderSigPlus = xOrderSigPlus, local_wcs=chip_wcs_local) # 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], pos_shear = self.addSLStoChipImageWithPSF(sdp=sdp, chip=chip, pos_img_local=[x_nominal, y_nominal],
# psf_model=psf_model, bandNo=i + 1, psf_model=psf_model, bandNo=i + 1,
# grating_split_pos=grating_split_pos, grating_split_pos=grating_split_pos,
# local_wcs=chip_wcs_local, pos_img=pos_img) local_wcs=chip_wcs_local, pos_img=pos_img)
del sdp del sdp
# print(self.y_nominal, starImg.center.y, starImg.ymin) # print(self.y_nominal, starImg.center.y, starImg.ymin)
......
observation_sim/mock_objects/MockObject.py
View file @ 64eb7956
...@@ -11,8 +11,6 @@ from observation_sim.mock_objects._util import integrate_sed_bandpass, getNormFa ...@@ -11,8 +11,6 @@ from observation_sim.mock_objects._util import integrate_sed_bandpass, getNormFa
getABMAG getABMAG
from observation_sim.mock_objects.SpecDisperser import SpecDisperser from observation_sim.mock_objects.SpecDisperser import SpecDisperser
from observation_sim.instruments.chip import chip_utils
class MockObject(object): class MockObject(object):
def __init__(self, param, logger=None): def __init__(self, param, logger=None):
...@@ -61,13 +59,13 @@ class MockObject(object): ...@@ -61,13 +59,13 @@ class MockObject(object):
flux = self.getFluxFilter(filt) flux = self.getFluxFilter(filt)
return flux * tel.pupil_area * exptime return flux * tel.pupil_area * exptime
def getPosWorld(self, ra_offset=0., dec_offset=0.): def getPosWorld(self):
ra = self.param["ra"] + ra_offset ra = self.param["ra"]
dec = self.param["dec"] + dec_offset dec = self.param["dec"]
return galsim.CelestialCoord(ra=ra * galsim.degrees, dec=dec * galsim.degrees) 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, ra_offset=0., dec_offset=0.): 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(ra_offset, dec_offset)) self.posImg = img.wcs.toImage(self.getPosWorld())
self.localWCS = img.wcs.local(self.posImg) self.localWCS = img.wcs.local(self.posImg)
# Apply field distortion model # Apply field distortion model
if (fdmodel is not None) and (chip is not None): if (fdmodel is not None) and (chip is not None):
...@@ -120,6 +118,7 @@ class MockObject(object): ...@@ -120,6 +118,7 @@ class MockObject(object):
if self.logger: if self.logger:
self.logger.error(e) self.logger.error(e)
return 2, None return 2, None
# Set Galsim Parameters # Set Galsim Parameters
if self.getMagFilter(filt) <= 15: if self.getMagFilter(filt) <= 15:
folding_threshold = 5.e-4 folding_threshold = 5.e-4
...@@ -196,26 +195,33 @@ class MockObject(object): ...@@ -196,26 +195,33 @@ class MockObject(object):
# DEBUG # DEBUG
######################################################### #########################################################
# print("before convolveGaussXorders, img_s:", img_s) # print("before convolveGaussXorders, img_s:", img_s)
# nan_ids = np.isnan(img_s) nan_ids = np.isnan(img_s)
# if img_s[nan_ids].shape[0] > 0: if img_s[nan_ids].shape[0] > 0:
# # img_s[nan_ids] = 0 # img_s[nan_ids] = 0
# print("DEBUG: before convolveGaussXorders specImg nan num is", print("DEBUG: before convolveGaussXorders specImg nan num is",
# img_s[nan_ids].shape[0]) img_s[nan_ids].shape[0])
######################################################### #########################################################
# img_s, orig_off = convolveGaussXorders(img_s, xOrderSigPlus[k]) img_s, orig_off = convolveGaussXorders(img_s, xOrderSigPlus[k])
orig_off = 0
origin_order_x = v[1] - orig_off origin_order_x = v[1] - orig_off
origin_order_y = v[2] - orig_off origin_order_y = v[2] - orig_off
######################################################### #########################################################
# DEBUG # DEBUG
######################################################### #########################################################
# print("DEBUG: orig_off is", orig_off) # print("DEBUG: orig_off is", orig_off)
# nan_ids = np.isnan(img_s) nan_ids = np.isnan(img_s)
# if img_s[nan_ids].shape[0] > 0: if img_s[nan_ids].shape[0] > 0:
# img_s[nan_ids] = 0 img_s[nan_ids] = 0
# print("DEBUG: specImg nan num is", img_s[nan_ids].shape[0]) print("DEBUG: specImg nan num is", img_s[nan_ids].shape[0])
######################################################### #########################################################
stamp = galsim.ImageF(img_s) 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.wcs = local_wcs
stamp.setOrigin(origin_order_x, origin_order_y) stamp.setOrigin(origin_order_x, origin_order_y)
...@@ -224,125 +230,57 @@ class MockObject(object): ...@@ -224,125 +230,57 @@ class MockObject(object):
if bounds.area() == 0: if bounds.area() == 0:
continue continue
chip.img.setOrigin(0, 0) chip.img.setOrigin(0, 0)
chip.img[bounds] = chip.img[bounds]+stamp[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) chip.img.setOrigin(chip.bound.xmin, chip.bound.ymin)
del stamp del stamp
del spec_orders 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): 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() spec_orders = sdp.compute_spec_orders()
pos_shear = galsim.Shear(e=0., beta=(np.pi/2)*galsim.radians)
if chip.slsPSFOptim:
for k, v in spec_orders.items(): for k, v in spec_orders.items():
img_s = v[0] 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)
nan_ids = np.isnan(img_s) psf_img = psf.drawImage(nx=100, ny=100, wcs=local_wcs)
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)
orig_off = [0,0]
origin_order_x = v[1] - orig_off[0]
origin_order_y = v[2] - orig_off[1]
specImg = galsim.ImageF(img_s)
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
# orders = {'A': 'order1', 'B': 'order0', 'C': 'order2', 'D': 'order-1', 'E': 'order-2'}
orders = {'A': 'order1', 'B': 'order0', 'C': 'order0', 'D': 'order0', 'E': 'order0'}
gratingN = chip_utils.getChipSLSGratingID(chip.chipID)[1]
if pos_img_local[0] < grating_split_pos:
gratingN = chip_utils.getChipSLSGratingID(chip.chipID)[0]
chip.img_stack[gratingN][orders[k]]['w' + str(bandNo)].setOrigin(0, 0)
chip.img_stack[gratingN][orders[k]]['w' + str(bandNo)][bounds] = chip.img_stack[gratingN][orders[k]]['w' + str(bandNo)][bounds] + specImg[bounds]
chip.img_stack[gratingN][orders[k]]['w' + str(bandNo)].setOrigin(chip.bound.xmin, chip.bound.ymin)
else:
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 psf_img_m = psf_img.array
# specImg.setOrigin(origin_order_x, origin_order_y)
# print('DEBUG: BEGIN -----------',bandNo,k) #########################################################
# DEBUG
#########################################################
# ids_p = psf_img_m < 0
# psf_img_m[ids_p] = 0
img_s = v[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) nan_ids = np.isnan(img_s)
if img_s[nan_ids].shape[0] > 0: if img_s[nan_ids].shape[0] > 0:
img_s[nan_ids] = 0 img_s[nan_ids] = 0
print("DEBUG: specImg nan num is", img_s[nan_ids].shape[0]) print("DEBUG: specImg nan num is", img_s[nan_ids].shape[0])
######################################################### #########################################################
origin_order_x = v[1]
origin_order_y = v[2]
specImg = galsim.ImageF(img_s)
specImg.wcs = local_wcs
specImg.setOrigin(origin_order_x, origin_order_y)
try:
specImg = psf_model.get_PSF_AND_convolve_withsubImg(chip, cutImg=specImg, pos_img_local=pos_img_local, bandNo=bandNo, 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
img_s, orig_off = convolveImg(img_s, psf_img_m) img_s, orig_off = convolveImg(img_s, psf_img_m)
origin_order_x = v[1] - orig_off[0] origin_order_x = v[1] - orig_off[0]
origin_order_y = v[2] - orig_off[1] origin_order_y = v[2] - orig_off[1]
specImg = galsim.ImageF(img_s) 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.wcs = local_wcs
specImg.setOrigin(origin_order_x, origin_order_y) specImg.setOrigin(origin_order_x, origin_order_y)
...@@ -421,76 +359,32 @@ class MockObject(object): ...@@ -421,76 +359,32 @@ class MockObject(object):
# psf, pos_shear = psf_model.get_PSF(chip=chip, pos_img=pos_img, bandpass=bandpass, # psf, pos_shear = psf_model.get_PSF(chip=chip, pos_img=pos_img, bandpass=bandpass,
# folding_threshold=folding_threshold) # 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)
# star = galsim.DeltaFunction(gsparams=gsp) starImg = star.drawImage(
# star = star.withFlux(tel.pupil_area * exptime) nx=60, ny=60, wcs=chip_wcs_local, offset=offset)
#psf list :["A","B","C","D","E"]
starImg_List = []
try:
pos_img_local = [0,0]
x_start = chip.x_cen/chip.pix_size - chip.npix_x / 2.
y_start = chip.y_cen/chip.pix_size - chip.npix_y / 2.
pos_img_local[0] = pos_img.x - x_start
pos_img_local[1] = pos_img.y - y_start
nnx = 0
nny = 0
for order in ["A","B"]:
psf, pos_shear = psf_model.get_PSF(
chip, pos_img_local=pos_img_local, bandNo=i+1, galsimGSObject=True, g_order=order, grating_split_pos=grating_split_pos)
# star_p = galsim.Convolve(psf, star)
star_p = psf.withFlux(tel.pupil_area * exptime)
if nnx == 0:
starImg = star_p.drawImage(wcs=chip_wcs_local, offset=offset)
nnx = starImg.xmax - starImg.xmin + 1
nny = starImg.ymax - starImg.ymin + 1
else:
starImg = star_p.drawImage(nx = nnx, ny = nny, wcs=chip_wcs_local, offset=offset)
# n1 = np.sum(np.isinf(starImg.array))
# n2 = np.sum(np.isnan(starImg.array))
# if n1>0 or n2 > 0:
# print("DEBUG: MockObject, inf:%d, nan:%d"%(n1, n2))
starImg.setOrigin(0, 0)
starImg_List.append(starImg)
for order in ["C","D","E"]:
starImg_List.append(starImg)
except:
psf, pos_shear = psf_model.get_PSF(chip=chip, pos_img=pos_img)
# star_p = galsim.Convolve(psf, star)
star_p = psf.withFlux(tel.pupil_area * exptime)
starImg = star_p.drawImage(wcs=chip_wcs_local, offset=offset)
starImg.setOrigin(0, 0)
for order in ["A","B","C","D","E"]:
starImg_List.append(starImg)
# 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), origin_star = [y_nominal - (starImg.center.y - starImg.ymin),
x_nominal - (starImg.center.x - starImg.xmin)] x_nominal - (starImg.center.x - starImg.xmin)]
starImg.setOrigin(0, 0)
gal_origin = [origin_star[0], origin_star[1]] gal_origin = [origin_star[0], origin_star[1]]
gal_end = [origin_star[0] + starImg.array.shape[0] - gal_end = [origin_star[0] + starImg.array.shape[0] -
1, origin_star[1] + starImg.array.shape[1] - 1] 1, origin_star[1] + starImg.array.shape[1] - 1]
if gal_origin[1] < grating_split_pos_chip < gal_end[1]: if gal_origin[1] < grating_split_pos_chip < gal_end[1]:
subSlitPos = int(grating_split_pos_chip - gal_origin[1] + 1) subSlitPos = int(grating_split_pos_chip - gal_origin[1] + 1)
# part img disperse # part img disperse
star_p1s=[]
for starImg in starImg_List:
subImg_p1 = starImg.array[:, 0:subSlitPos] subImg_p1 = starImg.array[:, 0:subSlitPos]
star_p1 = galsim.Image(subImg_p1) star_p1 = galsim.Image(subImg_p1)
star_p1.setOrigin(0, 0)
star_p1s.append(star_p1)
origin_p1 = origin_star origin_p1 = origin_star
star_p1.setOrigin(0, 0)
xcenter_p1 = min(x_nominal, grating_split_pos_chip - 1) - 0 xcenter_p1 = min(x_nominal, grating_split_pos_chip - 1) - 0
ycenter_p1 = y_nominal - 0 ycenter_p1 = y_nominal - 0
sdp_p1 = SpecDisperser(orig_img=star_p1s, xcenter=xcenter_p1, sdp_p1 = SpecDisperser(orig_img=star_p1, xcenter=xcenter_p1,
ycenter=ycenter_p1, origin=origin_p1, ycenter=ycenter_p1, origin=origin_p1,
tar_spec=normalSED, tar_spec=normalSED,
band_start=brange[0], band_end=brange[1], band_start=brange[0], band_end=brange[1],
...@@ -498,25 +392,20 @@ class MockObject(object): ...@@ -498,25 +392,20 @@ class MockObject(object):
isAlongY=0, isAlongY=0,
flat_cube=flat_cube) flat_cube=flat_cube)
self.addSLStoChipImage(sdp=sdp_p1, chip=chip, xOrderSigPlus=xOrderSigPlus, local_wcs=chip_wcs_local) # 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], 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, psf_model=psf_model, bandNo=i+1, grating_split_pos=grating_split_pos,
# local_wcs=chip_wcs_local, pos_img=pos_img) local_wcs=chip_wcs_local, pos_img=pos_img)
star_p2s=[]
for starImg in starImg_List:
subImg_p2 = starImg.array[:, subImg_p2 = starImg.array[:,
subSlitPos + 1:starImg.array.shape[1]] subSlitPos + 1:starImg.array.shape[1]]
star_p2 = galsim.Image(subImg_p2) star_p2 = galsim.Image(subImg_p2)
star_p2.setOrigin(0, 0) star_p2.setOrigin(0, 0)
star_p2s.append(star_p2)
origin_p2 = [origin_star[0], grating_split_pos_chip] origin_p2 = [origin_star[0], grating_split_pos_chip]
xcenter_p2 = max(x_nominal, grating_split_pos_chip - 1) - 0 xcenter_p2 = max(x_nominal, grating_split_pos_chip - 1) - 0
ycenter_p2 = y_nominal - 0 ycenter_p2 = y_nominal - 0
sdp_p2 = SpecDisperser(orig_img=star_p2s, xcenter=xcenter_p2, sdp_p2 = SpecDisperser(orig_img=star_p2, xcenter=xcenter_p2,
ycenter=ycenter_p2, origin=origin_p2, ycenter=ycenter_p2, origin=origin_p2,
tar_spec=normalSED, tar_spec=normalSED,
band_start=brange[0], band_end=brange[1], band_start=brange[0], band_end=brange[1],
...@@ -524,38 +413,38 @@ class MockObject(object): ...@@ -524,38 +413,38 @@ class MockObject(object):
isAlongY=0, isAlongY=0,
flat_cube=flat_cube) flat_cube=flat_cube)
self.addSLStoChipImage(sdp=sdp_p2, chip=chip, xOrderSigPlus=xOrderSigPlus, local_wcs=chip_wcs_local) # 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], 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, psf_model=psf_model, bandNo=i + 1, grating_split_pos=grating_split_pos,
# local_wcs=chip_wcs_local, pos_img=pos_img) local_wcs=chip_wcs_local, pos_img=pos_img)
del sdp_p1 del sdp_p1
del sdp_p2 del sdp_p2
elif grating_split_pos_chip <= gal_origin[1]: elif grating_split_pos_chip <= gal_origin[1]:
sdp = SpecDisperser(orig_img=starImg_List, xcenter=x_nominal - 0, sdp = SpecDisperser(orig_img=starImg, xcenter=x_nominal - 0,
ycenter=y_nominal - 0, origin=origin_star, ycenter=y_nominal - 0, origin=origin_star,
tar_spec=normalSED, tar_spec=normalSED,
band_start=brange[0], band_end=brange[1], band_start=brange[0], band_end=brange[1],
conf=chip.sls_conf[1], conf=chip.sls_conf[1],
isAlongY=0, isAlongY=0,
flat_cube=flat_cube) flat_cube=flat_cube)
self.addSLStoChipImage(sdp=sdp, chip=chip, xOrderSigPlus=xOrderSigPlus, local_wcs=chip_wcs_local) # 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], 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, psf_model=psf_model, bandNo=i + 1, grating_split_pos=grating_split_pos,
# local_wcs=chip_wcs_local, pos_img=pos_img) local_wcs=chip_wcs_local, pos_img=pos_img)
del sdp del sdp
elif grating_split_pos_chip >= gal_end[1]: elif grating_split_pos_chip >= gal_end[1]:
sdp = SpecDisperser(orig_img=starImg_List, xcenter=x_nominal - 0, sdp = SpecDisperser(orig_img=starImg, xcenter=x_nominal - 0,
ycenter=y_nominal - 0, origin=origin_star, ycenter=y_nominal - 0, origin=origin_star,
tar_spec=normalSED, tar_spec=normalSED,
band_start=brange[0], band_end=brange[1], band_start=brange[0], band_end=brange[1],
conf=chip.sls_conf[0], conf=chip.sls_conf[0],
isAlongY=0, isAlongY=0,
flat_cube=flat_cube) flat_cube=flat_cube)
self.addSLStoChipImage(sdp=sdp, chip=chip, xOrderSigPlus=xOrderSigPlus, local_wcs=chip_wcs_local) # 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], 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, psf_model=psf_model, bandNo=i + 1, grating_split_pos=grating_split_pos,
# local_wcs=chip_wcs_local, pos_img=pos_img) local_wcs=chip_wcs_local, pos_img=pos_img)
del sdp del sdp
# del psf # del psf
return 1, pos_shear return 1, pos_shear
......
observation_sim/mock_objects/SpecDisperser/SpecDisperser.py
View file @ 64eb7956
...@@ -65,30 +65,10 @@ class SpecDisperser(object): ...@@ -65,30 +65,10 @@ class SpecDisperser(object):
# self.img_x = orig_img.shape[1] # self.img_x = orig_img.shape[1]
# self.img_y = orig_img.shape[0] # self.img_y = orig_img.shape[0]
# 5 orders, A, B , self.thumb_img = np.abs(orig_img.array)
orderName=["A","B","C","D","E"] self.thumb_x = orig_img.center.x
self.orig_img_orders = OrderedDict() self.thumb_y = orig_img.center.y
if isinstance(orig_img, list): self.img_sh = orig_img.array.shape
orig_img_list = orig_img
list_len = len(orig_img_list)
if list_len < 5:
for i in np.arange(5-list_len):
orig_img_list.append(orig_img_list[list_len-1])
for i, k in enumerate(orig_img_list):
self.orig_img_orders[orderName[i]] = k
if isinstance(orig_img, galsim.Image):
for i in np.arange(5):
self.orig_img_orders[orderName[i]] = orig_img
orig_img_one = self.orig_img_orders["A"]
self.thumb_img = np.abs(orig_img_one.array)
self.thumb_x = orig_img_one.center.x
self.thumb_y = orig_img_one.center.y
self.img_sh = orig_img_one.array.shape
self.id = gid self.id = gid
...@@ -98,13 +78,10 @@ class SpecDisperser(object): ...@@ -98,13 +78,10 @@ class SpecDisperser(object):
self.isAlongY = isAlongY self.isAlongY = isAlongY
self.flat_cube = flat_cube self.flat_cube = flat_cube
if self.isAlongY == 1: if self.isAlongY == 1:
for order in orderName: self.thumb_img, self.thumb_x, self.thumb_y = rotate90(array_orig=self.thumb_img, xc=orig_img.center.x,
self.orig_img_orders[order], self.thumb_x, self.thumb_y = rotate90(array_orig=self.orig_img_orders[order], xc=orig_img_one.center.x, yc=orig_img.center.y, isClockwise=1)
yc=orig_img_one.center.y, isClockwise=1)
# self.thumb_img, self.thumb_x, self.thumb_y = rotate90(array_orig=self.thumb_img, xc=orig_img_one.center.x,
# yc=orig_img_one.center.y, isClockwise=1)
self.img_sh = self.orig_img_orders[order].array.T.shape self.img_sh = orig_img.array.T.shape
self.xcenter = ycenter self.xcenter = ycenter
self.ycenter = xcenter self.ycenter = xcenter
...@@ -134,16 +111,10 @@ class SpecDisperser(object): ...@@ -134,16 +111,10 @@ class SpecDisperser(object):
def compute_spec(self, beam): def compute_spec(self, beam):
# if beam == "B":
# return self.thumb_img, self.origin[1], self.origin[0], None, None, None
from .disperse_c import interp from .disperse_c import interp
from .disperse_c import disperse from .disperse_c import disperse
# from MockObject.disperse_c import disperse # from MockObject.disperse_c import disperse
self.thumb_img = np.abs(self.orig_img_orders[beam].array)
self.thumb_x = self.orig_img_orders[beam].center.x
self.thumb_y = self.orig_img_orders[beam].center.y
self.img_sh = self.orig_img_orders[beam].array.shape
dx = self.grating_conf.dxlam[beam] dx = self.grating_conf.dxlam[beam]
xoff = 0 xoff = 0
ytrace_beam, lam_beam = self.grating_conf.get_beam_trace(x=self.xcenter, y=self.ycenter, dx=(dx + xoff), ytrace_beam, lam_beam = self.grating_conf.get_beam_trace(x=self.xcenter, y=self.ycenter, dx=(dx + xoff),
...@@ -198,8 +169,7 @@ class SpecDisperser(object): ...@@ -198,8 +169,7 @@ class SpecDisperser(object):
dyc = cast[int](np.floor(ytrace_beam+0.5)) dyc = cast[int](np.floor(ytrace_beam+0.5))
# dypix = cast[int](np.floor(ytrace_beam - dyc[0] + x0[0] + 0.5)) dypix = cast[int](np.floor(ytrace_beam - dyc[0] + x0[0] + 0.5))
dypix = dyc - dyc[0] + x0[0]
frac_ids = yfrac_beam < 0 frac_ids = yfrac_beam < 0
...@@ -278,8 +248,7 @@ class SpecDisperser(object): ...@@ -278,8 +248,7 @@ class SpecDisperser(object):
# beam_flat[k] = self.flat_cube[:, originOut_y + i, originOut_x + j] # beam_flat[k] = self.flat_cube[:, originOut_y + i, originOut_x + j]
status = disperse.disperse_grism_object(self.thumb_img.astype(np.float32), status = disperse.disperse_grism_object(self.thumb_img.astype(np.float32),
flat_index[nonz], flat_index[nonz], yfrac_beam[nonz],
yfrac_beam[nonz],
sensitivity_beam[nonz], sensitivity_beam[nonz],
modelf, x0, modelf, x0,
array(self.img_sh, array(self.img_sh,
...@@ -289,24 +258,11 @@ class SpecDisperser(object): ...@@ -289,24 +258,11 @@ class SpecDisperser(object):
lam_beam[lam_index][nonz]) lam_beam[lam_index][nonz])
model = modelf.reshape(beam_sh) model = modelf.reshape(beam_sh)
# n1 = np.sum(np.isinf(model))
# n2 = np.sum(np.isnan(model))
# n3 = np.sum(np.isinf(modelf))
# n4 = np.sum(np.isnan(modelf))
# if n1>0 or n2 > 0:
# print("DEBUG: SpecDisperser, inf:%d, nan:%d--------%d,%d"%(n1, n2, n3, n4))
# print(dypix)
# n1 = np.sum(np.isinf(self.thumb_img.astype(np.float32)))
# n2 = np.sum(np.isnan(self.thumb_img.astype(np.float32)))
# n3 = np.sum(np.isinf(yfrac_beam))
# n4 = np.sum(np.isnan(yfrac_beam))
# n5 = np.sum(np.isinf(sensitivity_beam))
# n6 = np.sum(np.isnan(sensitivity_beam))
# print("DEBUG: SpecDisperser, innput ---inf:%d, nan:%d, yfrac_beam:%d/%d, sensitivity_beam:%d/%d"%(n1, n2, n3, n4, n5, n6))
self.beam_flux[beam] = sum(modelf) self.beam_flux[beam] = sum(modelf)
if self.isAlongY == 1: if self.isAlongY == 1:
model, _, _ = rotate90(array_orig=model, isClockwise=0) model, _, _ = rotate90(array_orig=model, isClockwise=0)
return model, originOut_x, originOut_y, dxpix, dypix, lam_beam, ysens return model, originOut_x, originOut_y, dxpix, dypix, lam_beam, ysens
def writerSensitivityFile(self, conffile='', beam='', w=None, sens=None): def writerSensitivityFile(self, conffile='', beam='', w=None, sens=None):
......
observation_sim/mock_objects/SpecDisperser/disperse_c/disperse.pyx
View file @ 64eb7956
...@@ -21,29 +21,6 @@ cdef extern from "math.h": ...@@ -21,29 +21,6 @@ cdef extern from "math.h":
double sqrt(double x) double sqrt(double x)
double exp(double x) double exp(double x)
def check_nan2D(np.ndarray[FTYPE_t, ndim=2] arr):
cdef int i, j
cdef int nrows = arr.shape[0]
cdef int ncols = arr.shape[1]
# 遍历数组的每个元素并检查是否存在 NaN
for i in range(nrows):
for j in range(ncols):
if np.isnan(arr[i, j]) | np.isinf(arr[i, j]):
return True
return False
def check_nan1d(np.ndarray[DTYPE_t, ndim=1] arr):
cdef int i
cdef int n = arr.shape[0]
# 遍历数组的每个元素并检查是否存在 NaN
for i in range(n):
if np.isnan(arr[i]) | np.isinf(arr[i]):
return True
return False
@cython.boundscheck(False) @cython.boundscheck(False)
@cython.wraparound(False) @cython.wraparound(False)
@cython.embedsignature(True) @cython.embedsignature(True)
...@@ -77,18 +54,6 @@ def disperse_grism_object(np.ndarray[FTYPE_t, ndim=2] flam, ...@@ -77,18 +54,6 @@ def disperse_grism_object(np.ndarray[FTYPE_t, ndim=2] flam,
nk = len(idxl) nk = len(idxl)
nl = len(full) nl = len(full)
#if check_nan2D(flam):
# print("DEBUG: disperse, input Array 'flam' contains NaN.")
#if check_nan1d(ysens):
# print("DEBUG: disperse, input Array 'ysens' contains NaN.")
#if check_nan1d(yfrac):
# print("DEBUG: disperse, input Array 'yfrac' contains NaN.")
#if check_nan1d(full):
# print("DEBUG: disperse, input Array 'full' contains NaN before processing.")
if (flat is not None): if (flat is not None):
nlamb = len(wlambda) nlamb = len(wlambda)
nflat = len(flat) nflat = len(flat)
...@@ -130,15 +95,14 @@ def disperse_grism_object(np.ndarray[FTYPE_t, ndim=2] flam, ...@@ -130,15 +95,14 @@ def disperse_grism_object(np.ndarray[FTYPE_t, ndim=2] flam,
else: else:
for i in range(0-x0[1], x0[1]): for i in range(0-x0[1], x0[1]):
x_pos = x0[1]+i if (x0[1]+i < 0) | (x0[1]+i >= shd[1]):
if (x_pos < 0) | (x_pos >= shd[1]):
continue continue
for j in range(0-x0[0], x0[0]): for j in range(0-x0[0], x0[0]):
y_pos = x0[0]+j if (x0[0]+j < 0) | (x0[0]+j >= shd[0]):
if (y_pos < 0) | (y_pos >= shd[0]):
continue continue
fl_ij = flam[y_pos, x_pos] #/1.e-17
fl_ij = flam[x0[0]+j, x0[1]+i] #/1.e-17
if (fl_ij == 0): if (fl_ij == 0):
continue continue
...@@ -146,14 +110,11 @@ def disperse_grism_object(np.ndarray[FTYPE_t, ndim=2] flam, ...@@ -146,14 +110,11 @@ def disperse_grism_object(np.ndarray[FTYPE_t, ndim=2] flam,
k1 = idxl[k]+j*shg[1]+i k1 = idxl[k]+j*shg[1]+i
if (k1 >= 0) & (k1 < nl): if (k1 >= 0) & (k1 < nl):
full[k1] += ysens[k]*fl_ij*(1-yfrac[k]) full[k1] += ysens[k]*fl_ij*(1-yfrac[k])
k2 = idxl[k]+(j+1)*shg[1]+i k2 = idxl[k]+(j+1)*shg[1]+i
if (k2 >= 0) & (k2 < nl): if (k2 >= 0) & (k2 < nl):
full[k2] += ysens[k]*fl_ij*yfrac[k] full[k2] += ysens[k]*fl_ij*yfrac[k]
#if (check_nan1d(full)):
# print("DEBUG: disperse, output Array 'full' contains NaN after processing.+++++++++++++++++++++++++++")
return True return True
@cython.boundscheck(False) @cython.boundscheck(False)
......
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