diff --git a/ObservationSim/MockObject/SkybackgroundMap.py b/ObservationSim/MockObject/SkybackgroundMap.py index 4a6ce90e0d741090972fc5329b1709a691d71e61..b63431b249428994113b8b50d529eb5837fae5c1 100644 --- a/ObservationSim/MockObject/SkybackgroundMap.py +++ b/ObservationSim/MockObject/SkybackgroundMap.py @@ -64,34 +64,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, - 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] + 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, @@ -99,42 +144,88 @@ 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, - 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] - - sdp = SpecDisperser(orig_img=skyImg2, xcenter=skyImg2.center.x+split_pos, ycenter=skyImg2.center.y, origin=origin2, - 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 - fImg[bounds] = fImg[bounds] + ssImg[bounds] + y_len = skyMap.shape[0] + x_len = skyMap.shape[1] + delt_x = 500 + delt_y = 500 + + 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, 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] + + 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, 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] + + 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] + + 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] 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 1b776da87d6cfa4d76eb36a48773a26082efe9e6..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, flat_cube=None): + 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 @@ -95,6 +95,7 @@ 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): @@ -103,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 diff --git a/ObservationSim/MockObject/SpecDisperser/disperse_c/disperse.pyx b/ObservationSim/MockObject/SpecDisperser/disperse_c/disperse.pyx index a96756ad841c0c395a3b8b7f8698dff1ccef6ecc..be1faffda2bd6ccfbbbdc1c035e28b32e4eab3b6 100644 --- a/ObservationSim/MockObject/SpecDisperser/disperse_c/disperse.pyx +++ b/ObservationSim/MockObject/SpecDisperser/disperse_c/disperse.pyx @@ -41,7 +41,7 @@ def disperse_grism_object(np.ndarray[FTYPE_t, ndim=2] flam, 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 + ysens: sensitivity use pixel describe full： output result ，1 dim, y_beam * x_beam x0: the center of gal in image thumbnail (y, x) shd: shape of direct image @@ -53,16 +53,20 @@ def disperse_grism_object(np.ndarray[FTYPE_t, ndim=2] flam, nk = len(idxl) nl = len(full) - doflat = 0 + + nlamb = len(wlambda) + flat_eff_all = np.ones([nl, nlamb])*ysens + if (flat is not None): nflat = len(flat) - print('DEBUG: nflat', nflat,nl,len(flat[0])) if(nflat == nl): if(len(flat[0])==4): - doflat=1 - - lambd_2 = wlambda*wlambda - lambd_3 = wlambda*wlambda*wlambda + lambda_co = np.zeros([4,nlamb]) + lambda_co[0] = lambda_co[0] + 1 + lambda_co[1] = wlambda*1 + lambda_co[2] = wlambda*lambda_co[1] + lambda_co[3] = wlambda*lambda_co[2] + flat_eff_all = np.dot(flat, lambda_co)*ysens for i in range(0-x0[1], x0[1]): if (x0[1]+i < 0) | (x0[1]+i >= shd[1]): @@ -79,19 +83,13 @@ def disperse_grism_object(np.ndarray[FTYPE_t, ndim=2] flam, for k in range(nk): k1 = idxl[k]+j*shg[1]+i if (k1 >= 0) & (k1 < nl): - flux_factor = 1. - if (doflat==1): - a = flat[k1] - flux_factor = a[0] + a[1]*wlambda[k] + a[2]*lambd_2[k] + a[3]*lambd_3[k] - full[k1] += ysens[k]*fl_ij*yfrac[k]*flux_factor + flux_factor = flat_eff_all[k1,k] + full[k1] += fl_ij*yfrac[k]*flux_factor k2 = idxl[k]+(j-1)*shg[1]+i if (k2 >= 0) & (k2 < nl): - flux_factor = 1. - if (doflat==1): - a = flat[k2] - flux_factor = a[0] + a[1]*wlambda[k] + a[2]*lambd_2[k] + a[3]*lambd_3[k] - full[k2] += ysens[k]*fl_ij*(1-yfrac[k])*flux_factor + flux_factor = flat_eff_all[k2,k] + full[k2] += fl_ij*(1-yfrac[k])*flux_factor return True diff --git a/ObservationSim/PSF/PSFGauss.py b/ObservationSim/PSF/PSFGauss.py index 7d09c3c67096e318d3935189a4e28b91c6b9deee..234d5d08f8c7e0f2ac7026115cfb185b923ce0d1 100644 --- a/ObservationSim/PSF/PSFGauss.py +++ b/ObservationSim/PSF/PSFGauss.py @@ -55,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): """