get_pointing_accuracy.py 5.83 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
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)