Commit ea03057d authored by Fang Yuedong's avatar Fang Yuedong
Browse files

Merge branch 'master' into 'current_stable_for_tests'

Master

See merge request !25
parents b6519ab0 d0c206b0
......@@ -92,16 +92,21 @@ class Observation(object):
input_time_str=time_str
)
ra_cen, dec_cen = ra_cen[0], dec_cen[0]
ra_offset, dec_offset = pointing.ra - ra_cen, pointing.dec - dec_cen
else:
ra_cen = pointing.ra
dec_cen = pointing.dec
ra_offset, dec_offset = 0., 0.
# Prepare necessary chip properties for simulation
chip = self.prepare_chip_for_exposure(chip, ra_cen, dec_cen, pointing)
# Initialize SimSteps
sim_steps = SimSteps(overall_config=self.config,
chip_output=chip_output, all_filters=self.all_filters)
chip_output=chip_output,
all_filters=self.all_filters,
ra_offset=ra_offset,
dec_offset=dec_offset)
for step in pointing.obs_param["call_sequence"]:
if self.config["run_option"]["out_cat_only"]:
......
......@@ -434,7 +434,7 @@ def generatePrimaryHeader(xlen=9216, ylen=9232, pointing_id='00000001', pointing
co = coord.SkyCoord(ra, dec, unit='deg')
ra_hms = format(co.ra.hms.h, '02.0f') + format(co.ra.hms.m,
'02.0f') + format(co.ra.hms.s, '02.1f')
'02.0f') + format(co.ra.hms.s, '04.1f')
dec_hms = format(co.dec.dms.d, '02.0f') + format(abs(co.dec.dms.m),
'02.0f') + format(abs(co.dec.dms.s), '02.0f')
if dec >= 0:
......
......@@ -59,13 +59,13 @@ class MockObject(object):
flux = self.getFluxFilter(filt)
return flux * tel.pupil_area * exptime
def getPosWorld(self):
ra = self.param["ra"]
dec = self.param["dec"]
def getPosWorld(self, ra_offset=0., dec_offset=0.):
ra = self.param["ra"] + ra_offset
dec = self.param["dec"] + dec_offset
return galsim.CelestialCoord(ra=ra * galsim.degrees, dec=dec * galsim.degrees)
def getPosImg_Offset_WCS(self, img, fdmodel=None, chip=None, verbose=True, chip_wcs=None, img_header=None):
self.posImg = img.wcs.toImage(self.getPosWorld())
def getPosImg_Offset_WCS(self, img, fdmodel=None, chip=None, verbose=True, chip_wcs=None, img_header=None, ra_offset=0., dec_offset=0.):
self.posImg = img.wcs.toImage(self.getPosWorld(ra_offset, dec_offset))
self.localWCS = img.wcs.local(self.posImg)
# Apply field distortion model
if (fdmodel is not None) and (chip is not None):
......
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",
"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
......
......@@ -86,10 +86,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)
......
......@@ -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
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
# 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