diff --git a/ObservationSim/Config/ChipOutput.py b/ObservationSim/Config/ChipOutput.py index 36da1803de00642ced3cbaa94fba02cf79737715..d9253fd03936e848b026d2fe065a86611c909ec6 100755 --- a/ObservationSim/Config/ChipOutput.py +++ b/ObservationSim/Config/ChipOutput.py @@ -68,8 +68,13 @@ class ChipOutput(object): self.cat.write(self.hdr) def cat_add_obj(self, obj, pos_img, pos_shear): - ximg = pos_img.x - self.chip.bound.xmin + 1.0 - yimg = pos_img.y - self.chip.bound.ymin + 1.0 + # ximg = pos_img.x - self.chip.bound.xmin + 1.0 + # yimg = pos_img.y - self.chip.bound.ymin + 1.0 + # print('-------------debug-----------------') + # print('from global',ximg, yimg) + ximg = obj.real_pos.x + 1.0 + yimg = obj.real_pos.y + 1.0 + # print('from loacl',ximg, yimg) line = self.fmt%( obj.id, int(self.chipLabel), self.filt.filter_type, ximg, yimg, obj.ra, obj.dec, obj.ra_orig, obj.dec_orig, obj.z, obj.getMagFilter(self.filt), obj.type, diff --git a/ObservationSim/Config/Header/ImageHeader.py b/ObservationSim/Config/Header/ImageHeader.py index 75229b718b932cc73623ddb0b50cc02ac3ec2ff9..86a46cbc4f634ef8d7d1000e66665c3965f8d5bd 100644 --- a/ObservationSim/Config/Header/ImageHeader.py +++ b/ObservationSim/Config/Header/ImageHeader.py @@ -7,7 +7,7 @@ from astropy.io import fits import astropy.wcs as pywcs from collections import OrderedDict -from scipy import math +# from scipy import math import random import os @@ -87,6 +87,17 @@ def rotate_CD_matrix(cd, pa_aper): cd_rot = np.dot(mat, cd) return cd_rot +def calcaluteSLSRotSkyCoor(pix_xy = None,rot_angle = 1, xlen = 9216, ylen = 9232, w = None): + rad = np.deg2rad(rot_angle) + mat = np.zeros((2,2)) + mat[0,:] = np.array([np.cos(rad),-np.sin(rad)]) + mat[1,:] = np.array([np.sin(rad),np.cos(rad)]) + center = np.array([xlen/2, ylen/2]) + rot_pix = np.dot(mat, pix_xy-center) + center + skyCoor = w.wcs_pix2world(np.array([rot_pix]), 1) + + return skyCoor + # def Header_extention(xlen = 9216, ylen = 9232, gain = 1.0, readout = 5.0, dark = 0.02,saturation=90000, row_num = 1, col_num = 1): # @@ -153,7 +164,7 @@ def rotate_CD_matrix(cd, pa_aper): ##9232 9216 898 534 1309 60 -40 -23.4333 -def WCS_def(xlen = 9216, ylen = 9232, gapy = 898.0, gapx1 = 534, gapx2 = 1309, ra = 60, dec = -40, pa = -23.433,psize = 0.074, row_num = 1, col_num = 1): +def WCS_def(xlen = 9216, ylen = 9232, gapy = 898.0, gapx1 = 534, gapx2 = 1309, ra = 60, dec = -40, pa = -23.433,psize = 0.074, row_num = 1, col_num = 1, filter = 'GI'): """ Creat a wcs frame for CCST with multiple extensions @@ -228,6 +239,63 @@ def WCS_def(xlen = 9216, ylen = 9232, gapy = 898.0, gapx1 = 534, gapx2 = 1309, r r_dat['CD1_2'] = cd_rot[0,1] r_dat['CD2_1'] = cd_rot[1,0] r_dat['CD2_2'] = cd_rot[1,1] + + if filter in ['GU', 'GV', 'GI']: + from astropy import wcs + + w = wcs.WCS(naxis=2) + w.wcs.crpix = [r_dat['CRPIX1'], r_dat['CRPIX2']] + w.wcs.cd = cd_rot + w.wcs.crval = [ra_ref, dec_ref] + w.wcs.ctype = [r_dat['CTYPE1'], r_dat['CTYPE2']] + + # test_center_o = w.wcs_pix2world(np.array([[xlen / 2, ylen / 2]]), 1) + + sls_rot = 1 + if i > 2: + sls_rot = -sls_rot + + sn_x = 30 + sn_y = 30 + x_pixs = np.zeros(sn_y * sn_x) + y_pixs = np.zeros(sn_y * sn_x) + xpixs_line = np.linspace(1, xlen, sn_x) + ypixs_line = np.linspace(1, ylen, sn_y) + + sky_coors = [] + + for n1, y in enumerate(ypixs_line): + for n2, x in enumerate(xpixs_line): + i_pix = n1 * sn_x + n2 + x_pixs[i_pix] = x + y_pixs[i_pix] = y + + pix_coor = np.array([x, y]) + sc1 = calcaluteSLSRotSkyCoor(pix_xy=pix_coor, rot_angle=sls_rot, w=w) + # print(sc1[0,0],sc1[0,1]) + sky_coors.append((sc1[0, 0], sc1[0, 1])) + + from astropy.coordinates import SkyCoord + from astropy.wcs.utils import fit_wcs_from_points + + wcs_new = fit_wcs_from_points(xy=np.array([x_pixs, y_pixs]), + world_coords=SkyCoord(sky_coors, frame="icrs", unit="deg"), projection='TAN') + + # print(wcs_new) + # test_center = wcs_new.wcs_pix2world(np.array([[xlen / 2, ylen / 2]]), 1) + # + # print(test_center - test_center_o) + + r_dat['CD1_1'] = wcs_new.wcs.cd[0, 0] + r_dat['CD1_2'] = wcs_new.wcs.cd[0, 1] + r_dat['CD2_1'] = wcs_new.wcs.cd[1, 0] + r_dat['CD2_2'] = wcs_new.wcs.cd[1, 1] + r_dat['CRPIX1'] = wcs_new.wcs.crpix[0] + r_dat['CRPIX2'] = wcs_new.wcs.crpix[1] + + r_dat['CRVAL1'] = wcs_new.wcs.crval[0] + r_dat['CRVAL2'] = wcs_new.wcs.crval[1] + return r_dat @@ -347,7 +415,7 @@ def generateExtensionHeader(xlen = 9216, ylen = 9232,ra = 60, dec = -40, pa = -2 h_ext['CCDCHIP'] = 'ccd' + CCDID[k-1].rjust(2,'0') h_ext['POS_ANG'] = pa header_wcs = WCS_def(xlen=xlen, ylen=ylen, gapy=898.0, gapx1=534, gapx2=1309, ra=ra, dec=dec, pa=pa, psize=psize, - row_num=row_num, col_num=col_num) + row_num=row_num, col_num=col_num, filter = h_ext['FILTER']) h_ext['CRPIX1'] = header_wcs['CRPIX1'] h_ext['CRPIX2'] = header_wcs['CRPIX2'] diff --git a/ObservationSim/Instrument/Chip/Chip.py b/ObservationSim/Instrument/Chip/Chip.py index 69d0c1038919febec2958d8e3bb1b5893b1dca3d..96aa78fe5272fd598ebebaf00d6b1c739cfd21de 100755 --- a/ObservationSim/Instrument/Chip/Chip.py +++ b/ObservationSim/Instrument/Chip/Chip.py @@ -82,6 +82,8 @@ class Chip(FocalPlane): else: self.sensor = galsim.Sensor() + self.flat_cube = None # for spectroscopic flat field cube simulation + # def _getChipRowCol(self): # self.rowID = (self.chipID - 1) // 5 + 1 # self.colID = (self.chipID - 1) % 5 + 1 @@ -846,3 +848,20 @@ class Chip(FocalPlane): # sub_img.write("%s/%s_%s.fits" % (chip_output.subdir, img_name_root, rowcoltag)) # del sub_img return img + + def loadSLSFLATCUBE(self, flat_fn='flat_cube.fits'): + from astropy.io import fits + try: + with pkg_resources.files('ObservationSim.Instrument.data').joinpath(flat_fn) as data_path: + flat_fits = fits.open(data_path, ignore_missing_simple=True) + except AttributeError: + with pkg_resources.path('ObservationSim.Instrument.data', flat_fn) as data_path: + flat_fits = fits.open(data_path, ignore_missing_simple=True) + + fl = len(flat_fits) + fl_sh = flat_fits[0].data.shape + assert fl == 4, 'FLAT Field Cube is Not 4 layess!!!!!!!' + self.flat_cube = np.zeros([fl, fl_sh[0], fl_sh[1]]) + for i in np.arange(0, fl, 1): + self.flat_cube[i] = flat_fits[i].data + diff --git a/ObservationSim/Instrument/data/flatCube/flat_cube.fits b/ObservationSim/Instrument/data/flatCube/flat_cube.fits new file mode 100644 index 0000000000000000000000000000000000000000..cca7b145abfca4b8f012a728f588d318d3b734f4 Binary files /dev/null and b/ObservationSim/Instrument/data/flatCube/flat_cube.fits differ diff --git a/ObservationSim/MockObject/Galaxy.py b/ObservationSim/MockObject/Galaxy.py index 5b2cb4488a2689a60e6717df862c1a8deed6e22e..6965224810b2e29ea9587d028c5618ac02dd26e5 100755 --- a/ObservationSim/MockObject/Galaxy.py +++ b/ObservationSim/MockObject/Galaxy.py @@ -107,6 +107,18 @@ class Galaxy(MockObject): 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.real_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) + + real_wcs_local = self.real_wcs.local(self.real_pos) + for i in range(len(bandpass_list)): bandpass = bandpass_list[i] @@ -150,20 +162,29 @@ class Galaxy(MockObject): gal = galsim.Convolve(psf, gal) # Use (explicit) stamps to draw - stamp = gal.drawImage(wcs=self.localWCS, method='phot', offset=self.offset, save_photons=True) + # stamp = gal.drawImage(wcs=self.localWCS, method='phot', offset=self.offset, save_photons=True) + # xmax = max(xmax, stamp.xmax) + # ymax = max(ymax, stamp.ymax) + # photons = stamp.photons + # photons.x += self.x_nominal + # photons.y += self.y_nominal + # photons_list.append(photons) + + stamp = gal.drawImage(wcs=real_wcs_local, method='phot', offset=offset, save_photons=True) xmax = max(xmax, stamp.xmax) ymax = max(ymax, stamp.ymax) photons = stamp.photons - photons.x += self.x_nominal - photons.y += self.y_nominal + photons.x += x_nominal + photons.y += y_nominal photons_list.append(photons) # print('xmax = %d, ymax = %d '%(xmax, ymax)) - stamp = galsim.ImageF(int(xmax*1.1), int(ymax*1.1)) - stamp.wcs = self.localWCS - stamp.setCenter(self.x_nominal, self.y_nominal) - bounds = stamp.bounds & chip.img.bounds + stamp = galsim.ImageF(int(xmax * 1.1), int(ymax * 1.1)) + stamp.wcs = real_wcs_local + stamp.setCenter(x_nominal, y_nominal) + bounds = stamp.bounds & galsim.BoundsI(0, chip.npix_x - 1, 0, chip.npix_y - 1) + chip.img.setOrigin(0, 0) stamp[bounds] = chip.img[bounds] if not big_galaxy: @@ -180,10 +201,36 @@ class Galaxy(MockObject): else: sensor.accumulate(photons_list[i], stamp, resume=True) - # print(stamp.array.sum()) - # chip.img[bounds] += stamp[bounds] chip.img[bounds] = stamp[bounds] - # print("nphotons_sum = ", nphotons_sum) + + chip.img.setOrigin(chip.bound.xmin, chip.bound.ymin) + + + + # stamp = galsim.ImageF(int(xmax*1.1), int(ymax*1.1)) + # stamp.wcs = self.localWCS + # stamp.setCenter(self.x_nominal, self.y_nominal) + # bounds = stamp.bounds & chip.img.bounds + # stamp[bounds] = chip.img[bounds] + # + # if not big_galaxy: + # for i in range(len(photons_list)): + # if i == 0: + # chip.sensor.accumulate(photons_list[i], stamp) + # else: + # chip.sensor.accumulate(photons_list[i], stamp, resume=True) + # else: + # sensor = galsim.Sensor() + # for i in range(len(photons_list)): + # if i == 0: + # sensor.accumulate(photons_list[i], stamp) + # else: + # sensor.accumulate(photons_list[i], stamp, resume=True) + # + # # print(stamp.array.sum()) + # # chip.img[bounds] += stamp[bounds] + # chip.img[bounds] = stamp[bounds] + # # print("nphotons_sum = ", nphotons_sum) del photons_list del stamp return True, pos_shear @@ -201,6 +248,19 @@ class Galaxy(MockObject): 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.real_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) + + real_wcs_local = self.real_wcs.local(self.real_pos) + + big_galaxy = False if self.hlr_disk > 3.0: # Very big galaxy big_galaxy = True @@ -211,8 +271,11 @@ class Galaxy(MockObject): 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 = chip.bound.xmin + grating_split_pos + grating_split_pos_chip = 0 + grating_split_pos for i in range(len(bandpass_list)): bandpass = bandpass_list[i] @@ -236,10 +299,10 @@ class Galaxy(MockObject): gal = gal.shear(gal_shear) gal = galsim.Convolve(psf, gal) - starImg = gal.drawImage(wcs=self.localWCS) + starImg = gal.drawImage(wcs=real_wcs_local) - origin_star = [self.y_nominal - (starImg.center.y - starImg.ymin), - self.x_nominal - (starImg.center.x - starImg.xmin)] + origin_star = [y_nominal - (starImg.center.y - starImg.ymin), + x_nominal - (starImg.center.x - starImg.xmin)] 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] @@ -251,52 +314,56 @@ class Galaxy(MockObject): subImg_p1 = starImg.array[:, 0:subSlitPos] star_p1 = galsim.Image(subImg_p1) origin_p1 = origin_star - xcenter_p1 = min(self.x_nominal,grating_split_pos_chip-1) - chip.bound.xmin - ycenter_p1 = self.y_nominal-chip.bound.ymin + 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=bandpass.blue_limit * 10, band_end=bandpass.red_limit * 10, conf=chip.sls_conf[0], - isAlongY=0) + isAlongY=0, + flat_cube=flat_cube) - self.addSLStoChipImage(sdp=sdp_p1, chip=chip, xOrderSigPlus = xOrderSigPlus) + self.addSLStoChipImage(sdp=sdp_p1, chip=chip, xOrderSigPlus = xOrderSigPlus, local_wcs=real_wcs_local) subImg_p2 = starImg.array[:, subSlitPos+1:starImg.array.shape[1]] star_p2 = galsim.Image(subImg_p2) origin_p2 = [origin_star[0], grating_split_pos_chip] - xcenter_p2 = max(self.x_nominal, grating_split_pos_chip - 1) - chip.bound.xmin - ycenter_p2 = self.y_nominal - chip.bound.ymin + 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=bandpass.blue_limit * 10, band_end=bandpass.red_limit * 10, conf=chip.sls_conf[1], - isAlongY=0) + isAlongY=0, + flat_cube=flat_cube) - self.addSLStoChipImage(sdp=sdp_p2, chip=chip, xOrderSigPlus = xOrderSigPlus) + self.addSLStoChipImage(sdp=sdp_p2, chip=chip, xOrderSigPlus = xOrderSigPlus, local_wcs=real_wcs_local) del sdp_p1 del sdp_p2 elif grating_split_pos_chip<=gal_origin[1]: - sdp = SpecDisperser(orig_img=starImg, xcenter=self.x_nominal - chip.bound.xmin, - ycenter=self.y_nominal - chip.bound.ymin, origin=origin_star, + sdp = SpecDisperser(orig_img=starImg, xcenter=x_nominal - 0, + ycenter=y_nominal - 0, origin=origin_star, tar_spec=normalSED, band_start=bandpass.blue_limit * 10, band_end=bandpass.red_limit * 10, conf=chip.sls_conf[1], - isAlongY=0) - self.addSLStoChipImage(sdp=sdp, chip=chip, xOrderSigPlus = xOrderSigPlus) + isAlongY=0, + flat_cube=flat_cube) + self.addSLStoChipImage(sdp=sdp, chip=chip, xOrderSigPlus = xOrderSigPlus, local_wcs=real_wcs_local) del sdp elif grating_split_pos_chip>=gal_end[1]: - sdp = SpecDisperser(orig_img=starImg, xcenter=self.x_nominal - chip.bound.xmin, - ycenter=self.y_nominal - chip.bound.ymin, origin=origin_star, + sdp = SpecDisperser(orig_img=starImg, xcenter=x_nominal - 0, + ycenter=y_nominal - 0, origin=origin_star, tar_spec=normalSED, band_start=bandpass.blue_limit * 10, band_end=bandpass.red_limit * 10, conf=chip.sls_conf[0], - isAlongY=0) - self.addSLStoChipImage(sdp=sdp, chip=chip, xOrderSigPlus = xOrderSigPlus) + isAlongY=0, + flat_cube=flat_cube) + self.addSLStoChipImage(sdp=sdp, chip=chip, xOrderSigPlus = xOrderSigPlus, local_wcs=real_wcs_local) del sdp # print(self.y_nominal, starImg.center.y, starImg.ymin) diff --git a/ObservationSim/MockObject/MockObject.py b/ObservationSim/MockObject/MockObject.py index 4bdb2a9a9859f32590c7d9140e5c044109dc316e..4efea6700d51d5aa30915eeb3b327a3f5fd91dbd 100755 --- a/ObservationSim/MockObject/MockObject.py +++ b/ObservationSim/MockObject/MockObject.py @@ -4,9 +4,11 @@ import astropy.constants as cons from astropy.table import Table from ObservationSim.MockObject._util import magToFlux, VC_A, convolveGaussXorders -from ObservationSim.MockObject._util import integrate_sed_bandpass, getNormFactorForSpecWithABMAG, getObservedSED, getABMAG +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 @@ -17,7 +19,7 @@ class MockObject(object): self.type = "star" elif self.param["star"] == 2: self.type = "quasar" - + self.id = self.param["id"] self.ra = self.param["ra"] self.dec = self.param["dec"] @@ -51,15 +53,15 @@ class MockObject(object): 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] + return self.param["mag_%s" % filt.filter_type] # (TEST) stamp size # return 13.0 - + def getFluxFilter(self, filt): - return self.param["flux_%s"%filt.filter_type] + return self.param["flux_%s" % filt.filter_type] def getNumPhotons(self, flux, tel, exptime=150.): - pupil_area = tel.pupil_area * (100.)**2 # m^2 to cm^2 + pupil_area = tel.pupil_area * (100.) ** 2 # m^2 to cm^2 return flux * pupil_area * exptime def getElectronFluxFilt(self, filt, tel, exptime=150.): @@ -73,27 +75,42 @@ class MockObject(object): def getPosWorld(self): ra = self.param["ra"] 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): + def getPosImg_Offset_WCS(self, img, fdmodel=None, chip=None, verbose=True, img_header=None): self.posImg = img.wcs.toImage(self.getPosWorld()) self.localWCS = img.wcs.local(self.posImg) 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) + print("x = %.2f, y = %.2f\n" % (self.posImg.x, self.posImg.y), flush=True) self.posImg = 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) + 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) - return self.posImg, self.offset, self.localWCS + + from astropy import wcs + + if img_header is not None: + self.real_wcs = galsim.FitsWCS(header=img_header) + else: + self.real_wcs = None + + return self.posImg, self.offset, self.localWCS, self.real_wcs + + 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 drawObject(self, img, final, flux=None, filt=None, tel=None, exptime=150.): """ Draw (point like) object on img. @@ -112,7 +129,7 @@ class MockObject(object): # Draw with Photon Shoot stamp = final.drawImage(wcs=self.localWCS, method='phot', offset=self.offset) - + stamp.setCenter(self.x_nominal, self.y_nominal) if np.sum(np.isnan(stamp.array)) >= 1: stamp.setZero() @@ -123,7 +140,8 @@ class MockObject(object): img[bounds] += stamp[bounds] return img, stamp, isUpdated - def drawObj_multiband(self, tel, pos_img, psf_model, bandpass_list, filt, chip, nphotons_tot=None, g1=0, g2=0, exptime=150.): + def drawObj_multiband(self, tel, pos_img, psf_model, bandpass_list, filt, chip, nphotons_tot=None, g1=0, g2=0, + exptime=150.): if nphotons_tot == None: nphotons_tot = self.getElectronFluxFilt(filt, tel, exptime) # print("nphotons_tot = ", nphotons_tot) @@ -146,6 +164,18 @@ class MockObject(object): 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.real_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) + + real_wcs_local = self.real_wcs.local(self.real_pos) + for i in range(len(bandpass_list)): bandpass = bandpass_list[i] try: @@ -155,8 +185,8 @@ class MockObject(object): self.logger.error(e) # return False continue - - ratio = sub/full + + ratio = sub / full if not (ratio == -1 or (ratio != ratio)): nphotons = ratio * nphotons_tot @@ -165,26 +195,33 @@ class MockObject(object): continue nphotons_sum += nphotons # print("nphotons_sub-band_%d = %.2f"%(i, nphotons)) - psf, pos_shear = psf_model.get_PSF(chip=chip, pos_img=pos_img, bandpass=bandpass, folding_threshold=folding_threshold) + 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) - stamp = star.drawImage(wcs=self.localWCS, method='phot', offset=self.offset, save_photons=True) + # stamp = star.drawImage(wcs=self.localWCS, method='phot', offset=self.offset, save_photons=True) + # xmax = max(xmax, stamp.xmax) + # ymax = max(ymax, stamp.ymax) + # photons = stamp.photons + # photons.x += self.x_nominal + # photons.y += self.y_nominal + # photons_list.append(photons) + + stamp = star.drawImage(wcs=real_wcs_local, method='phot', offset=offset, save_photons=True) xmax = max(xmax, stamp.xmax) ymax = max(ymax, stamp.ymax) photons = stamp.photons - photons.x += self.x_nominal - photons.y += self.y_nominal + photons.x += x_nominal + photons.y += y_nominal photons_list.append(photons) - # Test stamp size - # print(xmax, ymax) - - stamp = galsim.ImageF(int(xmax*1.1), int(ymax*1.1)) - stamp.wcs = self.localWCS - stamp.setCenter(self.x_nominal, self.y_nominal) - bounds = stamp.bounds & chip.img.bounds + stamp = galsim.ImageF(int(xmax * 1.1), int(ymax * 1.1)) + stamp.wcs = real_wcs_local + stamp.setCenter(x_nominal, y_nominal) + bounds = stamp.bounds & galsim.BoundsI(0, chip.npix_x - 1, 0, chip.npix_y - 1) + chip.img.setOrigin(0, 0) stamp[bounds] = chip.img[bounds] for i in range(len(photons_list)): if i == 0: @@ -193,20 +230,55 @@ class MockObject(object): chip.sensor.accumulate(photons_list[i], stamp, resume=True) chip.img[bounds] = stamp[bounds] + + chip.img.setOrigin(chip.bound.xmin, chip.bound.ymin) + # Test stamp size + # print(xmax, ymax) + + # stamp = galsim.ImageF(int(xmax*1.1), int(ymax*1.1)) + # stamp.wcs = self.localWCS + # stamp.setCenter(self.x_nominal, self.y_nominal) + # bounds = stamp.bounds & chip.img.bounds + # stamp[bounds] = chip.img[bounds] + # for i in range(len(photons_list)): + # if i == 0: + # chip.sensor.accumulate(photons_list[i], stamp) + # else: + # chip.sensor.accumulate(photons_list[i], stamp, resume=True) + # + # chip.img[bounds] = stamp[bounds] # print(chip.img.array.sum()) # print("nphotons_sum = ", nphotons_sum) del photons_list del stamp return True, pos_shear - - def addSLStoChipImage(self, sdp = None, chip = None, xOrderSigPlus = None): + 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 @@ -215,21 +287,23 @@ class MockObject(object): 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 = self.localWCS + stamp.wcs = local_wcs stamp.setOrigin(origin_order_x, origin_order_y) - bounds = stamp.bounds & chip.img.bounds + 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 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): - + norm_thr_rang_ids = normFilter['SENSITIVITY'] > 0.001 sedNormFactor = getNormFactorForSpecWithABMAG(ABMag=self.param['mag_use_normal'], spectrum=self.sed, norm_thr=normFilter, @@ -249,19 +323,34 @@ class MockObject(object): 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.real_wcs) - xOrderSigPlus = {'A':1.3909419820029296,'B':1.4760376591236062,'C':4.035447379743442,'D':5.5684364343742825,'E':16.260021029735388} - grating_split_pos_chip = chip.bound.xmin + grating_split_pos + 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) + + real_wcs_local = self.real_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 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) + 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) star = galsim.Convolve(psf, star) - starImg = star.drawImage(nx=100, ny=100, wcs=self.localWCS) + starImg = star.drawImage(nx=100, ny=100, wcs=real_wcs_local) - origin_star = [self.y_nominal - (starImg.center.y - starImg.ymin), - self.x_nominal - (starImg.center.x - starImg.xmin)] + origin_star = [y_nominal - (starImg.center.y - starImg.ymin), + x_nominal - (starImg.center.x - starImg.xmin)] 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] @@ -273,52 +362,56 @@ class MockObject(object): subImg_p1 = starImg.array[:, 0:subSlitPos] star_p1 = galsim.Image(subImg_p1) origin_p1 = origin_star - xcenter_p1 = min(self.x_nominal,grating_split_pos_chip-1) - chip.bound.xmin - ycenter_p1 = self.y_nominal-chip.bound.ymin + 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=bandpass.blue_limit * 10, band_end=bandpass.red_limit * 10, - conf=chip.sls_conf[0], - isAlongY=0) + ycenter=ycenter_p1, origin=origin_p1, + tar_spec=normalSED, + band_start=bandpass.blue_limit * 10, band_end=bandpass.red_limit * 10, + conf=chip.sls_conf[0], + isAlongY=0, + flat_cube=flat_cube) - self.addSLStoChipImage(sdp=sdp_p1, chip=chip, xOrderSigPlus = xOrderSigPlus) + self.addSLStoChipImage(sdp=sdp_p1, chip=chip, xOrderSigPlus=xOrderSigPlus, local_wcs=real_wcs_local) - subImg_p2 = starImg.array[:, subSlitPos+1:starImg.array.shape[1]] + subImg_p2 = starImg.array[:, subSlitPos + 1:starImg.array.shape[1]] star_p2 = galsim.Image(subImg_p2) origin_p2 = [origin_star[0], grating_split_pos_chip] - xcenter_p2 = max(self.x_nominal, grating_split_pos_chip - 1) - chip.bound.xmin - ycenter_p2 = self.y_nominal - chip.bound.ymin + 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=bandpass.blue_limit * 10, band_end=bandpass.red_limit * 10, conf=chip.sls_conf[1], - isAlongY=0) + isAlongY=0, + flat_cube=flat_cube) - self.addSLStoChipImage(sdp=sdp_p2, chip=chip,xOrderSigPlus = xOrderSigPlus) + self.addSLStoChipImage(sdp=sdp_p2, chip=chip, xOrderSigPlus=xOrderSigPlus, local_wcs=real_wcs_local) del sdp_p1 del sdp_p2 - elif grating_split_pos_chip<=gal_origin[1]: - sdp = SpecDisperser(orig_img=starImg, xcenter=self.x_nominal - chip.bound.xmin, - ycenter=self.y_nominal - chip.bound.ymin, origin=origin_star, + 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=bandpass.blue_limit * 10, band_end=bandpass.red_limit * 10, conf=chip.sls_conf[1], - isAlongY=0) - self.addSLStoChipImage(sdp=sdp, chip=chip, xOrderSigPlus = xOrderSigPlus) + isAlongY=0, + flat_cube=flat_cube) + self.addSLStoChipImage(sdp=sdp, chip=chip, xOrderSigPlus=xOrderSigPlus, local_wcs=real_wcs_local) del sdp - elif grating_split_pos_chip>=gal_end[1]: - sdp = SpecDisperser(orig_img=starImg, xcenter=self.x_nominal - chip.bound.xmin, - ycenter=self.y_nominal - chip.bound.ymin, origin=origin_star, + 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=bandpass.blue_limit * 10, band_end=bandpass.red_limit * 10, conf=chip.sls_conf[0], - isAlongY=0) - self.addSLStoChipImage(sdp=sdp, chip=chip, xOrderSigPlus = xOrderSigPlus) + isAlongY=0, + flat_cube=flat_cube) + self.addSLStoChipImage(sdp=sdp, chip=chip, xOrderSigPlus=xOrderSigPlus, local_wcs=real_wcs_local) del sdp del psf return True, pos_shear diff --git a/ObservationSim/MockObject/SkybackgroundMap.py b/ObservationSim/MockObject/SkybackgroundMap.py index 0557090714140efe18384a07a1ecb7c2b61be694..d1da14da333e8099b1f4d14d7504de94698a12be 100644 --- a/ObservationSim/MockObject/SkybackgroundMap.py +++ b/ObservationSim/MockObject/SkybackgroundMap.py @@ -10,6 +10,8 @@ import galsim import os +import time + try: import importlib.resources as pkg_resources except ImportError: @@ -20,7 +22,7 @@ except ImportError: ###calculate sky map by sky SED def calculateSkyMap_split_g(skyMap=None, blueLimit=4200, redLimit=6500, skyfn='sky_emiss_hubble_50_50_A.dat', conf=[''], pixelSize=0.074, isAlongY=0, - split_pos=3685): + split_pos=3685, flat_cube = None): # skyMap = np.ones([yLen, xLen], dtype='float32') # # if isAlongY == 1: @@ -64,33 +66,79 @@ def calculateSkyMap_split_g(skyMap=None, blueLimit=4200, redLimit=6500, skyfn='s # conf=conf1) # sdp.compute_spec_orders() - - sdp = SpecDisperser(orig_img=skyImg1, xcenter=skyImg1.center.x, ycenter=skyImg1.center.y, origin=origin1, - tar_spec=spec, - band_start=tbstart, band_end=tbend, - conf=conf2) - - 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] + 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, ignoreBeam=['D','E']) + + 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] + + # sdp = SpecDisperser(orig_img=skyImg1, xcenter=skyImg1.center.x, ycenter=skyImg1.center.y, origin=origin1, + # tar_spec=spec, + # band_start=tbstart, band_end=tbend, + # conf=conf2, + # flat_cube=flat_cube, ignoreBeam=['D','E']) + # + # 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: - skyImg1 = galsim.Image(skyImg.array[:, 0:split_pos]) - origin1 = [0, 0] - skyImg2 = galsim.Image(skyImg.array[:, split_pos:-1]) - origin2 = [0, split_pos] + # skyImg1 = galsim.Image(skyImg.array[:, 0:split_pos]) + # origin1 = [0, 0] + # skyImg2 = galsim.Image(skyImg.array[:, split_pos:]) + # origin2 = [0, split_pos] # 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, @@ -98,41 +146,101 @@ def calculateSkyMap_split_g(skyMap=None, blueLimit=4200, redLimit=6500, skyfn='s # conf=conf1) # sdp.compute_spec_orders() - - sdp = SpecDisperser(orig_img=skyImg1, xcenter=skyImg1.center.x, ycenter=skyImg1.center.y, origin=origin1, - tar_spec=spec, - band_start=tbstart, band_end=tbend, - conf=conf1) - - 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] - + 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) - sdp = SpecDisperser(orig_img=skyImg2, xcenter=skyImg2.center.x, ycenter=skyImg2.center.y, origin=origin2, - tar_spec=spec, - band_start=tbstart, band_end=tbend, - conf=conf2) - - 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 - fImg[bounds] = fImg[bounds] + ssImg[bounds] + 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 = rotate90(array_orig=fImg.array, xc=0, yc=0, isClockwise=0) diff --git a/ObservationSim/MockObject/SpecDisperser/SpecDisperser.py b/ObservationSim/MockObject/SpecDisperser/SpecDisperser.py index 9acee114d253c3cf559568dc3f2f9ebc13dbdb0d..ffe1ab12edd6d7756bbe8114f4c3ec8c0fb763f5 100644 --- a/ObservationSim/MockObject/SpecDisperser/SpecDisperser.py +++ b/ObservationSim/MockObject/SpecDisperser/SpecDisperser.py @@ -51,7 +51,7 @@ def rotate90(array_orig=None, xc=0, yc=0, isClockwise=0): 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): + 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 @@ -76,9 +76,10 @@ class SpecDisperser(object): 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) + yc=orig_img.center.y, isClockwise=1) self.img_sh = orig_img.array.T.shape self.xcenter = ycenter @@ -94,6 +95,8 @@ class SpecDisperser(object): self.grating_conf = aXeConf(conf) self.grating_conf.get_beams() self.grating_conf_file = conf + self.ignoreBeam = ignoreBeam + def compute_spec_orders(self): @@ -101,7 +104,8 @@ class SpecDisperser(object): beam_names = self.grating_conf.beams for beam in beam_names: - #for beam in ['A','B']: + if beam in self.ignoreBeam: + continue all_orders[beam] = self.compute_spec(beam) return all_orders @@ -125,7 +129,7 @@ class SpecDisperser(object): 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']) @@ -168,21 +172,10 @@ class SpecDisperser(object): nonz = sensitivity_beam != 0 - status = disperse.disperse_grism_object(self.thumb_img, - flat_index[nonz], yfrac_beam[nonz], - sensitivity_beam[nonz], - modelf, x0, - array(self.img_sh, dtype=int64), - array(beam_sh, dtype=int64)) - - model = modelf.reshape(beam_sh) - self.beam_flux[beam] = sum(modelf) - origin_in = zeros_like(self.origin) dx0_in = dx[0] dy0_in = dyc[0] if self.isAlongY == 1: - model, tmx, tmy = rotate90(array_orig=model, isClockwise=0) origin_in[0] = self.origin[0] origin_in[1] = self.origin[1] - len_spec_y dx0_in = -dyc[0] @@ -195,13 +188,73 @@ class SpecDisperser(object): originOut_x = origin_in[1] + dx0_in - 1 originOut_y = origin_in[0] + dy0_in - 1 + 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, + 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 - 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' + 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 = Table(array([w, sens]).T, names=('WAVELENGTH', 'SENSITIVITY')) senstivity_out.write(sens_file_name, format='fits') @@ -209,15 +262,16 @@ class SpecDisperser(object): 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) @@ -237,12 +291,12 @@ class aXeConf(): 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 @@ -316,23 +370,23 @@ class aXeConf(): def field_dependent(self, xi, yi, coeffs): """aXe field-dependent coefficients - + See the `aXe manual `_ 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 @@ -356,34 +410,34 @@ class aXeConf(): 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. + + For `dydx` polynomial orders 0, 1 or 2, integrate analytically. Higher orders must be integrated numerically. - - **Constant:** + + **Constant:** .. math:: dp = dx - **Linear:** + **Linear:** .. math:: dp = \sqrt{1+\mathrm{DYDX}[1]}\cdot dx - - **Quadratic:** + + **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$ ... @@ -432,30 +486,30 @@ class aXeConf(): 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 + 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 + 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 diff --git a/ObservationSim/MockObject/SpecDisperser/disperse_c/disperse.pyx b/ObservationSim/MockObject/SpecDisperser/disperse_c/disperse.pyx index bb728d5f0281cbcc3dbccdf9595256fb07c574dd..f5be0b2794ae24d94c1f17f7e50062bf0f1d26b4 100644 --- a/ObservationSim/MockObject/SpecDisperser/disperse_c/disperse.pyx +++ b/ObservationSim/MockObject/SpecDisperser/disperse_c/disperse.pyx @@ -31,48 +31,90 @@ def disperse_grism_object(np.ndarray[FTYPE_t, ndim=2] flam, 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[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: - ysense: sensitivity use pixel describe + yfrac: + ysens: sensitivity use pixel describe full： output result ，1 dim, y_beam * x_beam - x0: the center of gal in image thumbnail + 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) - - 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]): + + 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 - fl_ij = flam[x0[0]+j, x0[1]+i] #/1.e-17 - if (fl_ij == 0): + 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*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*(1-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 k in range(nk): - k1 = idxl[k]+j*shg[1]+i - if (k1 >= 0) & (k1 < nl): - full[k1] += ysens[k]*fl_ij*yfrac[k] - - k2 = idxl[k]+(j-1)*shg[1]+i - if (k2 >= 0) & (k2 < nl): - full[k2] += ysens[k]*fl_ij*(1-yfrac[k]) - + + 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*yfrac[k] + + k2 = idxl[k]+(j-1)*shg[1]+i + if (k2 >= 0) & (k2 < nl): + full[k2] += ysens[k]*fl_ij*(1-yfrac[k]) + return True @cython.boundscheck(False) diff --git a/ObservationSim/MockObject/_util.py b/ObservationSim/MockObject/_util.py index ffb978b910cc63d1bf281148554af29c54292d51..c6cb573d053f59bb2927f9af7c37f3f34c2b6570 100755 --- a/ObservationSim/MockObject/_util.py +++ b/ObservationSim/MockObject/_util.py @@ -545,10 +545,9 @@ def tag_sed(h5file, model_tag, teff=5000, logg=2, feh=0): flux = np.array(h5file["sed"][path][()]).ravel() return path, wave, flux -from astropy.modeling.models import Gaussian2D -from scipy import signal 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 diff --git a/ObservationSim/ObservationSim.py b/ObservationSim/ObservationSim.py index 4ca0669e0065792472db544b399870cc5095f18d..3d31a0d1e66179ef3513997bc9d3788d7f7fe7b6 100755 --- a/ObservationSim/ObservationSim.py +++ b/ObservationSim/ObservationSim.py @@ -7,6 +7,8 @@ import psutil from astropy.io import fits from datetime import datetime +import traceback + from ObservationSim.Config import config_dir, ChipOutput from ObservationSim.Config.Header import generatePrimaryHeader, generateExtensionHeader from ObservationSim.Instrument import Telescope, Filter, FilterParam, FocalPlane, Chip @@ -77,7 +79,7 @@ class Observation(object): chip_output.logger.info(':::::::::::::::::::::::::::END:::::::::::::::::::::::::::::::::::') if self.config["psf_setting"]["psf_model"] == "Gauss": - psf_model = PSFGauss(chip=chip) + psf_model = PSFGauss(chip=chip, psfRa=self.config["psf_setting"]["psf_rcont"]) elif self.config["psf_setting"]["psf_model"] == "Interp": psf_model = PSFInterp(chip=chip, PSF_data_file=self.path_dict["psf_dir"]) else: @@ -129,6 +131,7 @@ class Observation(object): # elif chip.survey_type == "spectroscopic": # sky_map = calculateSkyMap_split_g(xLen=chip.npix_x, yLen=chip.npix_y, blueLimit=filt.blue_limit, redLimit=filt.red_limit, skyfn=self.path_dict["sky_file"], conf=chip.sls_conf, pixelSize=chip.pix_scale, isAlongY=0) elif chip.survey_type == "spectroscopic": + # chip.loadSLSFLATCUBE(flat_fn='flat_cube.fits') flat_normal = np.ones_like(chip.img.array) if self.config["ins_effects"]["flat_fielding"] == True: # print("SLS flat preprocess,CHIP %d : Creating and applying Flat-Fielding"%chip.chipID, flush=True) @@ -148,12 +151,14 @@ class Observation(object): flat_normal = flat_normal*shuttimg flat_normal = np.array(flat_normal,dtype='float32') sky_map = calculateSkyMap_split_g( - skyMap=flat_normal, - blueLimit=filt.blue_limit, - redLimit=filt.red_limit, - conf=chip.sls_conf, - pixelSize=chip.pix_scale, - isAlongY=0) + skyMap=flat_normal, + blueLimit=filt.blue_limit, + redLimit=filt.red_limit, + conf=chip.sls_conf, + pixelSize=chip.pix_scale, + isAlongY=0, + flat_cube=chip.flat_cube) + # sky_map = np.ones([9216, 9232]) del flat_normal if pointing.pointing_type == 'MS': @@ -207,7 +212,8 @@ class Observation(object): ) except Exception as e: - print(e) + # print(e) + traceback.print_exc() chip_output.logger.error(e) continue @@ -258,8 +264,22 @@ class Observation(object): raise ValueError("Unknown shear input") # chip_output.logger.info("debug point #4") + header_wcs = generateExtensionHeader( + xlen=chip.npix_x, + ylen=chip.npix_y, + ra=ra_cen, + dec=dec_cen, + pa=pointing.img_pa.deg, + gain=chip.gain, + readout=chip.read_noise, + dark=chip.dark_noise, + saturation=90000, + psize=chip.pix_scale, + row_num=chip.rowID, + col_num=chip.colID, + extName='raw') - pos_img, offset, local_wcs = obj.getPosImg_Offset_WCS(img=chip.img, fdmodel=self.fd_model, chip=chip, verbose=False) + pos_img, offset, local_wcs, real_wcs = obj.getPosImg_Offset_WCS(img=chip.img, fdmodel=self.fd_model, chip=chip, verbose=False, img_header=header_wcs) if pos_img.x == -1 or pos_img.y == -1: # Exclude object which is outside the chip area (after field distortion) # print("obj missed!!") @@ -299,8 +319,7 @@ class Observation(object): g1=obj.g1, g2=obj.g2, exptime=pointing.exp_time, - normFilter=norm_filt, - ) + normFilter=norm_filt) # chip_output.logger.info("debug point #7") if isUpdated: # TODO: add up stats @@ -311,7 +330,8 @@ class Observation(object): # print("object omitted", flush=True) continue except Exception as e: - print(e) + # print(e) + traceback.print_exc() chip_output.logger.error(e) pass # Unload SED: diff --git a/ObservationSim/PSF/PSFGauss.py b/ObservationSim/PSF/PSFGauss.py index 508388055bec2e652418f5d41e69685b2d9c0191..234d5d08f8c7e0f2ac7026115cfb185b923ce0d1 100644 --- a/ObservationSim/PSF/PSFGauss.py +++ b/ObservationSim/PSF/PSFGauss.py @@ -6,12 +6,16 @@ from scipy.interpolate import interp1d from ObservationSim.PSF.PSFModel import PSFModel class PSFGauss(PSFModel): - def __init__(self, chip, fwhm=0.187, sigSpin=0., psfRa=0.15): + def __init__(self, chip, fwhm=0.187, sigSpin=0., psfRa=None): self.pix_size = chip.pix_scale self.chip = chip - self.fwhm = fwhm + if psfRa is None: + self.fwhm = fwhm + self.sigGauss = 0.15 + else: + self.fwhm = self.fwhmGauss(r=psfRa) + self.sigGauss = psfRa # 80% light radius self.sigSpin = sigSpin - self.sigGauss = psfRa # 80% light radius self.psf = galsim.Gaussian(flux=1.0,fwhm=fwhm) def perfGauss(self, r, sig): @@ -51,9 +55,8 @@ class PSFGauss(PSFModel): gaussImg = gaussImg.array size = np.size(gaussImg,axis=0) cxy = 0.5*(size-1) - flux, ferr, flag = sep.sum_circle(gaussImg,cxy,cxy,r/pscale,subpix=0) - frac = flux.tolist() - return frac + flux, ferr, flag = sep.sum_circle(gaussImg,[cxy],[cxy],[r/pscale],subpix=0) + return flux def fwhmGauss(self, r=0.15,fr=0.8,pscale=None): """ diff --git a/ObservationSim/PSF/PSFInterp.py b/ObservationSim/PSF/PSFInterp.py index 058acca032f174f0dcaf36d75f7f146937e7e29d..45bf211c8d648f9675cb6ab68c3b2f92ec6cef43 100644 --- a/ObservationSim/PSF/PSFInterp.py +++ b/ObservationSim/PSF/PSFInterp.py @@ -350,12 +350,21 @@ class PSFInterp(PSFModel): if galsimGSObject: imPSFt = np.zeros([257,257]) imPSFt[0:256, 0:256] = imPSF + # imPSFt[120:130, 0:256] = 1. img = galsim.ImageF(imPSFt, scale=pixSize) gsp = galsim.GSParams(folding_threshold=folding_threshold) - self.psf = galsim.InterpolatedImage(img, gsparams=gsp).rotate(pointing_pa*galsim.degrees) - - return self.PSFspin(x=px/0.01, y=py/0.01) + ############TEST: START + # Use sheared PSF to test the PSF orientation + # self.psf = galsim.InterpolatedImage(img, gsparams=gsp).shear(g1=0.8, g2=0.) + ############TEST: END + self.psf = galsim.InterpolatedImage(img, gsparams=gsp) + wcs = chip.img.wcs.local(pos_img) + scale = galsim.PixelScale(0.074) + self.psf = wcs.toWorld(scale.toImage(self.psf), image_pos=(pos_img)) + + # return self.PSFspin(x=px/0.01, y=py/0.01) + return self.psf, galsim.Shear(e=0., beta=(np.pi/2)*galsim.radians) return imPSF def PSFspin(self, x, y): diff --git a/config/config_NGP_dev.yaml b/config/config_NGP_dev.yaml index 27bcee0a2c8e72d34d5251badcfee2aede2694b0..9f5ffdbaedf48c0ba4eba4e25458b3d7aa8d7643 100644 --- a/config/config_NGP_dev.yaml +++ b/config/config_NGP_dev.yaml @@ -10,9 +10,9 @@ # Can add some of the command-line arguments here as well; # OK to pass either way or both, as long as they are consistent # work_dir: "/public/home/fangyuedong/sim_code_release/CSST/test/" -work_dir: "/share/home/fangyuedong/csst-simulation/workplace/" +work_dir: "/share/home/fangyuedong/test_psf_rot/csst-simulation/workplace/" data_dir: "/share/simudata/CSSOSDataProductsSims/data/" -run_name: "header_TEST" +run_name: "PSF_TEST" # (Optional) a file of point list # if you just want to run default pointing: @@ -34,7 +34,7 @@ run_option: out_cat_only: NO # Only simulate stars? - star_only: NO + star_only: YES # Only simulate galaxies? galaxy_only: NO @@ -48,7 +48,7 @@ obs_setting: # "Photometric": simulate photometric chips only # "Spectroscopic": simulate slitless spectroscopic chips only # "All": simulate full focal plane - survey_type: "Photometric" + survey_type: "All" # Exposure time [seconds] exp_time: 150. diff --git a/run_NGP.pbs b/run_NGP.pbs index 2763d98eaa45dda2f11a481db4561acd226d1c7f..f4d181ea980833987d570ea233891327598c42a3 100755 --- a/run_NGP.pbs +++ b/run_NGP.pbs @@ -1,7 +1,7 @@ #!/bin/bash #PBS -N SIMS -#PBS -l nodes=wcl-1:ppn=32+wcl-6:ppn=32 +#PBS -l nodes=wcl-1:ppn=8+wcl-2:ppn=8+wcl-3:ppn=8+wcl-4:ppn=8+wcl-5:ppn=8+wcl-6:ppn=8 #PBS -u fangyuedong ###PBS -j oe @@ -12,5 +12,5 @@ date echo $NP # mpirun -np $NP python /public/home/fangyuedong/test/CSST/run_sim.py config_NGP.yaml -c /public/home/fangyuedong/test/CSST/config -mpirun -np $NP python3 /share/home/fangyuedong/csst-simulation/run_sim.py config_NGP_dev.yaml -c /share/home/fangyuedong/csst-simulation/config +mpirun -np $NP python3 /share/home/fangyuedong/test_psf_rot/csst-simulation/run_sim.py config_NGP_dev.yaml -c /share/home/fangyuedong/test_psf_rot/csst-simulation/config diff --git a/setup.py b/setup.py index 7db1250b91c1c868f3c79bf139fdf0cd55e02df3..21918a34d95c39c58ee91a54cc35e2e0d1fe96cb 100644 --- a/setup.py +++ b/setup.py @@ -29,17 +29,17 @@ setup(name='CSSTSim', version='1.0.4', packages=find_packages(), install_requires=[ - 'numpy>=1.18.5', - 'galsim>=2.2.4', - 'pyyaml>=5.3.1', - 'astropy>=4.0.1', - 'scipy>=1.5.0', - 'mpi4py>=3.0.3', - 'sep>=1.0.3', - 'healpy>=1.14.0', - 'h5py>=2.10.0', - 'Cython>=0.29.21', - 'numba>=0.50.1' + # 'numpy>=1.18.5', + # 'galsim>=2.2.4', + # 'pyyaml>=5.3.1', + # 'astropy>=4.0.1', + # 'scipy>=1.5.0', + # 'mpi4py>=3.0.3', + # 'sep>=1.0.3', + # 'healpy>=1.14.0', + # 'h5py>=2.10.0', + # 'Cython>=0.29.21', + # 'numba>=0.50.1' ], package_data = { 'ObservationSim.Astrometry.lib': ['libshao.so'], @@ -49,7 +49,8 @@ setup(name='CSSTSim', 'ObservationSim.Instrument.data.filters': ['*.txt', '*.list', '*.dat'], 'ObservationSim.Instrument.data.throughputs': ['*.txt', '*.dat'], 'ObservationSim.Instrument.data.sls_conf': ['*.conf', '*.fits'], + 'ObservationSim.Instrument.data.flatCube': ['*.fits'], 'Catalog.data': ['*.fits'], }, ext_modules = cythonize(extensions), -) \ No newline at end of file +)