Skip to content
GitLab
Commit 45682570 authored by Zhang Xin's avatar Zhang Xin
Browse files

pep8

parent d10eb6a8
Pipeline #7114 failed with stage
in 0 seconds
Hide whitespace changes
Inline Side-by-side
observation_sim/mock_objects/SpecDisperser/SpecDisperser.py
View file @ 45682570
...@@ -336,7 +336,7 @@ class SpecDisperser(object): ...@@ -336,7 +336,7 @@ class SpecDisperser(object):
orders = {"A": "1st", "B": "0st", "C": "2st", "D": "-1st", "E": "-2st"} orders = {"A": "1st", "B": "0st", "C": "2st", "D": "-1st", "E": "-2st"}
sens_file_name = conffile[0:-5] + \ sens_file_name = conffile[0:-5] + \
"_sensitivity_" + orders[beam] + ".fits" "_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( senstivity_out = Table(
array([w, sens]).T, names=("WAVELENGTH", "SENSITIVITY")) array([w, sens]).T, names=("WAVELENGTH", "SENSITIVITY"))
senstivity_out.write(sens_file_name, format="fits") senstivity_out.write(sens_file_name, format="fits")
......
observation_sim/mock_objects/SpecDisperser/setup_c.py
View file @ 45682570
...@@ -8,16 +8,16 @@ import numpy ...@@ -8,16 +8,16 @@ import numpy
extensions = [ extensions = [
Extension("disperse_c.interp", ["disperse_c/interp.pyx"], Extension("disperse_c.interp", ["disperse_c/interp.pyx"],
include_dirs = [numpy.get_include()], include_dirs=[numpy.get_include()],
libraries=["m"]), libraries=["m"]),
Extension("disperse_c.disperse", ["disperse_c/disperse.pyx"], Extension("disperse_c.disperse", ["disperse_c/disperse.pyx"],
include_dirs = [numpy.get_include()], include_dirs=[numpy.get_include()],
libraries=["m"]), libraries=["m"]),
] ]
setup( setup(
name = "slssim_disperse", name="slssim_disperse",
ext_modules = cythonize(extensions), ext_modules=cythonize(extensions),
) )
observation_sim/psf/PSFInterpSLS.py
View file @ 45682570
...@@ -435,14 +435,15 @@ class PSFInterpSLS(PSFModel): ...@@ -435,14 +435,15 @@ class PSFInterpSLS(PSFModel):
# PSF_int_trans[ids_szero] = 0 # PSF_int_trans[ids_szero] = 0
# print(PSF_int_trans[ids_szero].shape[0],PSF_int_trans.shape) # print(PSF_int_trans[ids_szero].shape[0],PSF_int_trans.shape)
PSF_int_trans = PSF_int_trans/np.sum(PSF_int_trans) PSF_int_trans = PSF_int_trans/np.sum(PSF_int_trans)
###DEBGU # DEBGU
ids_szero = PSF_int_trans<0 ids_szero = PSF_int_trans < 0
n01 = PSF_int_trans[ids_szero].shape[0] n01 = PSF_int_trans[ids_szero].shape[0]
n1 = np.sum(np.isinf(PSF_int_trans)) n1 = np.sum(np.isinf(PSF_int_trans))
n2 = np.sum(np.isnan(PSF_int_trans)) n2 = np.sum(np.isnan(PSF_int_trans))
if n1>0 or n2>0: if n1 > 0 or n2 > 0:
print("DEBUG: PSFInterpSLS, inf:%d, nan:%d, 0 num:%d"%(n1, n2, n01)) print("DEBUG: PSFInterpSLS, inf:%d, nan:%d, 0 num:%d" %
(n1, n2, n01))
#### ####
# from astropy.io import fits # from astropy.io import fits
...@@ -537,36 +538,37 @@ class PSFInterpSLS(PSFModel): ...@@ -537,36 +538,37 @@ class PSFInterpSLS(PSFModel):
sumImg = np.sum(cutImg.array) sumImg = np.sum(cutImg.array)
tmp_img = cutImg*0 tmp_img = cutImg*0
for j in np.arange(npc): 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() Z_ = (pc_coeff[j].astype(np.float32)).flatten()
# print(pc_coeff[j].shape[0], pos_p[:,1].shape[0], pos_p[:,0].shape[0]) # print(pc_coeff[j].shape[0], pos_p[:,1].shape[0], pos_p[:,0].shape[0])
cx_len = int(chip.npix_x) cx_len = int(chip.npix_x)
cy_len = int(chip.npix_y) cy_len = int(chip.npix_y)
n_x = np.arange(0, cx_len, 1, dtype = int) n_x = np.arange(0, cx_len, 1, dtype=int)
n_y = np.arange(0, cy_len, 1, dtype = int) n_y = np.arange(0, cy_len, 1, dtype=int)
M, N = np.meshgrid(n_x, n_y) 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]), # 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) # method='nearest',fill_value=1.0)
b_img = galsim.Image(cx_len, cy_len) b_img = galsim.Image(cx_len, cy_len)
b_img.setOrigin(0,0) b_img.setOrigin(0, 0)
bounds = cutImg.bounds & b_img.bounds bounds = cutImg.bounds & b_img.bounds
if bounds.area() == 0: if bounds.area() == 0:
continue continue
# ys = cutImg.ymin # ys = cutImg.ymin
# if ys < 0: # if ys < 0:
# ys = 0 # ys = 0
# ye = cutImg.ymin+cutImg.nrow # ye = cutImg.ymin+cutImg.nrow
# if ye >= cy_len-1: # if ye >= cy_len-1:
# ye = cy_len-1 # ye = cy_len-1
# if ye - ys <=0: # if ye - ys <=0:
# continue # continue
# xs = cutImg.xmin # xs = cutImg.xmin
# if xs < 0: # if xs < 0:
# xs = 0 # xs = 0
# xe = cutImg.xmin+cutImg.ncol # xe = cutImg.xmin+cutImg.ncol
# if xe >= cx_len-1: # if xe >= cx_len-1:
# xe = cx_len-1 # xe = cx_len-1
# if xe - xs <=0: # if xe - xs <=0:
# continue # continue
...@@ -574,10 +576,10 @@ class PSFInterpSLS(PSFModel): ...@@ -574,10 +576,10 @@ class PSFInterpSLS(PSFModel):
ye = bounds.ymax+1 ye = bounds.ymax+1
xs = bounds.xmin xs = bounds.xmin
xe = bounds.xmax+1 xe = bounds.xmax+1
U = interpolate.griddata(X_, Z_, (M[ys:ye, xs:xe],N[ys:ye, xs:xe]), U = interpolate.griddata(X_, Z_, (M[ys:ye, xs:xe], N[ys:ye, xs:xe]),
method='nearest',fill_value=1.0) method='nearest', fill_value=1.0)
# t2=datetime.datetime.now() # t2=datetime.datetime.now()
# print("time interpolate:", t2-t1) # print("time interpolate:", t2-t1)
# if U.shape != cutImg.array.shape: # if U.shape != cutImg.array.shape:
...@@ -586,53 +588,53 @@ class PSFInterpSLS(PSFModel): ...@@ -586,53 +588,53 @@ class PSFInterpSLS(PSFModel):
img_tmp = cutImg img_tmp = cutImg
img_tmp[bounds] = img_tmp[bounds]*U img_tmp[bounds] = img_tmp[bounds]*U
psf = pcs[:, j].reshape(m_size, m_size) 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() # t3=datetime.datetime.now()
# print("time convole:", t3-t2) # print("time convole:", t3-t2)
del U del U
del img_tmp del img_tmp
if np.sum(tmp_img.array)==0: if np.sum(tmp_img.array) == 0:
tmp_img = cutImg tmp_img = cutImg
else: else:
tmp_img = tmp_img/np.sum(tmp_img.array)*sumImg tmp_img = tmp_img/np.sum(tmp_img.array)*sumImg
return tmp_img return tmp_img
def convolveFullImgWithPCAPSF(self, chip, folding_threshold=5.e-3): 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 = ['order-2','order-1','order0','order1','order2']
keys_L2 = ['order0','order1'] keys_L2 = ['order0', 'order1']
keys_L3 = ['w1','w2','w3','w4'] keys_L3 = ['w1', 'w2', 'w3', 'w4']
npca = 10 npca = 10
x_start = chip.x_cen/chip.pix_size - chip.npix_x / 2. x_start = chip.x_cen/chip.pix_size - chip.npix_x / 2.
y_start = chip.y_cen/chip.pix_size - chip.npix_y / 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 psfCo = self.grating1_data
if i > 0: if i > 0:
psfCo = self.grating2_data psfCo = self.grating2_data
for od in keys_L2: for od in keys_L2:
psfCo_L2 = psfCo['order1'] 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'] psfCo_L2 = psfCo['order0']
for w in keys_L3: for w in keys_L3:
img = chip.img_stack[gt][od][w] img = chip.img_stack[gt][od][w]
pcs = psfCo_L2['band'+w[1]]['band_data'][0].data 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 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])) # 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) sum_img = np.sum(img.array)
# coeff_mat = np.zeros([npca, chip.npix_y, chip.npix_x]) # coeff_mat = np.zeros([npca, chip.npix_y, chip.npix_x])
# for m in np.arange(chip.npix_y): # for m in np.arange(chip.npix_y):
# for n in np.arange(chip.npix_x): # for n in np.arange(chip.npix_x):
# px = n # px = n
# py = m # py = m
# dist2 = (pos_p[:, 1] - px)*(pos_p[:, 1] - px) + (pos_p[:, 0] - py)*(pos_p[:, 0] - py) # 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 = np.zeros([dist2.shape[0], 2])
# temp_sort_dist[:, 0] = np.arange(0, dist2.shape[0], 1) # temp_sort_dist[:, 0] = np.arange(0, dist2.shape[0], 1)
...@@ -660,41 +662,45 @@ class PSFInterpSLS(PSFModel): ...@@ -660,41 +662,45 @@ class PSFInterpSLS(PSFModel):
# coeff_mat[:, m, n] = coeff_int # coeff_mat[:, m, n] = coeff_int
m_size = int(pcs.shape[0]**0.5) 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): for j in np.arange(npca):
print(gt, od, w, j) 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() Z_ = (pc_coeff[j].astype(np.float32)).flatten()
# print(pc_coeff[j].shape[0], pos_p[:,1].shape[0], pos_p[:,0].shape[0]) # print(pc_coeff[j].shape[0], pos_p[:,1].shape[0], pos_p[:,0].shape[0])
sub_size = 4 sub_size = 4
cx_len = int(chip.npix_x/sub_size) cx_len = int(chip.npix_x/sub_size)
cy_len = int(chip.npix_y/sub_size) cy_len = int(chip.npix_y/sub_size)
n_x = np.arange(0, chip.npix_x, 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) n_y = np.arange(0, chip.npix_y, sub_size, dtype=int)
M, N = np.meshgrid(n_x, n_y) 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]), # 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) # method='nearest',fill_value=1.0)
U1 = interpolate.griddata(X_, Z_, (M, N), 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) U = np.zeros_like(chip.img.array, dtype=np.float32)
for mi in np.arange(cy_len): for mi in np.arange(cy_len):
for mj in np.arange(cx_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] U[mi*sub_size:(mi+1)*sub_size, mj *
t2=datetime.datetime.now() sub_size:(mj+1)*sub_size] = U1[mi, mj]
t2 = datetime.datetime.now()
print("time interpolate:", t2-t1) print("time interpolate:", t2-t1)
img_tmp = img.array*U img_tmp = img.array*U
psf = pcs[:, j].reshape(m_size, m_size) 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) print("time convole:", t3-t2)
del U del U
del U1 del U1
chip.img = chip.img + tmp_img*sum_img/np.sum(tmp_img) chip.img = chip.img + tmp_img*sum_img/np.sum(tmp_img)
del tmp_img del tmp_img
gc.collect() gc.collect()
......
tools/get_pointing.py
View file @ 45682570
...@@ -2,7 +2,7 @@ ...@@ -2,7 +2,7 @@
# get_pointing # get_pointing
# PURPOSE: # PURPOSE:
# Return the pointing of CSST from a given [ra, dec] # Return the pointing of CSST from a given [ra, dec]
# CALLING: # CALLING:
# pointing = findPointingbyChipID(chipID=tchip, ra=tra, dec=tdec) # pointing = findPointingbyChipID(chipID=tchip, ra=tra, dec=tdec)
# INPUTS: # INPUTS:
# chipID - chip-# of CSST, int # chipID - chip-# of CSST, int
...@@ -13,8 +13,8 @@ ...@@ -13,8 +13,8 @@
# HISTORY: # HISTORY:
# Written by Xin Zhang, 23 Apr. 2023 # Written by Xin Zhang, 23 Apr. 2023
# Included by csst-simulation, C.W. 25 Apr. 2023 # Included by csst-simulation, C.W. 25 Apr. 2023
# #
# #
from tkinter.tix import INTEGER from tkinter.tix import INTEGER
import astropy.coordinates as coord import astropy.coordinates as coord
...@@ -25,6 +25,7 @@ import galsim ...@@ -25,6 +25,7 @@ import galsim
import math import math
# from numba import jit # from numba import jit
class Chip(object): class Chip(object):
def __init__(self, chipID): def __init__(self, chipID):
self.chipID = chipID self.chipID = chipID
...@@ -41,9 +42,9 @@ class Chip(object): ...@@ -41,9 +42,9 @@ class Chip(object):
self.npix_x = 9216 self.npix_x = 9216
self.npix_y = 9232 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 """ Get the WCS of the image mosaic using Gnomonic/TAN projection
Parameter: Parameter:
...@@ -57,7 +58,8 @@ class Chip(object): ...@@ -57,7 +58,8 @@ class Chip(object):
WCS of the focal plane WCS of the focal plane
""" """
if logger is not None: 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): if (xcen == None) or (ycen == None):
xcen = self.cen_pix_x xcen = self.cen_pix_x
ycen = self.cen_pix_y ycen = self.cen_pix_y
...@@ -72,23 +74,24 @@ class Chip(object): ...@@ -72,23 +74,24 @@ class Chip(object):
dudy = +np.sin(img_rot.rad) * pix_scale dudy = +np.sin(img_rot.rad) * pix_scale
dvdx = -np.sin(img_rot.rad) * pix_scale dvdx = -np.sin(img_rot.rad) * pix_scale
dvdy = -np.cos(img_rot.rad) * pix_scale dvdy = -np.cos(img_rot.rad) * pix_scale
# dudx = +np.sin(img_rot.rad) * pix_scale # dudx = +np.sin(img_rot.rad) * pix_scale
# dudy = +np.cos(img_rot.rad) * pix_scale # dudy = +np.cos(img_rot.rad) * pix_scale
# dvdx = -np.cos(img_rot.rad) * pix_scale # dvdx = -np.cos(img_rot.rad) * pix_scale
# dvdy = +np.sin(img_rot.rad) * pix_scale # dvdy = +np.sin(img_rot.rad) * pix_scale
moscen = galsim.PositionD(x=xcen, y=ycen) 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) affine = galsim.AffineTransform(dudx, dudy, dvdx, dvdy, origin=moscen)
WCS = galsim.TanWCS(affine, sky_center, units=galsim.arcsec) WCS = galsim.TanWCS(affine, sky_center, units=galsim.arcsec)
return WCS return WCS
def getChipRowCol(self, chipID): def getChipRowCol(self, chipID):
rowID = ((chipID - 1) % 5) + 1 rowID = ((chipID - 1) % 5) + 1
colID = 6 - ((chipID - 1) // 5) colID = 6 - ((chipID - 1) // 5)
return rowID, colID return rowID, colID
def getChipCenter(self): def getChipCenter(self):
"""Calculate the edges in pixel for a given CCD chip on the focal plane """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. NOTE: There are 5*4 CCD chips in the focus plane for photometric observation.
...@@ -99,8 +102,6 @@ class Chip(object): ...@@ -99,8 +102,6 @@ class Chip(object):
A galsim BoundsD object A galsim BoundsD object
""" """
chipID = self.chipID chipID = self.chipID
rowID, colID = self.getChipRowCol(chipID) rowID, colID = self.getChipRowCol(chipID)
...@@ -125,37 +126,42 @@ class Chip(object): ...@@ -125,37 +126,42 @@ class Chip(object):
return galsim.PositionD(xcen, ycen) return galsim.PositionD(xcen, ycen)
def transRaDec2D(ra, dec): def transRaDec2D(ra, dec):
x1 = np.cos(dec / 57.2957795) * np.cos(ra / 57.2957795); x1 = np.cos(dec / 57.2957795) * np.cos(ra / 57.2957795)
y1 = np.cos(dec / 57.2957795) * np.sin(ra / 57.2957795); y1 = np.cos(dec / 57.2957795) * np.sin(ra / 57.2957795)
z1 = np.sin(dec / 57.2957795); z1 = np.sin(dec / 57.2957795)
return np.array([x1, y1, z1]) return np.array([x1, y1, z1])
def getobsPA(ra, dec): def getobsPA(ra, dec):
l1 = np.array([0,0,1]) l1 = np.array([0, 0, 1])
l2 = transRaDec2D(ra, dec) 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') polar_eq = polar_ec.transform_to('icrs')
# print(polar_eq.ra.value,polar_eq.dec.value) # print(polar_eq.ra.value,polar_eq.dec.value)
polar_d = transRaDec2D(polar_eq.ra.value, polar_eq.dec.value) polar_d = transRaDec2D(polar_eq.ra.value, polar_eq.dec.value)
l1l2cross = np.cross(l2,l1) l1l2cross = np.cross(l2, l1)
pdl2cross = np.cross(l2,polar_d) pdl2cross = np.cross(l2, polar_d)
angle = math.acos(np.dot(l1l2cross,pdl2cross)/(np.linalg.norm(l1l2cross)*np.linalg.norm(pdl2cross))) angle = math.acos(np.dot(l1l2cross, pdl2cross) /
(np.linalg.norm(l1l2cross)*np.linalg.norm(pdl2cross)))
angle = (angle)/math.pi*180 angle = (angle)/math.pi*180
angle = angle + 90 angle = angle + 90
if (ra<90 or ra> 270): if (ra < 90 or ra > 270):
angle=-angle angle = -angle
return angle return angle
# @jit() # @jit()
def getSelectPointingList(center = [60,-40], radius = 2):
def getSelectPointingList(center=[60, -40], radius=2):
points = np.loadtxt('sky.dat') points = np.loadtxt('sky.dat')
center = center # ra dec
center = center#ra dec radius = radius # degree
radius = radius # degree
radii_dec = 1 radii_dec = 1
radii_ra = 1/math.cos(math.pi*center[1]/180) radii_ra = 1/math.cos(math.pi*center[1]/180)
...@@ -163,7 +169,8 @@ def getSelectPointingList(center = [60,-40], radius = 2): ...@@ -163,7 +169,8 @@ def getSelectPointingList(center = [60,-40], radius = 2):
if radii_ra > 180: if radii_ra > 180:
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') c_equtor = c_eclip.transform_to('icrs')
# print(np.min((c_equtor.ra*u.degree).value), np.max((c_equtor.ra*u.degree).value)) # 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): ...@@ -174,13 +181,15 @@ def getSelectPointingList(center = [60,-40], radius = 2):
ra_range_lo = center[0]-radii_ra ra_range_lo = center[0]-radii_ra
ra_range_hi = center[0]+radii_ra ra_range_hi = center[0]+radii_ra
if ra_range_lo< 0: 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) 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: 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) ids1 = ((c_equtor.ra*u.degree).value >
ra_range_lo) | ((c_equtor.ra*u.degree).value < ra_range_hi-360)
else: 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 dec_range_lo = center[1]-radii_dec
if center[1]-radii_dec < -90: if center[1]-radii_dec < -90:
...@@ -189,7 +198,7 @@ def getSelectPointingList(center = [60,-40], radius = 2): ...@@ -189,7 +198,7 @@ def getSelectPointingList(center = [60,-40], radius = 2):
dec_range_hi = center[1]+radii_dec dec_range_hi = center[1]+radii_dec
if center[1]+radii_dec > 90: if center[1]+radii_dec > 90:
dec_range_hi = 90 dec_range_hi = 90
ids3 = (c_equtor[ids1].dec*u.degree).value > dec_range_lo ids3 = (c_equtor[ids1].dec*u.degree).value > dec_range_lo
ids4 = (c_equtor[ids1][ids3].dec*u.degree).value < dec_range_hi ids4 = (c_equtor[ids1][ids3].dec*u.degree).value < dec_range_hi
...@@ -198,25 +207,24 @@ def getSelectPointingList(center = [60,-40], radius = 2): ...@@ -198,25 +207,24 @@ def getSelectPointingList(center = [60,-40], radius = 2):
p_result = np.zeros([num, 5]) p_result = np.zeros([num, 5])
i = 0 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 ra = (p_.ra*u.degree).value
dec = (p_.dec*u.degree).value dec = (p_.dec*u.degree).value
# print(ra, dec) # print(ra, dec)
lon = p[2] lon = p[2]
lat = p[1] lat = p[1]
p_result[i,0] = ra p_result[i, 0] = ra
p_result[i,1] = dec p_result[i, 1] = dec
p_result[i,2] = lon p_result[i, 2] = lon
p_result[i,3] = lat p_result[i, 3] = lat
p_result[i,4] = getobsPA(ra, dec) + 90 p_result[i, 4] = getobsPA(ra, dec) + 90
i = i + 1 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_ """_summary_
Args: Args:
...@@ -226,9 +234,9 @@ def findPointingbyChipID(chipID = 8, ra = 60., dec = -40.): ...@@ -226,9 +234,9 @@ def findPointingbyChipID(chipID = 8, ra = 60., dec = -40.):
Returns: Returns:
_type_: [ra, dec, rotation angle] _type_: [ra, dec, rotation angle]
""" """
chip_center = [ra, dec] chip_center = [ra, dec]
p_list = getSelectPointingList(center = chip_center) p_list = getSelectPointingList(center=chip_center)
pchip = Chip(chipID) pchip = Chip(chipID)
p_num = p_list.shape[0] p_num = p_list.shape[0]
...@@ -238,10 +246,10 @@ def findPointingbyChipID(chipID = 8, ra = 60., dec = -40.): ...@@ -238,10 +246,10 @@ def findPointingbyChipID(chipID = 8, ra = 60., dec = -40.):
r_ra = ra r_ra = ra
r_dec = dec r_dec = dec
r_rot = 0. r_rot = 0.
for i in np.arange(0,p_num,1): for i in np.arange(0, p_num, 1):
ra_n = p_list[i,0] ra_n = p_list[i, 0]
dec_n = p_list[i,1] dec_n = p_list[i, 1]
rot = p_list[i,4]*galsim.degrees rot = p_list[i, 4]*galsim.degrees
chip_wcs = pchip.getTanWCS(ra_n, dec_n, rot) chip_wcs = pchip.getTanWCS(ra_n, dec_n, rot)
c_center = pchip.getChipCenter() c_center = pchip.getChipCenter()
...@@ -257,15 +265,14 @@ def findPointingbyChipID(chipID = 8, ra = 60., dec = -40.): ...@@ -257,15 +265,14 @@ def findPointingbyChipID(chipID = 8, ra = 60., dec = -40.):
r_ra = ra_n r_ra = ra_n
r_dec = dec_n r_dec = dec_n
r_rot = rot.deg r_rot = rot.deg
if min_d == max_value: 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__": if __name__ == "__main__":
tchip, tra, tdec = 13, 60., -40. tchip, tra, tdec = 13, 60., -40.
pointing = findPointingbyChipID(chipID=tchip, ra=tra, dec=tdec) pointing = findPointingbyChipID(chipID=tchip, ra=tra, dec=tdec)
print("[ra_center, dec_center, image_rot]: ", pointing) print("[ra_center, dec_center, image_rot]: ", pointing)
tools/get_pointing_accuracy.py
View file @ 45682570
from pylab import * from pylab import *
import math, sys, numpy as np import math
import sys
import numpy as np
import astropy.coordinates as coord import astropy.coordinates as coord
from astropy.coordinates import SkyCoord from astropy.coordinates import SkyCoord
from astropy import wcs, units as u from astropy import wcs, units as u
...@@ -14,8 +16,9 @@ def transRaDec2D(ra, dec): ...@@ -14,8 +16,9 @@ def transRaDec2D(ra, dec):
z1 = np.sin(dec / 57.2957795) z1 = np.sin(dec / 57.2957795)
return np.array([x1, y1, z1]) return np.array([x1, y1, z1])
def ecl2radec(lon_ecl, lat_ecl): def ecl2radec(lon_ecl, lat_ecl):
## convert from ecliptic coordinates to equatorial coordinates # convert from ecliptic coordinates to equatorial coordinates
c_ecl = SkyCoord( c_ecl = SkyCoord(
lon=lon_ecl * u.degree, lat=lat_ecl * u.degree, frame="barycentrictrueecliptic" lon=lon_ecl * u.degree, lat=lat_ecl * u.degree, frame="barycentrictrueecliptic"
) )
...@@ -25,25 +28,26 @@ def ecl2radec(lon_ecl, lat_ecl): ...@@ -25,25 +28,26 @@ def ecl2radec(lon_ecl, lat_ecl):
def radec2ecl(ra, dec): 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_eq = SkyCoord(ra=ra * u.degree, dec=dec * u.degree, frame="icrs")
c_ecl = c_eq.transform_to("barycentrictrueecliptic") c_ecl = c_eq.transform_to("barycentrictrueecliptic")
lon_ecl, lat_ecl = c_ecl.lon.degree, c_ecl.lat.degree lon_ecl, lat_ecl = c_ecl.lon.degree, c_ecl.lat.degree
return lon_ecl, lat_ecl 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. def cal_FoVcenter_1P_equatorial(ra_FieldCenter, dec_FieldCenter, chipID=1, pa=-23.5):
### [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] is the center ra, dec of calibration fileds, such as: NEP, NGC 6397, etc.
## Calculate PA angle # [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) chip = Chip(chipID)
h_ext = ImageHeader.generateExtensionHeader( h_ext = ImageHeader.generateExtensionHeader(
chip=chip, chip=chip,
xlen=chip.npix_x, xlen=chip.npix_x,
ylen=chip.npix_y, ylen=chip.npix_y,
ra=(ra_FieldCenter * u.degree).value, ra=(ra_FieldCenter * u.degree).value,
dec=(dec_FieldCenter * u.degree).value, dec=(dec_FieldCenter * u.degree).value,
pa=pa, pa=pa,
gain=chip.gain, gain=chip.gain,
...@@ -54,7 +58,7 @@ def cal_FoVcenter_1P_equatorial(ra_FieldCenter, dec_FieldCenter, chipID = 1, pa ...@@ -54,7 +58,7 @@ def cal_FoVcenter_1P_equatorial(ra_FieldCenter, dec_FieldCenter, chipID = 1, pa
pixel_size=chip.pix_size, pixel_size=chip.pix_size,
xcen=chip.x_cen, xcen=chip.x_cen,
ycen=chip.y_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 # 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 ...@@ -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 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. def cal_FoVcenter_1P_ecliptic(lon_ecl_FieldCenter, lat_ecl_FieldCenter, chipID=1, pa=-23.5):
### [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] 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( ra_FieldCenter, dec_FieldCenter = ecl2radec(
lon_ecl_FieldCenter, lat_ecl_FieldCenter lon_ecl_FieldCenter, lat_ecl_FieldCenter
) )
## Calculate PA angle # Calculate PA angle
chip = Chip(chipID) chip = Chip(chipID)
h_ext = ImageHeader.generateExtensionHeader( h_ext = ImageHeader.generateExtensionHeader(
chip=chip, chip=chip,
xlen=chip.npix_x, xlen=chip.npix_x,
ylen=chip.npix_y, ylen=chip.npix_y,
ra=(ra_FieldCenter * u.degree).value, ra=(ra_FieldCenter * u.degree).value,
dec=(dec_FieldCenter * u.degree).value, dec=(dec_FieldCenter * u.degree).value,
pa=pa, pa=pa,
gain=chip.gain, gain=chip.gain,
...@@ -103,7 +108,7 @@ def cal_FoVcenter_1P_ecliptic(lon_ecl_FieldCenter, lat_ecl_FieldCenter, chipID = ...@@ -103,7 +108,7 @@ def cal_FoVcenter_1P_ecliptic(lon_ecl_FieldCenter, lat_ecl_FieldCenter, chipID =
pixel_size=chip.pix_size, pixel_size=chip.pix_size,
xcen=chip.x_cen, xcen=chip.x_cen,
ycen=chip.y_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 # 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 = ...@@ -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 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) chip = Chip(chipID)
h_ext = ImageHeader.generateExtensionHeader( h_ext = ImageHeader.generateExtensionHeader(
chip=chip, chip=chip,
xlen=chip.npix_x, xlen=chip.npix_x,
ylen=chip.npix_y, ylen=chip.npix_y,
ra=(p_ra * u.degree).value, ra=(p_ra * u.degree).value,
dec=(p_dec * u.degree).value, dec=(p_dec * u.degree).value,
pa=pa, pa=pa,
gain=chip.gain, gain=chip.gain,
...@@ -142,24 +148,25 @@ def getChipCenterRaDec(chipID = 1, p_ra = 60., p_dec = -40.): ...@@ -142,24 +148,25 @@ def getChipCenterRaDec(chipID = 1, p_ra = 60., p_dec = -40.):
pixel_size=chip.pix_size, pixel_size=chip.pix_size,
xcen=chip.x_cen, xcen=chip.x_cen,
ycen=chip.y_cen, ycen=chip.y_cen,
) )
wcs_h = wcs.WCS(h_ext) wcs_h = wcs.WCS(h_ext)
pixs = np.array([[4608, 4616]]) pixs = np.array([[4608, 4616]])
world_point = wcs_h.wcs_pix2world(pixs, 0) world_point = wcs_h.wcs_pix2world(pixs, 0)
RA_chip, Dec_chip = world_point[0][0], world_point[0][1] RA_chip, Dec_chip = world_point[0][0], world_point[0][1]
return RA_chip, Dec_chip return RA_chip, Dec_chip
if __name__ == '__main__': if __name__ == '__main__':
ra_input, dec_input = 270.00000, 66.56000 # NEP ra_input, dec_input = 270.00000, 66.56000 # NEP
pa = 23.5 pa = 23.5
# chipid = 2 # 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, 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)) print("chip id is %d, chip center [ra,dec] is [%f, %f], pointing center calculated [ra,dec] is [%f, %f]" % (
#for check the result chipid, ra_input, dec_input, ra, dec))
# for check the result
# testRA, testDec = getChipCenterRaDec(chipID = chipid, p_ra = ra, p_dec = dec) # testRA, testDec = getChipCenterRaDec(chipID = chipid, p_ra = ra, p_dec = dec)
# print(ra_input-testRA, dec_input-testDec) # print(ra_input-testRA, dec_input-testDec)
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment