diff --git a/ObservationSim/Config/ChipOutput.py b/ObservationSim/Config/ChipOutput.py index 6abe5783162d9ae2a5f992613007cb52a9ade758..c4d2421b6e84321dc109f542c1d56735ac874da7 100755 --- a/ObservationSim/Config/ChipOutput.py +++ b/ObservationSim/Config/ChipOutput.py @@ -73,7 +73,7 @@ class ChipOutput(object): return hdr # def cat_add_obj(self, obj, pos_img, snr, pos_shear, g1, g2): - def cat_add_obj(self, obj, pos_img, pos_shear, g1, g2): + def cat_add_obj(self, obj, pos_img, pos_shear, g1, g2): ximg = pos_img.x - self.chip.bound.xmin + 1.0 yimg = pos_img.y - self.chip.bound.ymin + 1.0 e1, e2, g1, g2, e1OBS, e2OBS = obj.getObservedEll(g1, g2) @@ -90,4 +90,4 @@ class ChipOutput(object): line = self.fmt%(obj.id, int(self.chipLabel), self.filt.filter_type, ximg, yimg, obj.ra, obj.dec, obj.z, obj.getMagFilter(self.filt), obj.param["star"], 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, pos_shear.g1, pos_shear.g2, e1, e2, g1, g2, e1OBS, e2OBS, 0, 0.0, 0.0, obj.param['model_tag'], obj.param['teff'], obj.param['logg'],obj.param['feh']) # line = self.fmt%(obj.id, int(self.chipLabel), self.filt.filter_type, ximg, yimg, obj.ra, obj.dec, obj.z, obj.getMagFilter(self.filt), obj.param["star"], pos_shear.g1, pos_shear.g2, e1, e2, g1, g2, e1OBS, e2OBS) - self.cat.write(line) \ No newline at end of file + self.cat.write(line) diff --git a/ObservationSim/MockObject/Galaxy.py b/ObservationSim/MockObject/Galaxy.py index 3be0dbe4f386a9fa97885d1626a8a0658ecd7212..a7ec576821fd818b6f85a21f0b365f495f4d1244 100755 --- a/ObservationSim/MockObject/Galaxy.py +++ b/ObservationSim/MockObject/Galaxy.py @@ -3,7 +3,7 @@ import galsim import os, sys import astropy.constants as cons from astropy.table import Table -from ._util import eObs, integrate_sed_bandpass, getNormFactorForSpecWithABMAG, getObservedSED, getABMAG +from ._util import eObs, integrate_sed_bandpass, getNormFactorForSpecWithABMAG, getObservedSED, getABMAG,convolveGaussXorders from .SpecDisperser import SpecDisperser from .MockObject import MockObject from scipy import interpolate @@ -279,7 +279,7 @@ class Galaxy(MockObject): folding_threshold = 5.e-3 gsp = galsim.GSParams(folding_threshold=folding_threshold) # nphotons_sum = 0 - + xOrderSigPlus = {'A':1.3909419820029296,'B':1.4760376591236062,'C':4.035447379743442,'D':5.5684364343742825,'E':16.260021029735388} for i in range(len(bandpass_list)): bandpass = bandpass_list[i] @@ -314,8 +314,9 @@ class Galaxy(MockObject): 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] + img_s,orig_off = convolveGaussXorders(img_s,xOrderSigPlus[k]) + origin_order_x = v[1]-orig_off + origin_order_y = v[2]-orig_off specImg = galsim.ImageF(img_s) photons = galsim.PhotonArray.makeFromImage(specImg) photons.x += origin_order_x @@ -413,4 +414,4 @@ class Galaxy(MockObject): 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 \ No newline at end of file + return self.e1_total, self.e2_total, g1, g2, e1_obs, e2_obs diff --git a/ObservationSim/MockObject/MockObject.py b/ObservationSim/MockObject/MockObject.py index 13dcac8572043297a30c54f868d29b645313e4cf..2ac3fb106bc3dec200b240fc6d67786a27546050 100755 --- a/ObservationSim/MockObject/MockObject.py +++ b/ObservationSim/MockObject/MockObject.py @@ -2,7 +2,7 @@ import galsim import numpy as np import astropy.constants as cons from astropy.table import Table -from ._util import magToFlux, vc_A +from ._util import magToFlux, vc_A, convolveGaussXorders from ._util import integrate_sed_bandpass, getNormFactorForSpecWithABMAG, getObservedSED, getABMAG from .SpecDisperser import SpecDisperser @@ -190,6 +190,8 @@ class MockObject(object): normalSED = Table(np.array([self.sed['WAVELENGTH'], self.sed['FLUX'] * sedNormFactor]).T, names=('WAVELENGTH', 'FLUX')) + + xOrderSigPlus = {'A':1.3909419820029296,'B':1.4760376591236062,'C':4.035447379743442,'D':5.5684364343742825,'E':16.260021029735388} for i in range(len(bandpass_list)): bandpass = bandpass_list[i] psf, pos_shear = psf_model.get_PSF(chip=chip, pos_img=pos_img, bandpass=bandpass, folding_threshold=folding_threshold) @@ -211,8 +213,9 @@ class MockObject(object): 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] + img_s,orig_off = convolveGaussXorders(img_s,xOrderSigPlus[k]) + origin_order_x = v[1]-orig_off + origin_order_y = v[2]-orig_off specImg = galsim.ImageF(img_s) photons = galsim.PhotonArray.makeFromImage(specImg) photons.x += origin_order_x @@ -247,4 +250,4 @@ class MockObject(object): return snr_obj def getObservedEll(self, g1=0, g2=0): - return 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 \ No newline at end of file + return 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 diff --git a/ObservationSim/MockObject/SpecDisperser/SpecDisperser.py b/ObservationSim/MockObject/SpecDisperser/SpecDisperser.py index 3321b16c1631c0eb2a269f7e7cc3c73d0accc1ad..9a6a64df6aeceb439d97764cbc66516b9734f55b 100644 --- a/ObservationSim/MockObject/SpecDisperser/SpecDisperser.py +++ b/ObservationSim/MockObject/SpecDisperser/SpecDisperser.py @@ -101,6 +101,7 @@ class SpecDisperser(object): beam_names = self.grating_conf.beams for beam in beam_names: + #for beam in ['A','B']: all_orders[beam] = self.compute_spec(beam) return all_orders @@ -123,8 +124,8 @@ class SpecDisperser(object): 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.5)) - + 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']) @@ -196,7 +197,7 @@ class SpecDisperser(object): return model, originOut_x, originOut_y def writerSensitivityFile(self, conffile = '', beam = '', w = None, sens = None): - orders={'A':'1st','B':'0st','C':'2st'} + 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')) diff --git a/ObservationSim/MockObject/SpecDisperser/disperse_c/interp.pyx b/ObservationSim/MockObject/SpecDisperser/disperse_c/interp.pyx index 5d691856705985ad0e404843d2c3d4f595fe201d..32d0f3fd1ed58269f43a3004f12b1344f016ccf2 100644 --- a/ObservationSim/MockObject/SpecDisperser/disperse_c/interp.pyx +++ b/ObservationSim/MockObject/SpecDisperser/disperse_c/interp.pyx @@ -189,7 +189,7 @@ def interp_conserve_c(np.ndarray[DTYPE_t, ndim=1] x, np.ndarray[DTYPE_t, ndim=1] if (templmid[k+1] == tlam[i]) & (i < ntlam): h = tlam[i]-tlam[i-1]; numsum+=h*(tf[i]+tf[i-1]); - else: + elif (i < ntlam): i-=1; h = templmid[k+1]-tlam[i]; numsum+=h*(tempfmid[k+1]+tf[i]); diff --git a/ObservationSim/MockObject/_util.py b/ObservationSim/MockObject/_util.py index 00d13897c88cf013e0fdb74bbb109a008eeca800..f957e7b0b9367cc1e6178cbe43d259d9818c3cda 100755 --- a/ObservationSim/MockObject/_util.py +++ b/ObservationSim/MockObject/_util.py @@ -545,4 +545,22 @@ def tag_sed(h5file, model_tag, teff=5000, logg=2, feh=0): 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 \ No newline at end of file + return path, wave, flux + +from astropy.modeling.models import Gaussian2D +from scipy import signal +def convolveGaussXorders(img=None, sigma = 1): + + 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 + +