optimize skylight disperse;fig bug gauss psf error for seg lib

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)

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

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
 

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):
 		"""