Commit 1b4f4012 authored by Fang Yuedong's avatar Fang Yuedong
Browse files

Merge branch 'master' into 'release_v3.0'

version v3.1.0

See merge request !31
parents d8118d46 5544935f
......@@ -20,7 +20,30 @@ import cython
cdef extern from "math.h":
double sqrt(double x)
double exp(double x)
def check_nan2D(np.ndarray[FTYPE_t, ndim=2] arr):
cdef int i, j
cdef int nrows = arr.shape[0]
cdef int ncols = arr.shape[1]
# 遍历数组的每个元素并检查是否存在 NaN
for i in range(nrows):
for j in range(ncols):
if np.isnan(arr[i, j]) | np.isinf(arr[i, j]):
return True
return False
def check_nan1d(np.ndarray[DTYPE_t, ndim=1] arr):
cdef int i
cdef int n = arr.shape[0]
# 遍历数组的每个元素并检查是否存在 NaN
for i in range(n):
if np.isnan(arr[i]) | np.isinf(arr[i]):
return True
return False
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.embedsignature(True)
......@@ -53,6 +76,18 @@ def disperse_grism_object(np.ndarray[FTYPE_t, ndim=2] flam,
nk = len(idxl)
nl = len(full)
#if check_nan2D(flam):
# print("DEBUG: disperse, input Array 'flam' contains NaN.")
#if check_nan1d(ysens):
# print("DEBUG: disperse, input Array 'ysens' contains NaN.")
#if check_nan1d(yfrac):
# print("DEBUG: disperse, input Array 'yfrac' contains NaN.")
#if check_nan1d(full):
# print("DEBUG: disperse, input Array 'full' contains NaN before processing.")
if (flat is not None):
nlamb = len(wlambda)
......@@ -95,14 +130,15 @@ def disperse_grism_object(np.ndarray[FTYPE_t, ndim=2] flam,
else:
for i in range(0-x0[1], x0[1]):
if (x0[1]+i < 0) | (x0[1]+i >= shd[1]):
x_pos = x0[1]+i
if (x_pos < 0) | (x_pos >= shd[1]):
continue
for j in range(0-x0[0], x0[0]):
if (x0[0]+j < 0) | (x0[0]+j >= shd[0]):
y_pos = x0[0]+j
if (y_pos < 0) | (y_pos >= shd[0]):
continue
fl_ij = flam[x0[0]+j, x0[1]+i] #/1.e-17
fl_ij = flam[y_pos, x_pos] #/1.e-17
if (fl_ij == 0):
continue
......@@ -110,10 +146,13 @@ def disperse_grism_object(np.ndarray[FTYPE_t, ndim=2] flam,
k1 = idxl[k]+j*shg[1]+i
if (k1 >= 0) & (k1 < nl):
full[k1] += ysens[k]*fl_ij*(1-yfrac[k])
k2 = idxl[k]+(j+1)*shg[1]+i
if (k2 >= 0) & (k2 < nl):
full[k2] += ysens[k]*fl_ij*yfrac[k]
#if (check_nan1d(full)):
# print("DEBUG: disperse, output Array 'full' contains NaN after processing.+++++++++++++++++++++++++++")
return True
......
......@@ -5,3 +5,4 @@ from .Quasar import Quasar
from .Star import Star
from .Stamp import Stamp
from .FlatLED import FlatLED
from .ExtinctionMW import ExtinctionMW
......@@ -17,7 +17,7 @@ class PSFGauss(PSFModel):
self.fwhm = self.fwhmGauss(r=psfRa)
self.sigGauss = psfRa # 80% light radius
self.sigSpin = sigSpin
self.psf = galsim.Gaussian(flux=1.0, fwhm=fwhm)
self.psf = galsim.Gaussian(flux=1.0, fwhm=self.fwhm)
def perfGauss(self, r, sig):
"""
......@@ -122,4 +122,4 @@ class PSFGauss(PSFModel):
# return ell, beta, qr
PSFshear = galsim.Shear(e=ell, beta=beta*galsim.radians)
return self.psf.shear(PSFshear), PSFshear
return self.psf.shear(PSFshear), PSFshear
\ No newline at end of file
......@@ -20,8 +20,10 @@ import os
from astropy.io import fits
from astropy.modeling.models import Gaussian2D
from scipy import signal
from scipy import signal, interpolate
import datetime
import gc
# from jax import numpy as jnp
LOG_DEBUG = False # ***#
NPSF = 900 # ***# 30*30
......@@ -433,7 +435,16 @@ 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
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))
####
# from astropy.io import fits
# fits.writeto(str(bandNo) + '_' + g_order+ '_psf_o.fits', PSF_int_trans)
......@@ -479,6 +490,215 @@ class PSFInterpSLS(PSFModel):
return PSF_int_trans, PSF_int
def get_PSF_AND_convolve_withsubImg(self, chip, cutImg=None, pos_img_local=[1000, 1000], bandNo=1, g_order='A', grating_split_pos=3685):
"""
Get the PSF at a given image position
Parameters:
chip: A 'Chip' object representing the chip we want to extract PSF from.
pos_img: A 'galsim.Position' object representing the image position.
bandpass: A 'galsim.Bandpass' object representing the wavelength range.
pixSize: The pixels size of psf matrix
findNeighMode: 'treeFind' or 'hoclistFind'
Returns:
PSF: A 'galsim.GSObject'.
"""
order_IDs = {'A': '1', 'B': '0', 'C': '0', 'D': '0', 'E': '0'}
contam_order_sigma = {'C': 0.28032344707964174,
'D': 0.39900182912061344, 'E': 1.1988309797685412} # arcsec
x_start = chip.x_cen/chip.pix_size - chip.npix_x / 2.
y_start = chip.y_cen/chip.pix_size - chip.npix_y / 2.
# print(pos_img.x - x_start)
# pos_img_x = pos_img_local[0] + x_start
# pos_img_y = pos_img_local[1] + y_start
# pos_img = galsim.PositionD(pos_img_x, pos_img_y)
# centerPos_local = cutImg.ncol/2.
if pos_img_local[0] < grating_split_pos:
psf_data = self.grating1_data
else:
psf_data = self.grating2_data
grating_order = order_IDs[g_order]
# if grating_order in ['-2','-1','2']:
# grating_order = '1'
# if grating_order in ['0', '1']:
psf_order = psf_data['order'+grating_order]
psf_order_b = psf_order['band'+str(bandNo)]
psf_b_dat = psf_order_b['band_data']
# pos_p = psf_b_dat[1].data
pos_p = psf_b_dat[1].data/chip.pix_size - np.array([y_start, x_start])
pc_coeff = psf_b_dat[2].data
pcs = psf_b_dat[0].data
npc = 10
m_size = int(pcs.shape[0]**0.5)
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)
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)
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)
bounds = cutImg.bounds & b_img.bounds
if bounds.area() == 0:
continue
# ys = cutImg.ymin
# if ys < 0:
# ys = 0
# ye = cutImg.ymin+cutImg.nrow
# if ye >= cy_len-1:
# ye = cy_len-1
# if ye - ys <=0:
# continue
# xs = cutImg.xmin
# if xs < 0:
# xs = 0
# xe = cutImg.xmin+cutImg.ncol
# if xe >= cx_len-1:
# xe = cx_len-1
# if xe - xs <=0:
# continue
ys = bounds.ymin
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)
# t2=datetime.datetime.now()
# print("time interpolate:", t2-t1)
# if U.shape != cutImg.array.shape:
# print('DEBUG:SHAPE',cutImg.ncol,cutImg.nrow,cutImg.xmin, cutImg.ymin)
# continue
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)
# t3=datetime.datetime.now()
# print("time convole:", t3-t2)
del U
del img_tmp
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_L2 = ['order-2','order-1','order0','order1','order2']
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):
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']:
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])
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)
# temp_sort_dist[:, 1] = dist2
# # print(temp_sort_dist)
# dits2_sortlist = sorted(temp_sort_dist, key=lambda x: x[1])
# # print(dits2_sortlist)
# nearest4p = np.zeros([4, 3])
# pc_coeff_4p = np.zeros([npca, 4])
# for i in np.arange(4):
# smaller_ids = int(dits2_sortlist[i][0])
# nearest4p[i, 0] = pos_p[smaller_ids, 1]
# nearest4p[i, 1] = pos_p[smaller_ids, 0]
# # print(pos_p[smaller_ids, 1],pos_p[smaller_ids, 0])
# nearest4p[i, 2] = dits2_sortlist[i][1]
# pc_coeff_4p[:, i] = pc_coeff[npca, smaller_ids]
# # idw_dist = 1/(np.sqrt((px-nearest4p[:, 0]) * (px-nearest4p[:, 0]) + (
# # py-nearest4p[:, 1]) * (py-nearest4p[:, 1])))
# idw_dist = 1/(np.sqrt(nearest4p[:, 2]))
# coeff_int = np.zeros(npca)
# for i in np.arange(4):
# coeff_int = coeff_int + pc_coeff_4p[:, i]*idw_dist[i]
# coeff_mat[:, m, n] = coeff_int
m_size = int(pcs.shape[0]**0.5)
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)
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)
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)
U1 = interpolate.griddata(X_, Z_, (M, N),
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()
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)
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()
# pixSize = np.rad2deg(self.pixsize*1e-3/28)*3600 #set psf pixsize
#
# # assert self.iccd == int(chip.getChipLabel(chipID=chip.chipID)), 'ERROR: self.iccd != chip.chipID'
......
import os
class SimSteps:
def __init__(self, overall_config, chip_output, all_filters):
def __init__(self, overall_config, chip_output, all_filters, ra_offset=0., dec_offset=0.):
self.overall_config = overall_config
self.chip_output = chip_output
self.all_filters = all_filters
self.ra_offset = ra_offset
self.dec_offset = dec_offset
from .prepare_headers import prepare_headers, updateHeaderInfo
from .add_sky_background import add_sky_background_sci, add_sky_flat_calibration, add_sky_background
......@@ -15,6 +18,7 @@ class SimSteps:
from .readout_output import add_prescan_overscan, add_readout_noise, apply_gain, quantization_and_output
from .add_LED_flat import add_LED_Flat
SIM_STEP_TYPES = {
"scie_obs": "add_objects",
"sky_background": "add_sky_background",
......@@ -31,6 +35,6 @@ SIM_STEP_TYPES = {
"readout_noise": "add_readout_noise",
"gain": "apply_gain",
"quantization_and_output": "quantization_and_output",
"led_calib_model":"add_LED_Flat",
"sky_flatField":"add_sky_flat_calibration",
}
\ No newline at end of file
"led_calib_model": "add_LED_Flat",
"sky_flatField": "add_sky_flat_calibration",
}
......@@ -145,7 +145,7 @@ def add_objects(self, chip, filt, tel, pointing, catalog, obs_param):
# Get position of object on the focal plane
pos_img, _, _, _, fd_shear = obj.getPosImg_Offset_WCS(
img=chip.img, fdmodel=fd_model, chip=chip, verbose=False, chip_wcs=chip_wcs, img_header=self.h_ext)
img=chip.img, fdmodel=fd_model, chip=chip, verbose=False, chip_wcs=chip_wcs, img_header=self.h_ext, ra_offset=self.ra_offset, dec_offset=self.dec_offset)
# [TODO] For now, only consider objects which their centers (after field distortion) are projected within the focal plane
# Otherwise they will be considered missed objects
......@@ -194,7 +194,8 @@ def add_objects(self, chip, filt, tel, pointing, catalog, obs_param):
if isUpdated == 1:
# TODO: add up stats
self.chip_output.cat_add_obj(obj, pos_img, pos_shear)
self.chip_output.cat_add_obj(
obj, pos_img, pos_shear, ra_offset=self.ra_offset, dec_offset=self.dec_offset)
pass
elif isUpdated == 0:
missed_obj += 1
......@@ -216,7 +217,27 @@ def add_objects(self, chip, filt, tel, pointing, catalog, obs_param):
# Unload SED:
obj.unload_SED()
del obj
gc.collect()
# gc.collect()
if chip.survey_type == "spectroscopic" and not self.overall_config["run_option"]["out_cat_only"] and chip.slsPSFOptim:
# from observation_sim.instruments.chip import chip_utils as chip_utils
# gn = chip_utils.getChipSLSGratingID(chip.chipID)[0]
# img1 = np.zeros([2,chip.img.array.shape[0],chip.img.array.shape[1]])
# for id1 in np.arange(2):
# gn = chip_utils.getChipSLSGratingID(chip.chipID)[id1]
# img_i = 0
# for id2 in ['0','1']:
# o_n = "order"+id2
# for id3 in ['1','2','3','4']:
# w_n = "w"+id3
# img1[img_i] = img1[img_i] + chip.img_stack[gn][o_n][w_n].array
# img_i = img_i + 1
# from astropy.io import fits
# fits.writeto('order0.fits',img1[0],overwrite=True)
# fits.writeto('order1.fits',img1[1],overwrite=True)
psf_model.convolveFullImgWithPCAPSF(chip)
del psf_model
gc.collect()
......
......@@ -27,7 +27,7 @@ def add_poisson_and_dark(self, chip, filt, tel, pointing, catalog, obs_param):
InputDark=None)
else:
chip.img, _ = chip_utils.add_poisson(img=chip.img,
chip=self,
chip=chip,
exptime=exptime,
poisson_noise=chip.poisson_noise,
dark_noise=0.)
......@@ -81,5 +81,5 @@ def add_bias(self, chip, filt, tel, pointing, catalog, obs_param):
nsecx=chip.nsecx,
seed=self.overall_config["random_seeds"]["seed_biasNonUniform"]+chip.chipID)
elif obs_param["bias_16channel"] == False:
chip.img += self.bias_level
chip.img += chip.bias_level
return chip, filt, tel, pointing
......@@ -44,6 +44,7 @@ def apply_gain(self, chip, filt, tel, pointing, catalog, obs_param):
seed=self.overall_config["random_seeds"]["seed_gainNonUniform"]+chip.chipID)
elif obs_param["gain_16channel"] == False:
chip.img /= chip.gain
chip.gain_channel = np.ones(chip.nsecy*chip.nsecx)*chip.gain
return chip, filt, tel, pointing
......@@ -86,10 +87,10 @@ def quantization_and_output(self, chip, filt, tel, pointing, catalog, obs_param)
fname = os.path.join(self.chip_output.subdir,
self.h_prim['FILENAME'] + '.fits')
f_name_size = 68
if (len(self.h_prim['FILENAME']) > f_name_size):
self.updateHeaderInfo(header_flag='prim', keys=['FILENAME'], values=[
self.h_prim['FILENAME'][0:f_name_size]])
# f_name_size = 68
# if (len(self.h_prim['FILENAME']) > f_name_size):
# self.updateHeaderInfo(header_flag='prim', keys=['FILENAME'], values=[
# self.h_prim['FILENAME'][0:f_name_size]])
hdu1 = fits.PrimaryHDU(header=self.h_prim)
......
......@@ -76,7 +76,7 @@ with open("requirements.txt", "r") as f:
]
setup(name='csst_msc_sim',
version='3.0.0',
version='3.1.0',
packages=find_packages(),
# install_requires=[
# # 'numpy>=1.18.5',
......
#!/bin/bash
#SBATCH -J CSSTSim
#SBATCH -N 1
#SBATCH -N 2
#SBATCH --ntasks-per-node=6
#SBATCH -p debug
#SBATCH --mem=60G
#SBATCH --mem=96G
module load mpi/hpcx/2.4.1/gcc-7.3.1
date
......@@ -12,4 +12,4 @@ date
#限定单节点任务数
srun hostname -s | sort -n | awk -F"-" '{print $2}' | uniq > pnodes
mpirun -mca pml ucx -x UCX_NET_DEVICES=mlx5_0:1 -machinefile pnodes -np 6 --map-by node python3 /public/home/fangyuedong/project/csst_msc_sim/run_sim.py --config_file config_overall.yaml --catalog C9_Catalog -c /public/home/fangyuedong/project/csst_msc_sim/config
\ No newline at end of file
mpirun -mca pml ucx -x UCX_NET_DEVICES=mlx5_0:1 -machinefile pnodes -np 12 --map-by node python3 /public/home/fangyuedong/project/csst_msc_sim/run_sim.py --config_file config_overall.yaml --catalog C9_Catalog -c /public/home/fangyuedong/project/csst_msc_sim/config
\ No newline at end of file
......@@ -63,10 +63,15 @@ class Chip(object):
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.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
......@@ -139,12 +144,11 @@ def getobsPA(ra, dec):
angle = math.acos(np.dot(l1l2cross,pdl2cross)/(np.linalg.norm(l1l2cross)*np.linalg.norm(pdl2cross)))
angle = (angle)/math.pi*180
# if (ra>90 and ra< 270):
# angle=-angle
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')
......@@ -261,7 +265,7 @@ def findPointingbyChipID(chipID = 8, ra = 60., dec = -40.):
if __name__ == "__main__":
tchip, tra, tdec = 8, 60., -40.
tchip, tra, tdec = 13, 60., -40.
pointing = findPointingbyChipID(chipID=tchip, ra=tra, dec=tdec)
print("[ra_center, dec_center, image_rot]: ", pointing)
from pylab import *
import math, sys, numpy as np
import astropy.coordinates as coord
from astropy.coordinates import SkyCoord
from astropy import wcs, units as u
from observation_sim.config.header import ImageHeader
from observation_sim.instruments import Telescope, Chip, FilterParam, Filter, FocalPlane
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 ecl2radec(lon_ecl, lat_ecl):
## convert from ecliptic coordinates to equatorial coordinates
c_ecl = SkyCoord(
lon=lon_ecl * u.degree, lat=lat_ecl * u.degree, frame="barycentrictrueecliptic"
)
c_eq = c_ecl.transform_to("icrs")
ra_eq, dec_eq = c_eq.ra.degree, c_eq.dec.degree
return ra_eq, dec_eq
def radec2ecl(ra, dec):
## 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
chip = Chip(chipID)
h_ext = ImageHeader.generateExtensionHeader(
chip=chip,
xlen=chip.npix_x,
ylen=chip.npix_y,
ra=(ra_FieldCenter * u.degree).value,
dec=(dec_FieldCenter * u.degree).value,
pa=pa,
gain=chip.gain,
readout=chip.read_noise,
dark=chip.dark_noise,
saturation=90000,
pixel_scale=chip.pix_scale,
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
wcs_h = wcs.WCS(h_ext)
pixs = np.array([[4608, 4616]])
world_point = wcs_h.wcs_pix2world(pixs, 0)
ra_ChipCenter, dec_ChipCenter = world_point[0][0], world_point[0][1]
d_ra = ra_FieldCenter - ra_ChipCenter
d_dec = dec_FieldCenter - dec_ChipCenter
ra_PointCenter = ra_FieldCenter + d_ra
dec_PointCenter = dec_FieldCenter + d_dec
lon_ecl_PointCenter, lat_ecl_PointCenter = radec2ecl(
ra_PointCenter, dec_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.
### [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
chip = Chip(chipID)
h_ext = ImageHeader.generateExtensionHeader(
chip=chip,
xlen=chip.npix_x,
ylen=chip.npix_y,
ra=(ra_FieldCenter * u.degree).value,
dec=(dec_FieldCenter * u.degree).value,
pa=pa,
gain=chip.gain,
readout=chip.read_noise,
dark=chip.dark_noise,
saturation=90000,
pixel_scale=chip.pix_scale,
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
wcs_h = wcs.WCS(h_ext)
pixs = np.array([[4608, 4616]])
world_point = wcs_h.wcs_pix2world(pixs, 0)
ra_ChipCenter, dec_ChipCenter = world_point[0][0], world_point[0][1]
d_ra = ra_FieldCenter - ra_ChipCenter
d_dec = dec_FieldCenter - dec_ChipCenter
ra_PointCenter = ra_FieldCenter + d_ra
dec_PointCenter = dec_FieldCenter + d_dec
lon_ecl_PointCenter, lat_ecl_PointCenter = radec2ecl(
ra_PointCenter, dec_PointCenter
)
return ra_PointCenter, dec_PointCenter, lon_ecl_PointCenter, lat_ecl_PointCenter
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,
dec=(p_dec * u.degree).value,
pa=pa,
gain=chip.gain,
readout=chip.read_noise,
dark=chip.dark_noise,
saturation=90000,
pixel_scale=chip.pix_scale,
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):
ra, dec, lon_ecl, lat_ecl = cal_FoVcenter_1P_equatorial(
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
# testRA, testDec = getChipCenterRaDec(chipID = chipid, p_ra = ra, p_dec = dec)
# print(ra_input-testRA, dec_input-testDec)
......@@ -119,10 +119,15 @@ def getTanWCS(ra, dec, img_rot, pix_scale=0.074):
"""
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
dudx = -np.cos(img_rot.rad) * pix_scale
dudy = -np.sin(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
dvdy = -np.cos(img_rot.rad) * pix_scale
moscen = galsim.PositionD(x=xcen, y=ycen)
sky_center = galsim.CelestialCoord(
......
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