旋转角测量a

懒得起名字22

已于 2022-02-16 10:43:30 修改

阅读量227

点赞数

CC 4.0 BY-SA版权

文章标签： python 算法

于 2021-06-22 14:57:41 首次发布

本文链接：https://blog.csdn.net/qq_40692494/article/details/118107299

from astropy.io import fits
from astropy.wcs import WCS
import numpy as np
from astropy import units as u
from astropy.coordinates import SkyCoord
from astropy.utils.data import get_pkg_data_filename
import math

fn = get_pkg_data_filename('tmp/TS-990-500-D1_21_1-offset_1024-w_4096-h_4096-Mono16.fits')

f = fits.open(fn)
w = WCS(f[0].header)
c1 = w.pixel_to_world(1, 1)
c2 = w.pixel_to_world(4095, 4095)

c1 = SkyCoord(c1, frame='fk5')
ra1 = c1.ra.degree
print(ra1)
dec1 = c1.dec.degree
print(dec1)

k=math.pi/180.0

c2 = SkyCoord(c2, frame='fk5')
ra2 = c2.ra.degree
# print(ra2)
dec2 = c2.dec.degree
# print(dec2)
tan1 = (dec2-dec1)/((ra2-ra1)*math.cos((dec2+dec1)*k/2))
theta1 = math.atan(tan1)*180/np.pi


fn2 = get_pkg_data_filename('tmp/TS-1015-500_18_1-offset_1024-w_4096-h_4096-Mono16.fits')

f2 = fits.open(fn2)
w2 = WCS(f2[0].header)
c3 = w2.pixel_to_world(1, 1)
c4 = w2.pixel_to_world(4095, 4095)

c3 = SkyCoord(c3)
ra3 = c3.ra.degree
dec3 = c3.dec.degree

c4 = SkyCoord(c4)
ra4 = c4.ra.degree
dec4 = c4.dec.degree

tan2 = (dec4-dec3)/((ra4-ra3)*math.cos((dec2+dec1)*k/2))
theta2 = math.atan(tan2)*180/np.pi

deta = theta2-theta1
print(deta)