diff --git a/observation_sim/ObservationSim.py b/observation_sim/ObservationSim.py index 1e82700de0f7e447194fcb3c8b29361b244f1b59..247636f61194ce09e32e985cc7199e088b5be058 100755 --- a/observation_sim/ObservationSim.py +++ b/observation_sim/ObservationSim.py @@ -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"]: diff --git a/observation_sim/config/header/ImageHeader.py b/observation_sim/config/header/ImageHeader.py index 7fc3444c1e48d712673ad67c2dadd3cf0e1eb6b6..25a2867df16aff0fb6d5c25942d9f757b9803371 100644 --- a/observation_sim/config/header/ImageHeader.py +++ b/observation_sim/config/header/ImageHeader.py @@ -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: diff --git a/observation_sim/mock_objects/MockObject.py b/observation_sim/mock_objects/MockObject.py index 12e36124de22019c6382b3381fefdb8394eecc65..ef660caea1f7a057fff78abe44f1202a9e3745d0 100755 --- a/observation_sim/mock_objects/MockObject.py +++ b/observation_sim/mock_objects/MockObject.py @@ -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): diff --git a/observation_sim/sim_steps/__init__.py b/observation_sim/sim_steps/__init__.py index 65a8e21b1636885b371cabe32268dab60a1d4715..e8e2661c7d92f787e27fb0216bfd13131a70b976 100644 --- a/observation_sim/sim_steps/__init__.py +++ b/observation_sim/sim_steps/__init__.py @@ -1,10 +1,13 @@ 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", +} diff --git a/observation_sim/sim_steps/add_objects.py b/observation_sim/sim_steps/add_objects.py index 9f41d3621c8457adbd66db88bac38f2ccd64d505..b50c7afcdf209622e13cd49dd7c8d38b243a79af 100644 --- a/observation_sim/sim_steps/add_objects.py +++ b/observation_sim/sim_steps/add_objects.py @@ -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 diff --git a/observation_sim/sim_steps/readout_output.py b/observation_sim/sim_steps/readout_output.py index 8ef872fd6d9348f3f46491486010fa41af935c6f..c2cb8a16ab8054a55ad7456fa01987e5d42a8d34 100644 --- a/observation_sim/sim_steps/readout_output.py +++ b/observation_sim/sim_steps/readout_output.py @@ -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) diff --git a/tools/get_pointing.py b/tools/get_pointing.py index 073d160b59321ede75351dd8c850f3a1022be287..39216a16c7c2ed7c9cb7481a799c7211288e07aa 100755 --- a/tools/get_pointing.py +++ b/tools/get_pointing.py @@ -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) diff --git a/tools/get_pointing_accurcy.py b/tools/get_pointing_accurcy.py new file mode 100644 index 0000000000000000000000000000000000000000..33d8575298c806b4631fe348b54c9bc690002260 --- /dev/null +++ b/tools/get_pointing_accurcy.py @@ -0,0 +1,165 @@ + +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) + diff --git a/tools/target_location_check.py b/tools/target_location_check.py index 34293be892b00d330c8fac94037a3502e2b6055e..11f22bac7e8eb8d74606c4e5d003ac92e5e23633 100644 --- a/tools/target_location_check.py +++ b/tools/target_location_check.py @@ -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(