Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
import os, sys
import random
import numpy as np
import astropy.constants as cons
from astropy.table import Table
from scipy import interpolate
import astropy.io.fits as fitsio
import galsim
import gc
from ObservationSim.MockObject.MockObject import MockObject
from ObservationSim.MockObject._util import magToFlux,VC_A
from ObservationSim.MockObject._util import eObs, integrate_sed_bandpass, getNormFactorForSpecWithABMAG, getObservedSED, getABMAG,convolveGaussXorders
class Stamp(MockObject):
def __init__(self, param):
super().__init__(param)
def unload_SED(self):
"""(Test) free up SED memory
"""
del self.sed
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):
if nphotons_tot == None:
nphotons_tot = self.getElectronFluxFilt(filt, tel, exptime)
try:
full = integrate_sed_bandpass(sed=self.sed, bandpass=filt.bandpass_full)
except Exception as e:
print(e)
self.logger.error(e)
return False
nphotons_sum = 0
photons_list = []
xmax, ymax = 0, 0
if self.getMagFilter(filt) <= 15:
folding_threshold = 5.e-4
else:
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]
try:
sub = integrate_sed_bandpass(sed=self.sed, bandpass=bandpass)
except Exception as e:
print(e)
self.logger.error(e)
# return False
continue
ratio = sub/full
if not (ratio == -1 or (ratio != ratio)):
nphotons = ratio * nphotons_tot
else:
# return False
continue
nphotons_sum += nphotons
psf, pos_shear = psf_model.get_PSF(chip=chip, pos_img=pos_img, bandpass=bandpass, folding_threshold=folding_threshold)
_gal = self.param['image']
galIm = galsim.ImageF(_gal, scale=self.param['pixScale'])
gal = galsim.InterpolatedImage(galIm)
gal = gal.withFlux(nphotons)
#gal_shear = galsim.Shear(g1=g1, g2=g2)
#gal = gal.shear(gal_shear)
gal = galsim.Convolve(psf, gal)
if fd_shear is not None:
gal = gal.shear(fd_shear)
stamp = gal.drawImage(wcs=real_wcs_local, method='phot', offset=self.offset, save_photons=True)
xmax = max(xmax, stamp.xmax - stamp.xmin)
ymax = max(ymax, stamp.ymax - stamp.ymin)
photons = stamp.photons
photons.x += x_nominal
photons.y += y_nominal
photons_list.append(photons)
del gal
# print('xmax = %d, ymax = %d '%(xmax, ymax))
stamp = galsim.ImageF(int(xmax*1.1), int(ymax*1.1))
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)
if bounds.area() > 0:
chip.img.setOrigin(0, 0)
stamp[bounds] = chip.img[bounds]
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)
chip.img[bounds] = stamp[bounds]
chip.img.setOrigin(chip.bound.xmin, chip.bound.ymin)
del photons_list
del stamp
gc.collect()
return True, pos_shear