fix skymap sls center pos

ObservationSim/MockObject/SkybackgroundMap.py

@@ -58,7 +58,6 @@ def calculateSkyMap_split_g(skyMap=None, blueLimit=4200, redLimit=6500, skyfn='s
     if split_pos >= skyImg.array.shape[directParm]:
         skyImg1 = galsim.Image(skyImg.array)
         origin1 = [0, 0]
-        skyImg1.setOrigin(origin1)
         # 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,
         #                                   origin=origin1,
@@ -91,10 +90,8 @@ def calculateSkyMap_split_g(skyMap=None, blueLimit=4200, redLimit=6500, skyfn='s
 
         skyImg1 = galsim.Image(skyImg.array[:, 0:split_pos])
         origin1 = [0, 0]
-        skyImg1.setOrigin(origin1)
         skyImg2 = galsim.Image(skyImg.array[:, split_pos:-1])
         origin2 = [0, split_pos]
-        skyImg2.setOrigin(origin2)
 
         # 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,
@@ -122,8 +119,7 @@ def calculateSkyMap_split_g(skyMap=None, blueLimit=4200, redLimit=6500, skyfn='s
                 continue
             fImg[bounds] = fImg[bounds] + ssImg[bounds]
 
-        
-        sdp = SpecDisperser(orig_img=skyImg2, xcenter=skyImg2.center.x, ycenter=skyImg2.center.y, origin=origin2,
+        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,

ObservationSim/MockObject/SpecDisperser/SpecDisperser.py

@@ -190,16 +190,46 @@ class SpecDisperser(object):
             beam_flat = None
         else:
             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]
+
+            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],

ObservationSim/MockObject/SpecDisperser/disperse_c/disperse.pyx

@@ -61,7 +61,8 @@ def disperse_grism_object(np.ndarray[FTYPE_t, ndim=2] flam,
             if(len(flat[0])==4):
                 doflat=1
 
-
+    lambd_2 = wlambda*wlambda
+    lambd_3 = wlambda*wlambda*wlambda
 
     for i in range(0-x0[1], x0[1]):
         if (x0[1]+i < 0) | (x0[1]+i >= shd[1]):
@@ -81,8 +82,7 @@ def disperse_grism_object(np.ndarray[FTYPE_t, ndim=2] flam,
                     flux_factor = 1.
                     if (doflat==1):
                         a = flat[k1]
-                        lamb = wlambda[k]
-                        flux_factor = a[0] + a[1]*lamb + a[2]*lamb*lamb + a[3]*lamb*lamb
+                        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
 
                 k2 = idxl[k]+(j-1)*shg[1]+i
@@ -90,8 +90,7 @@ def disperse_grism_object(np.ndarray[FTYPE_t, ndim=2] flam,
                     flux_factor = 1.
                     if (doflat==1):
                         a = flat[k2]
-                       lamb = wlambda[k]
-                       flux_factor = a[0] + a[1]*lamb + a[2]*lamb*lamb + a[3]*lamb*lamb
+                        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
     
     return True