# NAME: # get_pointing # PURPOSE: # Return the pointing of CSST from a given [ra, dec] # CALLING: # pointing = findPointingbyChipID(chipID=tchip, ra=tra, dec=tdec) # INPUTS: # chipID - chip-# of CSST, int # ra - Right ascension in decimal degrees, float # dec- Declination in decimal degrees, float # OUTPUTS: # pointing - [ra_center, dec_center, image_rot] # HISTORY: # Written by Xin Zhang, 23 Apr. 2023 # Included by csst-simulation, C.W. 25 Apr. 2023 # # 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 import math # from numba import jit 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) def transRaDec2D(ra, dec): x1 = np.cos(dec / 57.2957795) * np.cos(ra / 57.2957795); y1 = np.cos(dec / 57.2957795) * np.sin(ra / 57.2957795); z1 = np.sin(dec / 57.2957795); return np.array([x1, y1, z1]) def getobsPA(ra, dec): l1 = np.array([0,0,1]) l2 = transRaDec2D(ra, dec) polar_ec = coord.SkyCoord(0*u.degree, 90*u.degree,frame='barycentrictrueecliptic') polar_eq = polar_ec.transform_to('icrs') # print(polar_eq.ra.value,polar_eq.dec.value) polar_d = transRaDec2D(polar_eq.ra.value, polar_eq.dec.value) l1l2cross = np.cross(l2,l1) pdl2cross = np.cross(l2,polar_d) angle = math.acos(np.dot(l1l2cross,pdl2cross)/(np.linalg.norm(l1l2cross)*np.linalg.norm(pdl2cross))) angle = (angle)/math.pi*180 # if (ra>90 and ra< 270): # angle=-angle return angle # @jit() def getSelectPointingList(center = [60,-40], radius = 2): points = np.loadtxt('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') # 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).value360+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_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 num = points[ids1][ids3][ids4].shape[0] p_result = np.zeros([num, 5]) i = 0 for p,p_ in zip(points[ids1][ids3][ids4],c_equtor[ids1][ids3][ids4]): 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] = getobsPA(ra, dec) + 90 i = i + 1 return p_result def findPointingbyChipID(chipID = 8, ra = 60., dec = -40.): """_summary_ Args: chipID (int, optional): Chip ID. ra (_type_, optional): Chip center ra. dec (_type_, optional): Chip center dec. Returns: _type_: [ra, dec, rotation angle] """ chip_center = [ra, dec] p_list = getSelectPointingList(center = chip_center) pchip = Chip(chipID) p_num = p_list.shape[0] max_value = 1000000000 min_d = max_value r_ra = ra r_dec = dec r_rot = 0. 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 r_rot = rot.deg if min_d == max_value: print("RA:%f,Dec:%f不在指向范围内，请于巡天规划序列比对！！！！！"%(ra, dec)) return [r_ra, r_dec , r_rot] if __name__ == "__main__": tchip, tra, tdec = 8, 60., -40. pointing = findPointingbyChipID(chipID=tchip, ra=tra, dec=tdec) print("[ra_center, dec_center, image_rot]: ", pointing)