close flat cube

ObservationSim/MockObject/SkybackgroundMap.py

@@ -97,7 +97,7 @@ def calculateSkyMap_split_g(skyMap=None, blueLimit=4200, redLimit=6500, skyfn='s
                                 tar_spec=spec,
                                 band_start=tbstart, band_end=tbend,
                                 conf=conf2,
-                                flat_cube=None, ignoreBeam=['D','E'])
+                                flat_cube=flat_cube, ignoreBeam=['D','E'])
 
                 spec_orders = sdp.compute_spec_orders()
 
@@ -169,7 +169,7 @@ def calculateSkyMap_split_g(skyMap=None, blueLimit=4200, redLimit=6500, skyfn='s
             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)
+                # 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]
@@ -179,7 +179,7 @@ def calculateSkyMap_split_g(skyMap=None, blueLimit=4200, redLimit=6500, skyfn='s
                                 tar_spec=spec,
                                 band_start=tbstart, band_end=tbend,
                                 conf=conf1,
-                                flat_cube=None, ignoreBeam=['B','C','D','E'])
+                                flat_cube=flat_cube)
 
                 spec_orders = sdp.compute_spec_orders()
 
@@ -212,7 +212,7 @@ def calculateSkyMap_split_g(skyMap=None, blueLimit=4200, redLimit=6500, skyfn='s
             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)
+                # print(i,j,sub_y_s, sub_y_e,sub_x_s,sub_x_e)
                 
                 T1 = time.time()
 
@@ -224,7 +224,7 @@ def calculateSkyMap_split_g(skyMap=None, blueLimit=4200, redLimit=6500, skyfn='s
                                 tar_spec=spec,
                                 band_start=tbstart, band_end=tbend,
                                 conf=conf2,
-                                flat_cube=None, ignoreBeam=['B','C','D','E'])
+                                flat_cube=flat_cube)
 
                 spec_orders = sdp.compute_spec_orders()
 

ObservationSim/MockObject/SpecDisperser/disperse_c/disperse.pyx

@@ -53,18 +53,23 @@ def disperse_grism_object(np.ndarray[FTYPE_t, ndim=2] flam,
 
     nk = len(idxl)
     nl = len(full)
+    yfrac_off = 1-yfrac
 
     if (flat is not None):
         nlamb = len(wlambda)
+        nflat = len(flat)
         flat_lamb_min = 2500
         flat_lamb_max = 10000
-        flat_eff_all = np.ones([nl, nlamb])*ysens
         lambda_co = np.zeros([len(flat[0]),nlamb])
-        lambda_co[0] = lambda_co[0] + 1
-        lambda_co[1] = (wlambda-flat_lamb_min)/(flat_lamb_max-flat_lamb_min)
-        lambda_co[2] = lambda_co[1]*lambda_co[1]
-        lambda_co[3] = lambda_co[1]*lambda_co[2]
-        flat_eff_all = np.dot(flat, lambda_co)*ysens
+        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]):
@@ -81,11 +86,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):
-                        full[k1] += fl_ij*yfrac[k]*flat_eff_all[k1,k]
+                        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):
-                        full[k2] += fl_ij*(1-yfrac[k])*flat_eff_all[k2,k]
+                        flat_ids = k2*nlamb+k
+                        full[k2] += fl_ij*yfrac_off[k]*flat_eff_all[flat_ids]
 
     else:
         for i in range(0-x0[1], x0[1]):
@@ -107,7 +114,7 @@ def disperse_grism_object(np.ndarray[FTYPE_t, ndim=2] flam,
 
                     k2 = idxl[k]+(j-1)*shg[1]+i
                     if (k2 >= 0) & (k2 < nl):
-                        full[k2] += fl_ij*(1-yfrac[k])*ysens[k]
+                        full[k2] += fl_ij*yfrac_off[k]*ysens[k]
 
     return True
 

ObservationSim/ObservationSim.py

@@ -131,7 +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')
+            # 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)
@@ -150,7 +150,6 @@ class Observation(object):
                                                     dt=1E-3)  # shutter effect normalized image for this chip
                 flat_normal = flat_normal*shuttimg
                 flat_normal = np.array(flat_normal,dtype='float32')
-            print('DEBUG: sky_map_flat-----------------------')
             sky_map = calculateSkyMap_split_g(
                 skyMap=flat_normal,
                 blueLimit=filt.blue_limit,
@@ -159,7 +158,6 @@ class Observation(object):
                 pixelSize=chip.pix_scale,
                 isAlongY=0,
                 flat_cube=chip.flat_cube)
-            print('DEBUG: sky_map_flat-----------------------')
             # sky_map = np.ones([9216, 9232])
             del flat_normal