From 456825709ab4635616fc88e2219522c41e8e1be8 Mon Sep 17 00:00:00 2001 From: Zhang Xin Date: Sat, 26 Oct 2024 11:23:37 +0800 Subject: [PATCH] pep8 --- .../SpecDisperser/SpecDisperser.py | 2 +- .../mock_objects/SpecDisperser/setup_c.py | 14 +-- observation_sim/psf/PSFInterpSLS.py | 84 +++++++------ tools/get_pointing.py | 111 ++++++++++-------- tools/get_pointing_accuracy.py | 57 +++++---- 5 files changed, 144 insertions(+), 124 deletions(-) diff --git a/observation_sim/mock_objects/SpecDisperser/SpecDisperser.py b/observation_sim/mock_objects/SpecDisperser/SpecDisperser.py index 2a8bfe4..a6073a7 100644 --- a/observation_sim/mock_objects/SpecDisperser/SpecDisperser.py +++ b/observation_sim/mock_objects/SpecDisperser/SpecDisperser.py @@ -336,7 +336,7 @@ class SpecDisperser(object): orders = {"A": "1st", "B": "0st", "C": "2st", "D": "-1st", "E": "-2st"} sens_file_name = conffile[0:-5] + \ "_sensitivity_" + orders[beam] + ".fits" - if not os.path.exists(sens_file_name) == True: + if os.path.exists(sens_file_name) if False: senstivity_out = Table( array([w, sens]).T, names=("WAVELENGTH", "SENSITIVITY")) senstivity_out.write(sens_file_name, format="fits") diff --git a/observation_sim/mock_objects/SpecDisperser/setup_c.py b/observation_sim/mock_objects/SpecDisperser/setup_c.py index feceb97..e5a4467 100644 --- a/observation_sim/mock_objects/SpecDisperser/setup_c.py +++ b/observation_sim/mock_objects/SpecDisperser/setup_c.py @@ -8,16 +8,16 @@ import numpy extensions = [ Extension("disperse_c.interp", ["disperse_c/interp.pyx"], - include_dirs = [numpy.get_include()], - libraries=["m"]), - + include_dirs=[numpy.get_include()], + libraries=["m"]), + Extension("disperse_c.disperse", ["disperse_c/disperse.pyx"], - include_dirs = [numpy.get_include()], - libraries=["m"]), + include_dirs=[numpy.get_include()], + libraries=["m"]), ] setup( - name = "slssim_disperse", - ext_modules = cythonize(extensions), + name="slssim_disperse", + ext_modules=cythonize(extensions), ) diff --git a/observation_sim/psf/PSFInterpSLS.py b/observation_sim/psf/PSFInterpSLS.py index 90ac027..50521b8 100644 --- a/observation_sim/psf/PSFInterpSLS.py +++ b/observation_sim/psf/PSFInterpSLS.py @@ -435,14 +435,15 @@ class PSFInterpSLS(PSFModel): # PSF_int_trans[ids_szero] = 0 # print(PSF_int_trans[ids_szero].shape[0],PSF_int_trans.shape) PSF_int_trans = PSF_int_trans/np.sum(PSF_int_trans) - ###DEBGU - ids_szero = PSF_int_trans<0 + # DEBGU + ids_szero = PSF_int_trans < 0 n01 = PSF_int_trans[ids_szero].shape[0] n1 = np.sum(np.isinf(PSF_int_trans)) n2 = np.sum(np.isnan(PSF_int_trans)) - if n1>0 or n2>0: - print("DEBUG: PSFInterpSLS, inf:%d, nan:%d, 0 num:%d"%(n1, n2, n01)) + if n1 > 0 or n2 > 0: + print("DEBUG: PSFInterpSLS, inf:%d, nan:%d, 0 num:%d" % + (n1, n2, n01)) #### # from astropy.io import fits @@ -537,36 +538,37 @@ class PSFInterpSLS(PSFModel): sumImg = np.sum(cutImg.array) tmp_img = cutImg*0 for j in np.arange(npc): - X_ = np.hstack((pos_p[:,1].flatten()[:, None], pos_p[:,0].flatten()[:, None]),dtype=np.float32) + X_ = np.hstack((pos_p[:, 1].flatten()[:, None], pos_p[:, 0].flatten()[ + :, None]), dtype=np.float32) Z_ = (pc_coeff[j].astype(np.float32)).flatten() # print(pc_coeff[j].shape[0], pos_p[:,1].shape[0], pos_p[:,0].shape[0]) cx_len = int(chip.npix_x) cy_len = int(chip.npix_y) - n_x = np.arange(0, cx_len, 1, dtype = int) - n_y = np.arange(0, cy_len, 1, dtype = int) + n_x = np.arange(0, cx_len, 1, dtype=int) + n_y = np.arange(0, cy_len, 1, dtype=int) M, N = np.meshgrid(n_x, n_y) # t1=datetime.datetime.now() # U = interpolate.griddata(X_, Z_, (M[0:cy_len, 0:cx_len],N[0:cy_len, 0:cx_len]), # method='nearest',fill_value=1.0) b_img = galsim.Image(cx_len, cy_len) - b_img.setOrigin(0,0) + b_img.setOrigin(0, 0) bounds = cutImg.bounds & b_img.bounds if bounds.area() == 0: continue # ys = cutImg.ymin - # if ys < 0: + # if ys < 0: # ys = 0 # ye = cutImg.ymin+cutImg.nrow - # if ye >= cy_len-1: + # if ye >= cy_len-1: # ye = cy_len-1 # if ye - ys <=0: # continue # xs = cutImg.xmin - # if xs < 0: + # if xs < 0: # xs = 0 # xe = cutImg.xmin+cutImg.ncol - # if xe >= cx_len-1: + # if xe >= cx_len-1: # xe = cx_len-1 # if xe - xs <=0: # continue @@ -574,10 +576,10 @@ class PSFInterpSLS(PSFModel): ye = bounds.ymax+1 xs = bounds.xmin xe = bounds.xmax+1 - U = interpolate.griddata(X_, Z_, (M[ys:ye, xs:xe],N[ys:ye, xs:xe]), - method='nearest',fill_value=1.0) + U = interpolate.griddata(X_, Z_, (M[ys:ye, xs:xe], N[ys:ye, xs:xe]), + method='nearest', fill_value=1.0) # t2=datetime.datetime.now() - + # print("time interpolate:", t2-t1) # if U.shape != cutImg.array.shape: @@ -586,53 +588,53 @@ class PSFInterpSLS(PSFModel): img_tmp = cutImg img_tmp[bounds] = img_tmp[bounds]*U psf = pcs[:, j].reshape(m_size, m_size) - tmp_img = tmp_img + signal.fftconvolve(img_tmp.array, psf, mode='same', axes=None) + tmp_img = tmp_img + \ + signal.fftconvolve(img_tmp.array, psf, mode='same', axes=None) # t3=datetime.datetime.now() # print("time convole:", t3-t2) del U del img_tmp - if np.sum(tmp_img.array)==0: + if np.sum(tmp_img.array) == 0: tmp_img = cutImg else: tmp_img = tmp_img/np.sum(tmp_img.array)*sumImg return tmp_img - def convolveFullImgWithPCAPSF(self, chip, folding_threshold=5.e-3): - keys_L1= chip_utils.getChipSLSGratingID(chip.chipID) + keys_L1 = chip_utils.getChipSLSGratingID(chip.chipID) # keys_L2 = ['order-2','order-1','order0','order1','order2'] - keys_L2 = ['order0','order1'] - keys_L3 = ['w1','w2','w3','w4'] + keys_L2 = ['order0', 'order1'] + keys_L3 = ['w1', 'w2', 'w3', 'w4'] npca = 10 x_start = chip.x_cen/chip.pix_size - chip.npix_x / 2. y_start = chip.y_cen/chip.pix_size - chip.npix_y / 2. - for i,gt in enumerate(keys_L1): + for i, gt in enumerate(keys_L1): psfCo = self.grating1_data if i > 0: psfCo = self.grating2_data for od in keys_L2: psfCo_L2 = psfCo['order1'] - if od in ['order-2','order-1','order0','order2']: + if od in ['order-2', 'order-1', 'order0', 'order2']: psfCo_L2 = psfCo['order0'] for w in keys_L3: img = chip.img_stack[gt][od][w] pcs = psfCo_L2['band'+w[1]]['band_data'][0].data - pos_p = psfCo_L2['band'+w[1]]['band_data'][1].data/chip.pix_size - np.array([y_start, x_start]) + pos_p = psfCo_L2['band'+w[1]]['band_data'][1].data / \ + chip.pix_size - np.array([y_start, x_start]) pc_coeff = psfCo_L2['band'+w[1]]['band_data'][2].data # print("DEBUG-----------",np.max(pos_p[:,1]),np.min(pos_p[:,1]), np.max(pos_p[:,0]),np.min(pos_p[:,0])) sum_img = np.sum(img.array) - - + # coeff_mat = np.zeros([npca, chip.npix_y, chip.npix_x]) # for m in np.arange(chip.npix_y): # for n in np.arange(chip.npix_x): # px = n # py = m - + # dist2 = (pos_p[:, 1] - px)*(pos_p[:, 1] - px) + (pos_p[:, 0] - py)*(pos_p[:, 0] - py) # temp_sort_dist = np.zeros([dist2.shape[0], 2]) # temp_sort_dist[:, 0] = np.arange(0, dist2.shape[0], 1) @@ -660,41 +662,45 @@ class PSFInterpSLS(PSFModel): # coeff_mat[:, m, n] = coeff_int m_size = int(pcs.shape[0]**0.5) - tmp_img = np.zeros_like(img.array,dtype=np.float32) + tmp_img = np.zeros_like(img.array, dtype=np.float32) for j in np.arange(npca): print(gt, od, w, j) - X_ = np.hstack((pos_p[:,1].flatten()[:, None], pos_p[:,0].flatten()[:, None]),dtype=np.float32) + X_ = np.hstack((pos_p[:, 1].flatten()[:, None], pos_p[:, 0].flatten()[ + :, None]), dtype=np.float32) Z_ = (pc_coeff[j].astype(np.float32)).flatten() # print(pc_coeff[j].shape[0], pos_p[:,1].shape[0], pos_p[:,0].shape[0]) sub_size = 4 cx_len = int(chip.npix_x/sub_size) cy_len = int(chip.npix_y/sub_size) - n_x = np.arange(0, chip.npix_x, sub_size, dtype = int) - n_y = np.arange(0, chip.npix_y, sub_size, dtype = int) + n_x = np.arange(0, chip.npix_x, sub_size, dtype=int) + n_y = np.arange(0, chip.npix_y, sub_size, dtype=int) M, N = np.meshgrid(n_x, n_y) - t1=datetime.datetime.now() + t1 = datetime.datetime.now() # U = interpolate.griddata(X_, Z_, (M[0:cy_len, 0:cx_len],N[0:cy_len, 0:cx_len]), # method='nearest',fill_value=1.0) U1 = interpolate.griddata(X_, Z_, (M, N), - method='nearest',fill_value=1.0) + method='nearest', fill_value=1.0) U = np.zeros_like(chip.img.array, dtype=np.float32) for mi in np.arange(cy_len): for mj in np.arange(cx_len): - U[mi*sub_size:(mi+1)*sub_size, mj*sub_size:(mj+1)*sub_size]=U1[mi,mj] - t2=datetime.datetime.now() - + U[mi*sub_size:(mi+1)*sub_size, mj * + sub_size:(mj+1)*sub_size] = U1[mi, mj] + t2 = datetime.datetime.now() + print("time interpolate:", t2-t1) img_tmp = img.array*U psf = pcs[:, j].reshape(m_size, m_size) - tmp_img = tmp_img + signal.fftconvolve(img_tmp, psf, mode='same', axes=None) + tmp_img = tmp_img + \ + signal.fftconvolve( + img_tmp, psf, mode='same', axes=None) - t3=datetime.datetime.now() + t3 = datetime.datetime.now() print("time convole:", t3-t2) del U del U1 - + chip.img = chip.img + tmp_img*sum_img/np.sum(tmp_img) del tmp_img gc.collect() diff --git a/tools/get_pointing.py b/tools/get_pointing.py index 39216a1..d5c1d17 100755 --- a/tools/get_pointing.py +++ b/tools/get_pointing.py @@ -2,7 +2,7 @@ # get_pointing # PURPOSE: # Return the pointing of CSST from a given [ra, dec] -# CALLING: +# CALLING: # pointing = findPointingbyChipID(chipID=tchip, ra=tra, dec=tdec) # INPUTS: # chipID - chip-# of CSST, int @@ -13,8 +13,8 @@ # 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 @@ -25,6 +25,7 @@ import galsim import math # from numba import jit + class Chip(object): def __init__(self, chipID): self.chipID = chipID @@ -41,9 +42,9 @@ class Chip(object): self.npix_x = 9216 self.npix_y = 9232 - self.pix_scale = 0.074 + self.pix_scale = 0.074 - def getTanWCS(self, ra, dec, img_rot, pix_scale = None, xcen=None, ycen=None, logger=None): + 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: @@ -57,7 +58,8 @@ class Chip(object): WCS of the focal plane """ if logger is not None: - logger.info(" Construct the wcs of the entire image mosaic using Gnomonic/TAN projection") + 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 @@ -72,23 +74,24 @@ class Chip(object): 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) + 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. @@ -99,8 +102,6 @@ class Chip(object): A galsim BoundsD object """ - - chipID = self.chipID rowID, colID = self.getChipRowCol(chipID) @@ -125,37 +126,42 @@ class Chip(object): 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); + 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]) + l1 = np.array([0, 0, 1]) l2 = transRaDec2D(ra, dec) - polar_ec = coord.SkyCoord(0*u.degree, 90*u.degree,frame='barycentrictrueecliptic') + 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))) + 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 + if (ra < 90 or ra > 270): + angle = -angle return angle # @jit() -def getSelectPointingList(center = [60,-40], radius = 2): + + +def getSelectPointingList(center=[60, -40], radius=2): points = np.loadtxt('sky.dat') - - center = center#ra dec - radius = radius # degree + center = center # ra dec + radius = radius # degree radii_dec = 1 radii_ra = 1/math.cos(math.pi*center[1]/180) @@ -163,7 +169,8 @@ def getSelectPointingList(center = [60,-40], radius = 2): if radii_ra > 180: radii_ra = 180 - c_eclip = coord.SkyCoord(points[:,2]*u.degree, points[:,1]*u.degree,frame='barycentrictrueecliptic') + 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)) @@ -174,13 +181,15 @@ def getSelectPointingList(center = [60,-40], radius = 2): 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) + 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_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) - + 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: @@ -189,7 +198,7 @@ def getSelectPointingList(center = [60,-40], radius = 2): 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 @@ -198,25 +207,24 @@ def getSelectPointingList(center = [60,-40], radius = 2): p_result = np.zeros([num, 5]) i = 0 - for p,p_ in zip(points[ids1][ids3][ids4],c_equtor[ids1][ids3][ids4]): + 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 + 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 + return p_result -def findPointingbyChipID(chipID = 8, ra = 60., dec = -40.): +def findPointingbyChipID(chipID=8, ra=60., dec=-40.): """_summary_ Args: @@ -226,9 +234,9 @@ def findPointingbyChipID(chipID = 8, ra = 60., dec = -40.): Returns: _type_: [ra, dec, rotation angle] - """ + """ chip_center = [ra, dec] - p_list = getSelectPointingList(center = chip_center) + p_list = getSelectPointingList(center=chip_center) pchip = Chip(chipID) p_num = p_list.shape[0] @@ -238,10 +246,10 @@ def findPointingbyChipID(chipID = 8, ra = 60., dec = -40.): 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 + 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() @@ -257,15 +265,14 @@ def findPointingbyChipID(chipID = 8, ra = 60., dec = -40.): r_ra = ra_n r_dec = dec_n r_rot = rot.deg - + if min_d == max_value: - print("RA:%f,Dec:%f不在指向范围内，请于巡天规划序列比对！！！！！"%(ra, dec)) + print("RA:%f,Dec:%f不在指向范围内，请于巡天规划序列比对！！！！！" % (ra, dec)) - return [r_ra, r_dec , r_rot] + 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) - diff --git a/tools/get_pointing_accuracy.py b/tools/get_pointing_accuracy.py index 33d8575..7c53134 100644 --- a/tools/get_pointing_accuracy.py +++ b/tools/get_pointing_accuracy.py @@ -1,6 +1,8 @@ from pylab import * -import math, sys, numpy as np +import math +import sys +import numpy as np import astropy.coordinates as coord from astropy.coordinates import SkyCoord from astropy import wcs, units as u @@ -14,8 +16,9 @@ def transRaDec2D(ra, dec): z1 = np.sin(dec / 57.2957795) return np.array([x1, y1, z1]) + def ecl2radec(lon_ecl, lat_ecl): - ## convert from ecliptic coordinates to equatorial coordinates + # convert from ecliptic coordinates to equatorial coordinates c_ecl = SkyCoord( lon=lon_ecl * u.degree, lat=lat_ecl * u.degree, frame="barycentrictrueecliptic" ) @@ -25,25 +28,26 @@ def ecl2radec(lon_ecl, lat_ecl): def radec2ecl(ra, dec): - ## convert from equatorial coordinates to ecliptic coordinates + # convert from equatorial coordinates to ecliptic coordinates c_eq = SkyCoord(ra=ra * u.degree, dec=dec * u.degree, frame="icrs") c_ecl = c_eq.transform_to("barycentrictrueecliptic") lon_ecl, lat_ecl = c_ecl.lon.degree, c_ecl.lat.degree return lon_ecl, lat_ecl -def cal_FoVcenter_1P_equatorial(ra_FieldCenter, dec_FieldCenter, chipID = 1, pa = -23.5): - ### [ra_FieldCenter, dec_FieldCenter] is the center ra, dec of calibration fileds, such as: NEP, NGC 6397, etc. - ### [ra_ChipCenter, dec_ChipCenter] is the center ra, dec of the Chip center. - ### [ra_PointCenter, dec_PointCenter] is the telescope pointing center. - ## Calculate PA angle +def cal_FoVcenter_1P_equatorial(ra_FieldCenter, dec_FieldCenter, chipID=1, pa=-23.5): + + # [ra_FieldCenter, dec_FieldCenter] is the center ra, dec of calibration fileds, such as: NEP, NGC 6397, etc. + # [ra_ChipCenter, dec_ChipCenter] is the center ra, dec of the Chip center. + # [ra_PointCenter, dec_PointCenter] is the telescope pointing center. + # Calculate PA angle chip = Chip(chipID) h_ext = ImageHeader.generateExtensionHeader( chip=chip, xlen=chip.npix_x, ylen=chip.npix_y, - ra=(ra_FieldCenter * u.degree).value, + ra=(ra_FieldCenter * u.degree).value, dec=(dec_FieldCenter * u.degree).value, pa=pa, gain=chip.gain, @@ -54,7 +58,7 @@ def cal_FoVcenter_1P_equatorial(ra_FieldCenter, dec_FieldCenter, chipID = 1, pa pixel_size=chip.pix_size, xcen=chip.x_cen, ycen=chip.y_cen, - ) + ) # assume that the telescope point to the sky center; then abtain the chip center under this situation; then calculate the differenc between the assumed telescope center and chip center; then add this difference to the sky center @@ -75,24 +79,25 @@ def cal_FoVcenter_1P_equatorial(ra_FieldCenter, dec_FieldCenter, chipID = 1, pa return ra_PointCenter, dec_PointCenter, lon_ecl_PointCenter, lat_ecl_PointCenter -def cal_FoVcenter_1P_ecliptic(lon_ecl_FieldCenter, lat_ecl_FieldCenter, chipID = 1, pa = -23.5): - ### [ra_FieldCenter, dec_FieldCenter] is the center ra, dec of calibration fileds, such as: NEP, NGC 6397, etc. - ### [ra_ChipCenter, dec_ChipCenter] is the center ra, dec of the Chip center. - ### [ra_PointCenter, dec_PointCenter] is the telescope pointing center. +def cal_FoVcenter_1P_ecliptic(lon_ecl_FieldCenter, lat_ecl_FieldCenter, chipID=1, pa=-23.5): + + # [ra_FieldCenter, dec_FieldCenter] is the center ra, dec of calibration fileds, such as: NEP, NGC 6397, etc. + # [ra_ChipCenter, dec_ChipCenter] is the center ra, dec of the Chip center. + # [ra_PointCenter, dec_PointCenter] is the telescope pointing center. ra_FieldCenter, dec_FieldCenter = ecl2radec( lon_ecl_FieldCenter, lat_ecl_FieldCenter ) - ## Calculate PA angle + # Calculate PA angle chip = Chip(chipID) h_ext = ImageHeader.generateExtensionHeader( chip=chip, xlen=chip.npix_x, ylen=chip.npix_y, - ra=(ra_FieldCenter * u.degree).value, + ra=(ra_FieldCenter * u.degree).value, dec=(dec_FieldCenter * u.degree).value, pa=pa, gain=chip.gain, @@ -103,7 +108,7 @@ def cal_FoVcenter_1P_ecliptic(lon_ecl_FieldCenter, lat_ecl_FieldCenter, chipID = pixel_size=chip.pix_size, xcen=chip.x_cen, ycen=chip.y_cen, - ) + ) # assume that the telescope point to the sky center; then abtain the chip center under this situation; then calculate the differenc between the assumed telescope center and chip center; then add this difference to the sky center @@ -124,14 +129,15 @@ def cal_FoVcenter_1P_ecliptic(lon_ecl_FieldCenter, lat_ecl_FieldCenter, chipID = return ra_PointCenter, dec_PointCenter, lon_ecl_PointCenter, lat_ecl_PointCenter -def getChipCenterRaDec(chipID = 1, p_ra = 60., p_dec = -40.): + +def getChipCenterRaDec(chipID=1, p_ra=60., p_dec=-40.): chip = Chip(chipID) h_ext = ImageHeader.generateExtensionHeader( chip=chip, xlen=chip.npix_x, ylen=chip.npix_y, - ra=(p_ra * u.degree).value, + ra=(p_ra * u.degree).value, dec=(p_dec * u.degree).value, pa=pa, gain=chip.gain, @@ -142,24 +148,25 @@ def getChipCenterRaDec(chipID = 1, p_ra = 60., p_dec = -40.): pixel_size=chip.pix_size, xcen=chip.x_cen, ycen=chip.y_cen, - ) + ) wcs_h = wcs.WCS(h_ext) pixs = np.array([[4608, 4616]]) world_point = wcs_h.wcs_pix2world(pixs, 0) RA_chip, Dec_chip = world_point[0][0], world_point[0][1] return RA_chip, Dec_chip + if __name__ == '__main__': ra_input, dec_input = 270.00000, 66.56000 # NEP pa = 23.5 # chipid = 2 - for chipid in np.arange(1,31,1): + for chipid in np.arange(1, 31, 1): ra, dec, lon_ecl, lat_ecl = cal_FoVcenter_1P_equatorial( - ra_input, dec_input,chipID=chipid, pa=pa) + ra_input, dec_input, chipID=chipid, pa=pa) - print("chip id is %d, chip center [ra,dec] is [%f, %f], pointing center calculated [ra,dec] is [%f, %f]"%(chipid, ra_input, dec_input, ra, dec)) - #for check the result + print("chip id is %d, chip center [ra,dec] is [%f, %f], pointing center calculated [ra,dec] is [%f, %f]" % ( + chipid, ra_input, dec_input, ra, dec)) + # for check the result # testRA, testDec = getChipCenterRaDec(chipID = chipid, p_ra = ra, p_dec = dec) # print(ra_input-testRA, dec_input-testDec) - -- GitLab