Commit 638b95d7 authored by Fang Yuedong's avatar Fang Yuedong
Browse files

Merge branch 'release_v3.0' into 'master'

# Conflicts:
#   observation_sim/mock_objects/Galaxy.py
#   observation_sim/mock_objects/MockObject.py
#   observation_sim/psf/PSFInterpSLS.py
parents 570a35a4 d8118d46
Pipeline #6547 failed with stage
in 0 seconds
...@@ -5,4 +5,3 @@ from .Quasar import Quasar ...@@ -5,4 +5,3 @@ from .Quasar import Quasar
from .Star import Star from .Star import Star
from .Stamp import Stamp from .Stamp import Stamp
from .FlatLED import FlatLED from .FlatLED import FlatLED
from .ExtinctionMW import ExtinctionMW
...@@ -17,7 +17,7 @@ class PSFGauss(PSFModel): ...@@ -17,7 +17,7 @@ class PSFGauss(PSFModel):
self.fwhm = self.fwhmGauss(r=psfRa) self.fwhm = self.fwhmGauss(r=psfRa)
self.sigGauss = psfRa # 80% light radius self.sigGauss = psfRa # 80% light radius
self.sigSpin = sigSpin self.sigSpin = sigSpin
self.psf = galsim.Gaussian(flux=1.0, fwhm=self.fwhm) self.psf = galsim.Gaussian(flux=1.0, fwhm=fwhm)
def perfGauss(self, r, sig): def perfGauss(self, r, sig):
""" """
......
...@@ -435,16 +435,7 @@ class PSFInterpSLS(PSFModel): ...@@ -435,16 +435,7 @@ 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
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 # from astropy.io import fits
# fits.writeto(str(bandNo) + '_' + g_order+ '_psf_o.fits', PSF_int_trans) # fits.writeto(str(bandNo) + '_' + g_order+ '_psf_o.fits', PSF_int_trans)
......
import os import os
class SimSteps: class SimSteps:
def __init__(self, overall_config, chip_output, all_filters, ra_offset=0., dec_offset=0.): def __init__(self, overall_config, chip_output, all_filters):
self.overall_config = overall_config self.overall_config = overall_config
self.chip_output = chip_output self.chip_output = chip_output
self.all_filters = all_filters self.all_filters = all_filters
self.ra_offset = ra_offset
self.dec_offset = dec_offset
from .prepare_headers import prepare_headers, updateHeaderInfo from .prepare_headers import prepare_headers, updateHeaderInfo
from .add_sky_background import add_sky_background_sci, add_sky_flat_calibration, add_sky_background from .add_sky_background import add_sky_background_sci, add_sky_flat_calibration, add_sky_background
...@@ -18,7 +15,6 @@ class SimSteps: ...@@ -18,7 +15,6 @@ class SimSteps:
from .readout_output import add_prescan_overscan, add_readout_noise, apply_gain, quantization_and_output from .readout_output import add_prescan_overscan, add_readout_noise, apply_gain, quantization_and_output
from .add_LED_flat import add_LED_Flat from .add_LED_flat import add_LED_Flat
SIM_STEP_TYPES = { SIM_STEP_TYPES = {
"scie_obs": "add_objects", "scie_obs": "add_objects",
"sky_background": "add_sky_background", "sky_background": "add_sky_background",
...@@ -35,6 +31,6 @@ SIM_STEP_TYPES = { ...@@ -35,6 +31,6 @@ SIM_STEP_TYPES = {
"readout_noise": "add_readout_noise", "readout_noise": "add_readout_noise",
"gain": "apply_gain", "gain": "apply_gain",
"quantization_and_output": "quantization_and_output", "quantization_and_output": "quantization_and_output",
"led_calib_model": "add_LED_Flat", "led_calib_model":"add_LED_Flat",
"sky_flatField": "add_sky_flat_calibration", "sky_flatField":"add_sky_flat_calibration",
} }
\ No newline at end of file
...@@ -145,7 +145,7 @@ def add_objects(self, chip, filt, tel, pointing, catalog, obs_param): ...@@ -145,7 +145,7 @@ def add_objects(self, chip, filt, tel, pointing, catalog, obs_param):
# Get position of object on the focal plane # Get position of object on the focal plane
pos_img, _, _, _, fd_shear = obj.getPosImg_Offset_WCS( 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, ra_offset=self.ra_offset, dec_offset=self.dec_offset) img=chip.img, fdmodel=fd_model, chip=chip, verbose=False, chip_wcs=chip_wcs, img_header=self.h_ext)
# [TODO] For now, only consider objects which their centers (after field distortion) are projected within the focal plane # [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 # Otherwise they will be considered missed objects
...@@ -194,8 +194,7 @@ def add_objects(self, chip, filt, tel, pointing, catalog, obs_param): ...@@ -194,8 +194,7 @@ def add_objects(self, chip, filt, tel, pointing, catalog, obs_param):
if isUpdated == 1: if isUpdated == 1:
# TODO: add up stats # TODO: add up stats
self.chip_output.cat_add_obj( self.chip_output.cat_add_obj(obj, pos_img, pos_shear)
obj, pos_img, pos_shear, ra_offset=self.ra_offset, dec_offset=self.dec_offset)
pass pass
elif isUpdated == 0: elif isUpdated == 0:
missed_obj += 1 missed_obj += 1
...@@ -217,27 +216,7 @@ def add_objects(self, chip, filt, tel, pointing, catalog, obs_param): ...@@ -217,27 +216,7 @@ def add_objects(self, chip, filt, tel, pointing, catalog, obs_param):
# Unload SED: # Unload SED:
obj.unload_SED() obj.unload_SED()
del obj 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 del psf_model
gc.collect() gc.collect()
......
...@@ -27,7 +27,7 @@ def add_poisson_and_dark(self, chip, filt, tel, pointing, catalog, obs_param): ...@@ -27,7 +27,7 @@ def add_poisson_and_dark(self, chip, filt, tel, pointing, catalog, obs_param):
InputDark=None) InputDark=None)
else: else:
chip.img, _ = chip_utils.add_poisson(img=chip.img, chip.img, _ = chip_utils.add_poisson(img=chip.img,
chip=chip, chip=self,
exptime=exptime, exptime=exptime,
poisson_noise=chip.poisson_noise, poisson_noise=chip.poisson_noise,
dark_noise=0.) dark_noise=0.)
...@@ -81,5 +81,5 @@ def add_bias(self, chip, filt, tel, pointing, catalog, obs_param): ...@@ -81,5 +81,5 @@ def add_bias(self, chip, filt, tel, pointing, catalog, obs_param):
nsecx=chip.nsecx, nsecx=chip.nsecx,
seed=self.overall_config["random_seeds"]["seed_biasNonUniform"]+chip.chipID) seed=self.overall_config["random_seeds"]["seed_biasNonUniform"]+chip.chipID)
elif obs_param["bias_16channel"] == False: elif obs_param["bias_16channel"] == False:
chip.img += chip.bias_level chip.img += self.bias_level
return chip, filt, tel, pointing return chip, filt, tel, pointing
...@@ -87,10 +87,10 @@ def quantization_and_output(self, chip, filt, tel, pointing, catalog, obs_param) ...@@ -87,10 +87,10 @@ def quantization_and_output(self, chip, filt, tel, pointing, catalog, obs_param)
fname = os.path.join(self.chip_output.subdir, fname = os.path.join(self.chip_output.subdir,
self.h_prim['FILENAME'] + '.fits') self.h_prim['FILENAME'] + '.fits')
# f_name_size = 68 f_name_size = 68
# if (len(self.h_prim['FILENAME']) > f_name_size): if (len(self.h_prim['FILENAME']) > f_name_size):
# self.updateHeaderInfo(header_flag='prim', keys=['FILENAME'], values=[ self.updateHeaderInfo(header_flag='prim', keys=['FILENAME'], values=[
# self.h_prim['FILENAME'][0:f_name_size]]) self.h_prim['FILENAME'][0:f_name_size]])
hdu1 = fits.PrimaryHDU(header=self.h_prim) hdu1 = fits.PrimaryHDU(header=self.h_prim)
......
#!/bin/bash #!/bin/bash
#SBATCH -J CSSTSim #SBATCH -J CSSTSim
#SBATCH -N 2 #SBATCH -N 1
#SBATCH --ntasks-per-node=6 #SBATCH --ntasks-per-node=6
#SBATCH -p debug #SBATCH -p debug
#SBATCH --mem=96G #SBATCH --mem=60G
module load mpi/hpcx/2.4.1/gcc-7.3.1 module load mpi/hpcx/2.4.1/gcc-7.3.1
date date
...@@ -12,4 +12,4 @@ date ...@@ -12,4 +12,4 @@ date
#限定单节点任务数 #限定单节点任务数
srun hostname -s | sort -n | awk -F"-" '{print $2}' | uniq > pnodes 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 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 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 \ No newline at end of file
...@@ -63,15 +63,10 @@ class Chip(object): ...@@ -63,15 +63,10 @@ class Chip(object):
ycen = self.cen_pix_y ycen = self.cen_pix_y
if pix_scale == None: if pix_scale == None:
pix_scale = self.pix_scale 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 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 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
...@@ -144,11 +139,12 @@ def getobsPA(ra, dec): ...@@ -144,11 +139,12 @@ def getobsPA(ra, dec):
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
if (ra<90 or ra> 270): # if (ra>90 and 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')
...@@ -265,7 +261,7 @@ def findPointingbyChipID(chipID = 8, ra = 60., dec = -40.): ...@@ -265,7 +261,7 @@ def findPointingbyChipID(chipID = 8, ra = 60., dec = -40.):
if __name__ == "__main__": if __name__ == "__main__":
tchip, tra, tdec = 13, 60., -40. tchip, tra, tdec = 8, 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)
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,15 +119,10 @@ def getTanWCS(ra, dec, img_rot, pix_scale=0.074): ...@@ -119,15 +119,10 @@ def getTanWCS(ra, dec, img_rot, pix_scale=0.074):
""" """
xcen, ycen = 0, 0 xcen, ycen = 0, 0
img_rot = img_rot * galsim.degrees 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 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 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) moscen = galsim.PositionD(x=xcen, y=ycen)
sky_center = galsim.CelestialCoord( 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