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 == 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.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)