Newer
Older
from tkinter.tix import INTEGER
import astropy.coordinates as coord
from astropy import units as u
from pylab import *
import numpy as np
import galsim
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
class Chip(object):
def __init__(self, chipID):
self.chipID = chipID
self.nchip_x = 6
self.nchip_y = 5
self.npix_tot_x = 59516
self.npix_tot_y = 49752
self.npix_gap_x = (534, 1309)
self.npix_gap_y = 898
self.cen_pix_x = 0
self.cen_pix_y = 0
self.npix_x = 9216
self.npix_y = 9232
self.pix_scale = 0.074
def getTanWCS(self, ra, dec, img_rot, pix_scale = None, xcen=None, ycen=None, logger=None):
""" Get the WCS of the image mosaic using Gnomonic/TAN projection
Parameter:
ra, dec: float
(RA, Dec) of pointing of optical axis
img_rot: galsim Angle object
Rotation of image
pix_scale: float
Pixel size in unit of as/pix
Returns:
WCS of the focal plane
"""
if logger is not None:
logger.info(" Construct the wcs of the entire image mosaic using Gnomonic/TAN projection")
if (xcen == None) or (ycen == None):
xcen = self.cen_pix_x
ycen = self.cen_pix_y
if pix_scale == None:
pix_scale = self.pix_scale
dudx = -np.cos(img_rot.rad) * pix_scale
dudy = -np.sin(img_rot.rad) * pix_scale
dvdx = -np.sin(img_rot.rad) * pix_scale
dvdy = +np.cos(img_rot.rad) * pix_scale
# dudx = +np.sin(img_rot.rad) * pix_scale
# dudy = +np.cos(img_rot.rad) * pix_scale
# dvdx = -np.cos(img_rot.rad) * pix_scale
# dvdy = +np.sin(img_rot.rad) * pix_scale
moscen = galsim.PositionD(x=xcen, y=ycen)
sky_center = galsim.CelestialCoord(ra=ra*galsim.degrees, dec=dec*galsim.degrees)
affine = galsim.AffineTransform(dudx, dudy, dvdx, dvdy, origin=moscen)
WCS = galsim.TanWCS(affine, sky_center, units=galsim.arcsec)
return WCS
def getChipRowCol(self, chipID):
rowID = ((chipID - 1) % 5) + 1
colID = 6 - ((chipID - 1) // 5)
return rowID, colID
def getChipCenter(self):
"""Calculate the edges in pixel for a given CCD chip on the focal plane
NOTE: There are 5*4 CCD chips in the focus plane for photometric observation.
Parameters:
chipID: int
the index of the chip
Returns:
A galsim BoundsD object
"""
chipID = self.chipID
rowID, colID = self.getChipRowCol(chipID)
gx1, gx2 = self.npix_gap_x
gy = self.npix_gap_y
# xlim of a given CCD chip
xrem = 2*(colID - 1) - (self.nchip_x - 1)
xcen = (self.npix_x//2 + gx1//2) * xrem
if chipID >= 26 or chipID == 21:
xcen = (self.npix_x//2 + gx1//2) * xrem - (gx2-gx1)
if chipID <= 5 or chipID == 10:
xcen = (self.npix_x//2 + gx1//2) * xrem + (gx2-gx1)
# nx0 = xcen - self.npix_x//2 + 1
# nx1 = xcen + self.npix_x//2
# ylim of a given CCD chip
yrem = (rowID - 1) - self.nchip_y // 2
ycen = (self.npix_y + gy) * yrem
# ny0 = ycen - self.npix_y//2 + 1
# ny1 = ycen + self.npix_y//2
return galsim.PositionD(xcen, ycen)
# @jit()
def getSelectPointingList(center = [60,-40], radius = 2):
points = np.loadtxt('skyMapOrSurveyList/sky.dat')
center = center#ra dec
radius = radius # degree
radii_dec = 1
radii_ra = 1/math.cos(math.pi*center[1]/180)
if radii_ra > 180:
radii_ra = 180
c_eclip = coord.SkyCoord(points[:,2]*u.degree, points[:,1]*u.degree,frame='barycentrictrueecliptic')
c_equtor = c_eclip.transform_to('icrs')
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
# print(np.min((c_equtor.ra*u.degree).value), np.max((c_equtor.ra*u.degree).value))
# c_equtor_sel = c_equtor
# points_sel = points
ra_range_lo = center[0]-radii_ra
ra_range_hi = center[0]+radii_ra
if ra_range_lo< 0:
ids1 = ((c_equtor.ra*u.degree).value<ra_range_hi) | ((c_equtor.ra*u.degree).value>360+ra_range_lo)
elif ra_range_hi > 360:
ids1 = ((c_equtor.ra*u.degree).value>ra_range_lo) | ((c_equtor.ra*u.degree).value<ra_range_hi-360)
else:
ids1 = ((c_equtor.ra*u.degree).value > ra_range_lo) & ((c_equtor.ra*u.degree).value < ra_range_hi)
dec_range_lo = center[1]-radii_dec
if center[1]-radii_dec < -90:
dec_range_lo = -90
dec_range_hi = center[1]+radii_dec
if center[1]+radii_dec > 90:
dec_range_hi = 90
ids3 = (c_equtor[ids1].dec*u.degree).value > dec_range_lo
ids4 = (c_equtor[ids1][ids3].dec*u.degree).value < dec_range_hi
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
ra = (p_.ra*u.degree).value
dec = (p_.dec*u.degree).value
# print(ra, dec)
lon = p[2]
lat = p[1]
p_result[i,0] = ra
p_result[i,1] = dec
p_result[i,2] = lon
p_result[i,3] = lat
p_result[i,4] = -113.4333
i = i + 1
return p_result
def findPointingbyChipID(chipID = 8, ra = 60., dec = -40.):
chip_center = [ra, dec]
p_list = getSelectPointingList(center = chip_center)
pchip = Chip(chipID)
p_num = p_list.shape[0]
min_d = 1000000000
r_ra = ra
r_dec = dec
for i in np.arange(0,p_num,1):
ra_n = p_list[i,0]
dec_n = p_list[i,1]
rot = p_list[i,4]*galsim.degrees
chip_wcs = pchip.getTanWCS(ra_n, dec_n, rot)
c_center = pchip.getChipCenter()
c_world = chip_wcs.toWorld(c_center)
ra_s = c_world.ra.deg
dec_s = c_world.dec.deg
# print(ra, dec, ra_s, dec_s)
d = (ra_s - ra)*(ra_s - ra) + (dec_s - dec)*(dec_s - dec)
if d < min_d:
min_d = d
r_ra = ra_n
r_dec = dec_n
return galsim.CelestialCoord(ra=r_ra*galsim.degrees,dec=r_dec*galsim.degrees)
if __name__ == "__main__":
pointing = findPointingbyChipID()
# chipID = 8
# pchip = Chip(chipID)
# ra = pointing.ra.deg
# dec = pointing.dec.deg
# rot = -113.433
# chip_wcs = pchip.getTanWCS(ra, dec, rot*galsim.degrees)
# c_center = pchip.getChipCenter()
# c_world = chip_wcs.toWorld(c_center)
# ra_s = c_world.ra.deg
# dec_s = c_world.dec.deg