Galaxy.py 18.8 KB
Newer Older
Fang Yuedong's avatar
Fang Yuedong committed
1
2
3
4
5
6
7
import numpy as np
import galsim
import os, sys
import astropy.constants as cons
from astropy.table import Table
from scipy import interpolate

8
9
10
11
from ObservationSim.MockObject._util import eObs, integrate_sed_bandpass, getNormFactorForSpecWithABMAG, getObservedSED, getABMAG,convolveGaussXorders
from ObservationSim.MockObject.SpecDisperser import SpecDisperser
from ObservationSim.MockObject.MockObject import MockObject

Fang Yuedong's avatar
Fang Yuedong committed
12
13
# import tracemalloc

Fang Yuedong's avatar
Fang Yuedong committed
14
class Galaxy(MockObject):
15
    def __init__(self, param, logger=None):
16
        super().__init__(param, logger=logger)
Fang Yuedong's avatar
Fang Yuedong committed
17

Fang Yuedong's avatar
Fang Yuedong committed
18
19
20
21
        if not hasattr(self, "disk_sersic_idx"):
            self.disk_sersic_idx = 1.
        if not hasattr(self, "bulge_sersic_idx"):
            self.bulge_sersic_idx = 4.
Fang Yuedong's avatar
Fang Yuedong committed
22
23
24
25
26
27

    def unload_SED(self):
        """(Test) free up SED memory
        """
        del self.sed

Fang Yuedong's avatar
Fang Yuedong committed
28
    def getGSObj_multiband(self, tel, psf_list, bandpass_list, filt, nphotons_tot=None, g1=0, g2=0, exptime=150., fd_shear=None):
Fang Yuedong's avatar
Fang Yuedong committed
29
30
31
32
33
34
35
36
37
38
39
        if len(psf_list) != len(bandpass_list):
            raise ValueError("!!!The number of PSF profiles and the number of bandpasses must be equal.")
        objs = []
        if nphotons_tot == None:
            nphotons_tot = self.getElectronFluxFilt(filt, tel, exptime)
        # print("nphotons_tot = ", nphotons_tot)

        try:
            full = integrate_sed_bandpass(sed=self.sed, bandpass=filt.bandpass_full)
        except Exception as e:
            print(e)
40
41
            if self.logger:
                self.logger.error(e)
Fang Yuedong's avatar
Fang Yuedong committed
42
43
44
45
46
47
48
            return -1
        for i in range(len(bandpass_list)):
            bandpass = bandpass_list[i]
            try:
                sub = integrate_sed_bandpass(sed=self.sed, bandpass=bandpass)
            except Exception as e:
                print(e)
49
50
                if self.logger:
                    self.logger.error(e)
Fang Yuedong's avatar
Fang Yuedong committed
51
52
53
54
55
56
57
58
59
                return -1
            
            ratio = sub/full
            if not (ratio == -1 or (ratio != ratio)):
                nphotons = ratio * nphotons_tot
            else:
                return -1

            psf = psf_list[i]
Fang Yuedong's avatar
Fang Yuedong committed
60
            disk = galsim.Sersic(n=self.disk_sersic_idx, half_light_radius=self.hlr_disk, flux=1.0)
Fang Yuedong's avatar
Fang Yuedong committed
61
62
            disk_shape = galsim.Shear(g1=self.e1_disk, g2=self.e2_disk)
            disk = disk.shear(disk_shape)
Fang Yuedong's avatar
Fang Yuedong committed
63
            bulge = galsim.Sersic(n=self.bulge_sersic_idx, half_light_radius=self.hlr_bulge, flux=1.0)
Fang Yuedong's avatar
Fang Yuedong committed
64
65
66
67
68
            bulge_shape = galsim.Shear(g1=self.e1_bulge, g2=self.e2_bulge)
            bulge = bulge.shear(bulge_shape)

            gal = self.bfrac * bulge + (1.0 - self.bfrac) * disk
            gal = gal.withFlux(nphotons)
69
70
71
            if fd_shear is not None:
                g1 += fd_shear.g1
                g2 += fd_shear.g2
Fang Yuedong's avatar
Fang Yuedong committed
72
73
74
            gal_shear = galsim.Shear(g1=g1, g2=g2)
            gal = gal.shear(gal_shear)
            gal = galsim.Convolve(psf, gal)
Fang Yuedong's avatar
Fang Yuedong committed
75
            
Fang Yuedong's avatar
Fang Yuedong committed
76
77
78
79
            objs.append(gal)
        final = galsim.Sum(objs)
        return final

Fang Yuedong's avatar
Fang Yuedong committed
80
    def drawObj_multiband(self, tel, pos_img, psf_model, bandpass_list, filt, chip, nphotons_tot=None, g1=0, g2=0, exptime=150., fd_shear=None):
Fang Yuedong's avatar
Fang Yuedong committed
81
82
83
84
85
86
87
88
        if nphotons_tot == None:
            nphotons_tot = self.getElectronFluxFilt(filt, tel, exptime)
        # print("nphotons_tot = ", nphotons_tot)

        try:
            full = integrate_sed_bandpass(sed=self.sed, bandpass=filt.bandpass_full)
        except Exception as e:
            print(e)
89
90
            if self.logger:
                self.logger.error(e)
Fang Yuedong's avatar
Fang Yuedong committed
91
92
93
            return 2, None
        
        # # [C6 TEST]
Fang Yuedong's avatar
Fang Yuedong committed
94
        # print('hlr_disk = %.4f, hlr_bulge = %.4f'%(self.hlr_disk, self.hlr_bulge))
Fang Yuedong's avatar
Fang Yuedong committed
95
96
        # tracemalloc.start()

Fang Yuedong's avatar
Fang Yuedong committed
97
        big_galaxy = False
Fang Yuedong's avatar
Fang Yuedong committed
98
        if self.hlr_disk > 3.0 or self.hlr_bulge > 3.0: # Very big galaxy
Fang Yuedong's avatar
Fang Yuedong committed
99
100
            big_galaxy = True

101
        # Set Galsim Parameters
Fang Yuedong's avatar
Fang Yuedong committed
102
103
104
105
106
107
        if self.getMagFilter(filt) <= 15 and (not big_galaxy):
            folding_threshold = 5.e-4
        else:
            folding_threshold = 5.e-3
        gsp = galsim.GSParams(folding_threshold=folding_threshold)

xin's avatar
xin committed
108
        self.real_pos = self.getRealPos(chip.img, global_x=self.posImg.x, global_y=self.posImg.y,
109
                                        img_real_wcs=self.chip_wcs)
xin's avatar
xin committed
110
111
112
113
114
115
116

        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)
117
118
119
        # Get real local wcs of object (deal with chip rotation w.r.t its center)
        chip_wcs_local = self.chip_wcs.local(self.real_pos)
        is_updated = 0
xin's avatar
xin committed
120

121
        # Model the galaxy as disk + bulge
122
123
124
125
126
127
        disk = galsim.Sersic(n=self.disk_sersic_idx, half_light_radius=self.hlr_disk, flux=1.0, gsparams=gsp)
        disk_shape = galsim.Shear(g1=self.e1_disk, g2=self.e2_disk)
        disk = disk.shear(disk_shape)
        bulge = galsim.Sersic(n=self.bulge_sersic_idx, half_light_radius=self.hlr_bulge, flux=1.0, gsparams=gsp)
        bulge_shape = galsim.Shear(g1=self.e1_bulge, g2=self.e2_bulge)
        bulge = bulge.shear(bulge_shape)
128
        
129
        # Get shear and deal with shear induced by field distortion
130
131
132
133
        if fd_shear:
            g1 += fd_shear.g1
            g2 += fd_shear.g2
        gal_shear = galsim.Shear(g1=g1, g2=g2)
134

135
        # Loop over all sub-bandpasses
Fang Yuedong's avatar
Fang Yuedong committed
136
137
138
139
140
141
        for i in range(len(bandpass_list)):
            bandpass = bandpass_list[i]
            try:
                sub = integrate_sed_bandpass(sed=self.sed, bandpass=bandpass)
            except Exception as e:
                print(e)
142
143
                if self.logger:
                    self.logger.error(e)
Fang Yuedong's avatar
Fang Yuedong committed
144
145
146
147
148
149
                continue
            ratio = sub/full
            if not (ratio == -1 or (ratio != ratio)):
                nphotons = ratio * nphotons_tot
            else:
                continue
150
151
            
            # nphotons_sum += nphotons
Fang Yuedong's avatar
Fang Yuedong committed
152
            # # [C6 TEST]
Fang Yuedong's avatar
Fang Yuedong committed
153
154
            # print("nphotons_sub-band_%d = %.2f"%(i, nphotons))

155
            # Get PSF model
Fang Yuedong's avatar
Fang Yuedong committed
156
            psf, pos_shear = psf_model.get_PSF(chip=chip, pos_img=pos_img, bandpass=bandpass, folding_threshold=folding_threshold)
157
158
            
            gal_temp = self.bfrac * bulge + (1.0 - self.bfrac) * disk
159
            gal_temp = gal_temp.shear(gal_shear)
Fang Yuedong's avatar
Fang Yuedong committed
160

161
162
163
164
165
166
167
168
            gal_temp = gal_temp.withFlux(nphotons)
            if not big_galaxy: # Not apply PSF for very big galaxy
                gal_temp = galsim.Convolve(psf, gal_temp)
            
            if i == 0:
                gal = gal_temp
            else:
                gal = gal + gal_temp
Fang Yuedong's avatar
Fang Yuedong committed
169

Fang Yuedong's avatar
Fang Yuedong committed
170
            # (TEST) Random knots
171
172
173
            # knots = galsim.RandomKnots(npoints=100, profile=disk)
            # kfrac = np.random.random()*(1.0 - self.bfrac)
            # gal = self.bfrac * bulge + (1.0 - self.bfrac - kfrac) * disk + kfrac * knots
Fang Yuedong's avatar
Fang Yuedong committed
174

Fang Yuedong's avatar
Fang Yuedong committed
175
        # # [C6 TEST]
Fang Yuedong's avatar
Fang Yuedong committed
176
        # print('xmax = %d, ymax = %d '%(xmax, ymax))
Fang Yuedong's avatar
Fang Yuedong committed
177
178
179
180
181
        # # Output memory usage
        # snapshot = tracemalloc.take_snapshot()
        # top_stats = snapshot.statistics('lineno')
        # for stat in top_stats[:10]:
        #     print(stat)
182

183
184
185
186
187
        # stamp = gal.drawImage(wcs=chip_wcs_local, method='phot', offset=offset, save_photons=True)
        stamp = gal.drawImage(wcs=chip_wcs_local, offset=offset)
        if np.sum(np.isnan(stamp.array)) > 0:
            # ERROR happens
            return 2, pos_shear
xin's avatar
xin committed
188
189
        stamp.setCenter(x_nominal, y_nominal)
        bounds = stamp.bounds & galsim.BoundsI(0, chip.npix_x - 1, 0, chip.npix_y - 1)
Fang Yuedong's avatar
Fang Yuedong committed
190
191
        if bounds.area() > 0:
            chip.img.setOrigin(0, 0)
192
193
            chip.img[bounds] += stamp[bounds]
            is_updated = 1
Fang Yuedong's avatar
Fang Yuedong committed
194
            chip.img.setOrigin(chip.bound.xmin, chip.bound.ymin)
195
196
197
            del stamp
        
        if is_updated == 0:
Fang Yuedong's avatar
Fang Yuedong committed
198
199
            # Return code 0: object photons missed this detector
            print("obj %s missed"%(self.id))
200
201
            if self.logger:
                self.logger.info("obj %s missed"%(self.id))
Fang Yuedong's avatar
Fang Yuedong committed
202
203
204
205
206
            return 0, pos_shear
        
        # # [C6 TEST]
        # print("nphotons_sum = ", nphotons_sum)
        return 1, pos_shear
Fang Yuedong's avatar
Fang Yuedong committed
207
208

    def drawObj_slitless(self, tel, pos_img, psf_model, bandpass_list, filt, chip, nphotons_tot=None, g1=0, g2=0,
Fang Yuedong's avatar
Fang Yuedong committed
209
210
211
212
213
214
215
216
217
218
219
                         exptime=150., normFilter=None, grating_split_pos=3685, fd_shear=None):
        if normFilter is not None:
            norm_thr_rang_ids = normFilter['SENSITIVITY'] > 0.001
            sedNormFactor = getNormFactorForSpecWithABMAG(ABMag=self.param['mag_use_normal'], spectrum=self.sed,
                                                        norm_thr=normFilter,
                                                        sWave=np.floor(normFilter[norm_thr_rang_ids][0][0]),
                                                        eWave=np.ceil(normFilter[norm_thr_rang_ids][-1][0]))
            if sedNormFactor == 0:
                return 2, None
        else:
            sedNormFactor = 1.
Fang Yuedong's avatar
Fang Yuedong committed
220
221
222
        normalSED = Table(np.array([self.sed['WAVELENGTH'], self.sed['FLUX'] * sedNormFactor]).T,
                          names=('WAVELENGTH', 'FLUX'))

xin's avatar
xin committed
223
        self.real_pos = self.getRealPos(chip.img, global_x=self.posImg.x, global_y=self.posImg.y,
224
                                        img_real_wcs=self.chip_wcs)
xin's avatar
xin committed
225
226
227
228
229
230
231
232

        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)

233
        chip_wcs_local = self.chip_wcs.local(self.real_pos)
xin's avatar
xin committed
234
235


Fang Yuedong's avatar
Fang Yuedong committed
236
        big_galaxy = False
Fang Yuedong's avatar
Fang Yuedong committed
237
        if self.hlr_disk > 3.0 or self.hlr_bulge > 3.0: # Very big galaxy
Fang Yuedong's avatar
Fang Yuedong committed
238
239
240
241
242
243
244
245
            big_galaxy = True

        if self.getMagFilter(filt) <= 15 and (not big_galaxy):
            folding_threshold = 5.e-4
        else:
            folding_threshold = 5.e-3
        gsp = galsim.GSParams(folding_threshold=folding_threshold)
        # nphotons_sum = 0
246
247
248

        flat_cube = chip.flat_cube

Zhang Xin's avatar
Zhang Xin committed
249
        xOrderSigPlus = {'A':1.3909419820029296,'B':1.4760376591236062,'C':4.035447379743442,'D':5.5684364343742825,'E':16.260021029735388}
xin's avatar
xin committed
250
        grating_split_pos_chip = 0 + grating_split_pos
Zhang Xin's avatar
Zhang Xin committed
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266

        branges = np.zeros([len(bandpass_list), 2])

        # print(hasattr(psf_model, 'bandranges'))

        if hasattr(psf_model, 'bandranges'):
            if psf_model.bandranges is None:
                return 2, None
            if len(psf_model.bandranges) != len(bandpass_list):
                return 2, None
            branges = psf_model.bandranges
        else:
            for i in range(len(bandpass_list)):
                branges[i, 0] = bandpass_list[i].blue_limit * 10
                branges[i, 1] = bandpass_list[i].red_limit * 10

Fang Yuedong's avatar
Fang Yuedong committed
267
        for i in range(len(bandpass_list)):
Zhang Xin's avatar
Zhang Xin committed
268
269
            # bandpass = bandpass_list[i]
            brange = branges[i]
Fang Yuedong's avatar
Fang Yuedong committed
270

Zhang Xin's avatar
Zhang Xin committed
271
            # psf, pos_shear = psf_model.get_PSF(chip=chip, pos_img=pos_img, bandpass=bandpass, folding_threshold=folding_threshold)
Fang Yuedong's avatar
Fang Yuedong committed
272
            disk = galsim.Sersic(n=self.disk_sersic_idx, half_light_radius=self.hlr_disk, flux=1.0, gsparams=gsp)
Fang Yuedong's avatar
Fang Yuedong committed
273
274
            disk_shape = galsim.Shear(g1=self.e1_disk, g2=self.e2_disk)
            disk = disk.shear(disk_shape)
Fang Yuedong's avatar
Fang Yuedong committed
275
            bulge = galsim.Sersic(n=self.bulge_sersic_idx, half_light_radius=self.hlr_bulge, flux=1.0, gsparams=gsp)
Fang Yuedong's avatar
Fang Yuedong committed
276
277
278
            bulge_shape = galsim.Shear(g1=self.e1_bulge, g2=self.e2_bulge)
            bulge = bulge.shear(bulge_shape)

279
            gal = self.bfrac * bulge + (1.0 - self.bfrac) * disk
Fang Yuedong's avatar
Fang Yuedong committed
280
281

            # (TEST) Random knots
282
283
284
            # knots = galsim.RandomKnots(npoints=100, profile=disk)
            # kfrac = np.random.random()*(1.0 - self.bfrac)
            # gal = self.bfrac * bulge + (1.0 - self.bfrac - kfrac) * disk + kfrac * knots
Fang Yuedong's avatar
Fang Yuedong committed
285

Fang Yuedong's avatar
Fang Yuedong committed
286
            gal = gal.withFlux(tel.pupil_area * exptime)
287
288
289
            if fd_shear:
                g1 += fd_shear.g1
                g2 += fd_shear.g2
Fang Yuedong's avatar
Fang Yuedong committed
290
291
            gal_shear = galsim.Shear(g1=g1, g2=g2)
            gal = gal.shear(gal_shear)
Zhang Xin's avatar
Zhang Xin committed
292
            # gal = galsim.Convolve(psf, gal)
Fang Yuedong's avatar
Fang Yuedong committed
293

Zhang Xin's avatar
Zhang Xin committed
294
295
296
297
            # if not big_galaxy: # Not apply PSF for very big galaxy
            #     gal = galsim.Convolve(psf, gal)
            #     # if fd_shear is not None:
            #     #     gal = gal.shear(fd_shear)
Fang Yuedong's avatar
Fang Yuedong committed
298

Zhang Xin's avatar
Zhang Xin committed
299
            starImg = gal.drawImage(wcs=chip_wcs_local, offset=offset,method = 'real_space')
Fang Yuedong's avatar
Fang Yuedong committed
300

xin's avatar
xin committed
301
302
            origin_star = [y_nominal - (starImg.center.y - starImg.ymin),
                           x_nominal - (starImg.center.x - starImg.xmin)]
Fang Yuedong's avatar
Fang Yuedong committed
303
            starImg.setOrigin(0, 0)
Xin Zhang's avatar
Xin Zhang committed
304
305
306
307
308
309
310
311
312
            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]

            if gal_origin[1] < grating_split_pos_chip < gal_end[1]:
                subSlitPos = int(grating_split_pos_chip - gal_origin[1] + 1)
                ## part img disperse

                subImg_p1 = starImg.array[:, 0:subSlitPos]
                star_p1 = galsim.Image(subImg_p1)
Fang Yuedong's avatar
Fang Yuedong committed
313
                star_p1.setOrigin(0, 0)
Xin Zhang's avatar
Xin Zhang committed
314
                origin_p1 = origin_star
xin's avatar
xin committed
315
316
                xcenter_p1 = min(x_nominal,grating_split_pos_chip-1) - 0
                ycenter_p1 = y_nominal-0
Xin Zhang's avatar
Xin Zhang committed
317
318
319
320

                sdp_p1 = SpecDisperser(orig_img=star_p1, xcenter=xcenter_p1,
                                    ycenter=ycenter_p1, origin=origin_p1,
                                    tar_spec=normalSED,
Zhang Xin's avatar
Zhang Xin committed
321
                                    band_start=brange[0], band_end=brange[1],
Xin Zhang's avatar
Xin Zhang committed
322
                                    conf=chip.sls_conf[0],
323
324
                                    isAlongY=0,
                                    flat_cube=flat_cube)
Xin Zhang's avatar
Xin Zhang committed
325

Zhang Xin's avatar
Zhang Xin committed
326
327
328
329
                # self.addSLStoChipImage(sdp=sdp_p1, chip=chip, xOrderSigPlus = xOrderSigPlus, local_wcs=chip_wcs_local)
                pos_shear = self.addSLStoChipImageWithPSF(sdp=sdp_p1, chip=chip, pos_img_local=[xcenter_p1, ycenter_p1],
                                                          psf_model=psf_model, bandNo=i + 1,
                                                          grating_split_pos=grating_split_pos,
Zhang Xin's avatar
Zhang Xin committed
330
                                                          local_wcs=chip_wcs_local, pos_img = pos_img)
Xin Zhang's avatar
Xin Zhang committed
331
332
333

                subImg_p2 = starImg.array[:, subSlitPos+1:starImg.array.shape[1]]
                star_p2 = galsim.Image(subImg_p2)
Fang Yuedong's avatar
Fang Yuedong committed
334
                star_p2.setOrigin(0, 0)
Xin Zhang's avatar
Xin Zhang committed
335
                origin_p2 = [origin_star[0], grating_split_pos_chip]
xin's avatar
xin committed
336
337
                xcenter_p2 = max(x_nominal, grating_split_pos_chip - 1) - 0
                ycenter_p2 = y_nominal - 0
Xin Zhang's avatar
Xin Zhang committed
338
339
340
341

                sdp_p2 = SpecDisperser(orig_img=star_p2, xcenter=xcenter_p2,
                                       ycenter=ycenter_p2, origin=origin_p2,
                                       tar_spec=normalSED,
Zhang Xin's avatar
Zhang Xin committed
342
                                       band_start=brange[0], band_end=brange[1],
Xin Zhang's avatar
Xin Zhang committed
343
                                       conf=chip.sls_conf[1],
344
345
                                       isAlongY=0,
                                       flat_cube=flat_cube)
Xin Zhang's avatar
Xin Zhang committed
346

Zhang Xin's avatar
Zhang Xin committed
347
348
349
350
                # self.addSLStoChipImage(sdp=sdp_p2, chip=chip, xOrderSigPlus = xOrderSigPlus, local_wcs=chip_wcs_local)
                pos_shear = self.addSLStoChipImageWithPSF(sdp=sdp_p2, chip=chip, pos_img_local=[xcenter_p2, ycenter_p2],
                                                          psf_model=psf_model, bandNo=i + 1,
                                                          grating_split_pos=grating_split_pos,
Zhang Xin's avatar
Zhang Xin committed
351
                                                          local_wcs=chip_wcs_local, pos_img = pos_img)
Xin Zhang's avatar
Xin Zhang committed
352
353
354
355

                del sdp_p1
                del sdp_p2
            elif grating_split_pos_chip<=gal_origin[1]:
xin's avatar
xin committed
356
357
                sdp = SpecDisperser(orig_img=starImg, xcenter=x_nominal - 0,
                                    ycenter=y_nominal - 0, origin=origin_star,
Xin Zhang's avatar
Xin Zhang committed
358
                                    tar_spec=normalSED,
Zhang Xin's avatar
Zhang Xin committed
359
                                    band_start=brange[0], band_end=brange[1],
Xin Zhang's avatar
Xin Zhang committed
360
                                    conf=chip.sls_conf[1],
361
362
                                    isAlongY=0,
                                    flat_cube=flat_cube)
Zhang Xin's avatar
Zhang Xin committed
363
364
365
366
                # self.addSLStoChipImage(sdp=sdp, chip=chip, xOrderSigPlus = xOrderSigPlus, local_wcs=chip_wcs_local)
                pos_shear = self.addSLStoChipImageWithPSF(sdp=sdp, chip=chip, pos_img_local=[x_nominal, y_nominal],
                                                          psf_model=psf_model, bandNo=i + 1,
                                                          grating_split_pos=grating_split_pos,
Zhang Xin's avatar
Zhang Xin committed
367
                                                          local_wcs=chip_wcs_local, pos_img = pos_img)
Xin Zhang's avatar
Xin Zhang committed
368
369
                del sdp
            elif grating_split_pos_chip>=gal_end[1]:
xin's avatar
xin committed
370
371
                sdp = SpecDisperser(orig_img=starImg, xcenter=x_nominal - 0,
                                    ycenter=y_nominal - 0, origin=origin_star,
Xin Zhang's avatar
Xin Zhang committed
372
                                    tar_spec=normalSED,
Zhang Xin's avatar
Zhang Xin committed
373
                                    band_start=brange[0], band_end=brange[1],
Xin Zhang's avatar
Xin Zhang committed
374
                                    conf=chip.sls_conf[0],
375
376
                                    isAlongY=0,
                                    flat_cube=flat_cube)
Zhang Xin's avatar
Zhang Xin committed
377
378
379
380
                # self.addSLStoChipImage(sdp=sdp, chip=chip, xOrderSigPlus = xOrderSigPlus, local_wcs=chip_wcs_local)
                pos_shear = self.addSLStoChipImageWithPSF(sdp=sdp, chip=chip, pos_img_local=[x_nominal, y_nominal],
                                                          psf_model=psf_model, bandNo=i + 1,
                                                          grating_split_pos=grating_split_pos,
Zhang Xin's avatar
Zhang Xin committed
381
                                                          local_wcs=chip_wcs_local, pos_img = pos_img)
Xin Zhang's avatar
Xin Zhang committed
382
                del sdp
Fang Yuedong's avatar
Fang Yuedong committed
383

Xin Zhang's avatar
Xin Zhang committed
384
            # print(self.y_nominal, starImg.center.y, starImg.ymin)
Zhang Xin's avatar
Zhang Xin committed
385
            # del psf
Fang Yuedong's avatar
Fang Yuedong committed
386
        return 1, pos_shear
Fang Yuedong's avatar
Fang Yuedong committed
387
388
389
390

    def getGSObj(self, psf, g1=0, g2=0, flux=None, filt=None, tel=None, exptime=150.):
        if flux == None:
            flux = self.getElectronFluxFilt(filt, tel, exptime)
Fang Yuedong's avatar
Fang Yuedong committed
391
        disk = galsim.Sersic(n=self.disk_sersic_idx, half_light_radius=self.hlr_disk, flux=1.0)
Fang Yuedong's avatar
Fang Yuedong committed
392
393
394
        disk_shape = galsim.Shear(g1=self.e1_disk, g2=self.e2_disk)
        disk = disk.shear(disk_shape)

Fang Yuedong's avatar
Fang Yuedong committed
395
        bulge = galsim.Sersic(n=self.bulge_sersic_idx, half_light_radius=self.hlr_bulge, flux=1.0)
Fang Yuedong's avatar
Fang Yuedong committed
396
397
398
399
400
401
402
403
404
405
406
407
        bulge_shape = galsim.Shear(g1=self.e1_bulge, g2=self.e2_bulge)
        bulge = bulge.shear(bulge_shape)

        gal = self.bfrac * bulge + (1.0 - self.bfrac) * disk
        gal = gal.withFlux(flux)
        gal_shear = galsim.Shear(g1=g1, g2=g2)
        gal = gal.shear(gal_shear)
        final = galsim.Convolve(psf, gal)
        return final

    def getObservedEll(self, g1=0, g2=0):
        e1_obs, e2_obs, e_obs, theta = eObs(self.e1_total, self.e2_total, g1, g2)
Zhang Xin's avatar
Zhang Xin committed
408
        return self.e1_total, self.e2_total, g1, g2, e1_obs, e2_obs