# NOTE: This is a stand-alone function, meaning that you do not need # to install the entire CSST image simulation pipeline. # For a given object's coordinate (Ra, Dec), the function will predict # the object's image position and corresponding filter in the focal plane # under a specified CSST pointing centered at (rap, decp). import galsim import numpy as np import argparse import matplotlib.pyplot as plt import os, sys def focalPlaneInf(ra_target, dec_target, ra_point, dec_point, image_rot=-113.4333, figout="zTargetOnCCD.pdf"): """ Input parameters: ra_target : right ascension of the target/input object; float, in unit of degrees; dec_target: declination of the target/input object; float, in unit of degrees; ra_point : right ascension of telescope pointing center; float, in unit of degrees; dec_point : declination of telescope pointing center; float, in unit of degrees; image_rot : orientation of the camera with respect the sky; float, in unit of degrees; NOTE: image_rot=-113.4333 is the default value in current CSST image simulation; figout : location of the target object in the focal plane; str -------------------------------------------------------------- Usage: 0) specify the coordinate (ra_target, dec_target) of your target and the pointing center (ra_point dec_point) of the telescope 1) open a terminal 2) type >> python TargetLocationCheck.py ra_target dec_target ra_point dec_point or type >> python TargetLocationCheck.py ra_target dec_target ra_point dec_point -image_rot=floatNum or type >> python TargetLocationCheck.py ra_target dec_target ra_point dec_point -image_rot=floatNum -figout=FigureName """ print("^_^ Input target coordinate: [Ra, Dec] = [%10.6f, %10.6f]"%(ra_target,dec_target)) print("^_^ Input telescope pointing center: [Ra, Dec] = [%10.6f, %10.6f]"%(ra_point,dec_point)) print("^_^ Input camera orientation: %12.6f degree(s)"%image_rot) print(" ") # load ccd parameters xsize, ysize, xchip, ychip, xgap, ygap, xnchip, ynchip = ccdParam() print("^_^ Pixel range of focal plane: x = [%5d, %5d], y = [%5d, %5d]"%(-xsize/2,xsize/2,-ysize/2,ysize/2)) # wcs wcs = getTanWCS(ra_point, dec_point, image_rot, pix_scale=0.074) skyObj = galsim.CelestialCoord(ra=ra_target*galsim.degrees,dec=dec_target*galsim.degrees) pixObj = wcs.toImage(skyObj) xpixObj = pixObj.x ypixObj = pixObj.y print("^_^ Image position of target: [xImage, yImage] = [%9.3f, %9.3f]"%(xpixObj,ypixObj)) # first determine if the target is in the focal plane xin = (xpixObj+xsize/2)*(xpixObj-xsize/2) yin = (ypixObj+ysize/2)*(ypixObj-ysize/2) if xin>0 or yin>0: raise ValueError("!!! Input target is out of the focal plane") # second determine the location of the target trigger = False for i in range(30): ichip = i+1 ischip = str("0%d"%ichip)[-2:] fId, fType = getChipFilter(ichip) ix0, ix1, iy0, iy1 = getChipLim(ichip) ixin = (xpixObj-ix0)*(xpixObj-ix1) iyin = (ypixObj-iy0)*(ypixObj-iy1) if ixin<=0 and iyin<=0: trigger = True idx = xpixObj - ix0 idy = ypixObj - iy0 print(" ---------------------------------------------") print(" ** Target locates in CHIP#%s with filter %s **"%(ischip,fType)) print(" ** Target position in the chip: [x, y] = [%7.2f, %7.2f]"%(idx, idy)) print(" ---------------------------------------------") break if not trigger: print("^|^ Target locates in CCD gap") # show the figure print(" Target on CCD layout is saved into %s"%figout) ccdLayout(xpixObj, ypixObj, figout=figout) return def ccdParam(): xt, yt = 59516, 49752 x0, y0 = 9216, 9232 xgap, ygap = (534,1309), 898 xnchip, ynchip = 6, 5 ccdSize = xt, yt, x0, y0, xgap, ygap, xnchip, ynchip return ccdSize def getTanWCS(ra, dec, img_rot, pix_scale=0.074): """ 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 """ xcen, ycen = 0, 0 img_rot = img_rot * galsim.degrees 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 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 getChipFilter(chipID): """ Return the filter index and type for a given chip #(chipID) """ filter_type_list = ["nuv","u", "g", "r", "i","z","y","GU", "GV", "GI"] # TODO: maybe a more elegent way other than hard coded? # e.g. use something like a nested dict: if chipID in [6, 15, 16, 25]: filter_type = "y" if chipID in [11, 20]: filter_type = "z" if chipID in [7, 24]: filter_type = "i" if chipID in [14, 17]: filter_type = "u" if chipID in [9, 22]: filter_type = "r" if chipID in [12, 13, 18, 19]: filter_type = "nuv" if chipID in [8, 23]: filter_type = "g" if chipID in [1, 10, 21, 30]: filter_type = "GI" if chipID in [2, 5, 26, 29]: filter_type = "GV" if chipID in [3, 4, 27, 28]: filter_type = "GU" filter_id = filter_type_list.index(filter_type) return filter_id, filter_type def getChipLim(chipID): """ 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 """ xt, yt, x0, y0, gx, gy, xnchip, ynchip = ccdParam() gx1, gx2 = gx rowID = ((chipID - 1) % 5) + 1 colID = 6 - ((chipID - 1) // 5) # xlim of a given CCD chip xrem = 2*(colID - 1) - (xnchip - 1) xcen = (x0//2 + gx1//2) * xrem if chipID >= 26 or chipID == 21: xcen = (x0//2 + gx1//2) * xrem - (gx2-gx1) if chipID <= 5 or chipID == 10: xcen = (x0//2 + gx1//2) * xrem + (gx2-gx1) nx0 = xcen - x0//2 + 1 nx1 = xcen + x0//2 # ylim of a given CCD chip yrem = (rowID - 1) - ynchip // 2 ycen = (y0 + gy) * yrem ny0 = ycen - y0//2 + 1 ny1 = ycen + y0//2 return nx0-1, nx1-1, ny0-1, ny1-1 def ccdLayout(xpixTar, ypixTar, figout="ccdLayout.pdf"): fig = plt.figure(figsize=(10.0,8.0)) ax = fig.add_axes([0.1,0.1,0.80,0.80]) # plot the layout of the ccd distribution for i in range(30): ichip = i+1 fId, fType = getChipFilter(ichip) ischip = str("0%d"%ichip)[-2:] ix0, ix1, iy0, iy1 = getChipLim(ichip) ax.plot([ix0,ix1],[iy0,iy0],"k-", linewidth=2.5) ax.plot([ix0,ix1],[iy1,iy1],"k-", linewidth=2.5) ax.plot([ix0,ix0],[iy0,iy1],"k-", linewidth=2.5) ax.plot([ix1,ix1],[iy0,iy1],"k-", linewidth=2.5) ax.text(ix0+500,iy0+1500,"%s#%s"%(fType, ischip), fontsize=12, color="grey") ax.plot(xpixTar, ypixTar, "r*", ms=12) ax.set_xlabel("$X\,[\mathrm{pixels}]$", fontsize=20) ax.set_ylabel("$Y\,[\mathrm{pixels}]$", fontsize=20) ax.invert_yaxis() ax.axis('off') plt.savefig(figout) def parseArguments(): # Create argument parser parser = argparse.ArgumentParser() # Positional arguments parser.add_argument("ra_target", type=float) parser.add_argument("dec_target", type=float) parser.add_argument("ra_point", type=float) parser.add_argument("dec_point", type=float) # Optional arguments parser.add_argument("-image_rot", type=float, default=-113.4333) parser.add_argument("-figout", type=str, default="zTargetOnCCD.pdf") # Parse arguments args = parser.parse_args() return args if __name__ == "__main__": # Parse the arguments args = parseArguments() # Run function focalPlaneInf(args.ra_target, args.dec_target, args.ra_point, args.dec_point, args.image_rot, args.figout)