Commit 1b917613 authored by xin's avatar xin
Browse files

updata galaxy, mockobject wcs

parent cd1351bc
Loading
Loading
Loading
Loading
+5 −0
Original line number Diff line number Diff line
@@ -70,6 +70,11 @@ class ChipOutput(object):
    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
        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, 
+85 −25
Original line number Diff line number Diff line
@@ -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.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
@@ -236,10 +296,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,8 +311,8 @@ 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,
@@ -261,13 +321,13 @@ class Galaxy(MockObject):
                                    conf=chip.sls_conf[0],
                                    isAlongY=0)

                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,
@@ -276,27 +336,27 @@ class Galaxy(MockObject):
                                       conf=chip.sls_conf[1],
                                       isAlongY=0)

                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)
                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)
                self.addSLStoChipImage(sdp=sdp, chip=chip, xOrderSigPlus = xOrderSigPlus, local_wcs=real_wcs_local)
                del sdp

            # print(self.y_nominal, starImg.center.y, starImg.ymin)
+35 −22
Original line number Diff line number Diff line
@@ -162,16 +162,16 @@ class MockObject(object):
            folding_threshold = 5.e-3
        gsp = galsim.GSParams(folding_threshold=folding_threshold)

        real_pos = self.getRealPos(chip.img, global_x=self.posImg.x, global_y=self.posImg.y, img_real_wcs=self.real_wcs)
        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 = real_pos.x + 0.5, real_pos.y + 0.5
        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(real_pos)
        real_wcs_local = self.real_wcs.local(self.real_pos)

        for i in range(len(bandpass_list)):
            bandpass = bandpass_list[i]
@@ -252,7 +252,7 @@ class MockObject(object):
        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]
@@ -267,15 +267,17 @@ 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

@@ -301,19 +303,30 @@ 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)

        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)


        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]
            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]
@@ -325,8 +338,8 @@ 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,
@@ -335,13 +348,13 @@ class MockObject(object):
                                    conf=chip.sls_conf[0],
                                    isAlongY=0)

                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,
@@ -350,27 +363,27 @@ class MockObject(object):
                                       conf=chip.sls_conf[1],
                                       isAlongY=0)

                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)
                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)
                self.addSLStoChipImage(sdp=sdp, chip=chip, xOrderSigPlus = xOrderSigPlus, local_wcs=real_wcs_local)
                del sdp
            del psf
        return True, pos_shear