get_pointing.py

# 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 is None) or (ycen is 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.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
    angle = angle + 90
    if (ra < 90 or 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).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

    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 = 13, 60., -40.
    pointing = findPointingbyChipID(chipID=tchip, ra=tra, dec=tdec)
    print("[ra_center, dec_center, image_rot]: ", pointing)