0

我正在尝试实现一个非常简单的天体测量代码。我手动找到了图片中几颗星星的坐标（RA/DEC 和 x/y 像素）。一切似乎都直截了当，但我仍然得到奇怪的结果，相差几度。

我正在尝试解决我拍摄的 CCD 图像的板常数，在那里我手动找到了图片中的星星坐标和位置，现在我想尝试找到 (0,0) 点的真实世界坐标。

我希望有人可以帮助我编写代码，或者可以告诉我如何正确执行。提前非常感谢！

这是我的代码：

import numpy as np
import os

def astrometry(star_pos, xpix, ypix, focallength, target_RA, target_DEC):
    pi = np.pi
    DegToRad = pi / 180
    RadToDeg = 180 / pi
    n = len(star_pos)
    (target_RA, target_DEC) = (target_RA, target_DEC)
    print(target_RA, target_DEC)
    # 1) Obtain star coordinates in pixel and RA/DEC
    x_pix = [row[0] for row in star_pos]
    y_pix = [row[1] for row in star_pos]

    ra_star = [row[2] for row in star_pos]
    dec_star = [row[3] for row in star_pos]

    # 2) Calculate the standard coordinates of the stars
    X_star = np.zeros(n)
    Y_star = np.zeros(n)
    for i in range(n):
        X_star[i] = -(np.cos(DegToRad*dec_star[i])*np.sin(DegToRad*(ra_star[i] - target_RA)))/(np.cos(DegToRad*target_DEC)*np.cos(DegToRad*dec_star[i])*np.cos(DegToRad*(ra_star[i]-target_RA)) + np.sin(DegToRad*target_DEC)*np.sin(DegToRad*dec_star[i]))
        Y_star[i] = -(np.sin(DegToRad*target_DEC)*np.cos(DegToRad*dec_star[i])*np.cos(DegToRad*(ra_star[i]-target_RA)) - np.cos(DegToRad*target_DEC)*np.sin(DegToRad*dec_star[i]))/(np.cos(DegToRad*target_DEC)*np.cos(DegToRad*dec_star[i])*np.cos(DegToRad*(ra_star[i]-target_RA)) + np.sin(DegToRad*target_DEC)*np.sin(DegToRad*dec_star[i]))

    # 3) Calculate the plate constants (Check my notes)
    def calc_plate_const(k,x,y,X):
        c_down = ((x[k+1]-x[k])*(y[k]*x[k+2]-y[k+2]*x[k])-(x[k+2]-x[k])*(y[k]*x[k+1]-y[k+1]*x[k]))
        c_up = (X[k]*x[k+1]*(y[k]*x[k+2]-y[k+2]*x[k])-X[k+1]*x[k]*(y[k]*x[k+2]-y[k+2]*x[k])-X[k]*x[k+2]*(y[k]*x[k+1]-y[k+1]*x[k])-X[k+2]*x[k]*(y[k]*x[k+1]-y[k+1]*x[k]))
        c = c_up/c_down
        print('c',c)
        b = ((X[k]*x[k+1]-X[k+1]*x[k])-(x[k+1]-x[k])*c)/(y[k]*x[k+1]-y[k+1]*x[k])
        print('b',b)
        a = (X[k]-b*y[k]-c)/x[k]
        print('a', a)
        return(a,b,c)

    (a,b,c) = calc_plate_const(0,x_pix,y_pix,X_star)
    (d,e,f) = calc_plate_const(0,x_pix,y_pix,Y_star)
    print(target_RA,target_DEC)
    # 4) Calculate the standard coordinates for the object
    # HIER object at (0,0)
    (x_ob, y_ob) = (0,0)
    X_ob = a*x_ob + b*y_ob + c
    Y_ob = d*x_ob + e*y_ob + f
    print('x', x_pix, x_ob, 'y', y_pix, y_ob)
    print('X', X_star, X_ob, 'Y', Y_star, Y_ob)
    print('RA', ra_star, 'DEC', dec_star)

    # 5) Calculate the RA/DEC of the objects standard coordinates
    a = target_RA + np.arctan(DegToRad*((-X_ob)/(np.cos(DegToRad*target_DEC)- Y_ob*np.sin(DegToRad*target_DEC))))
    d = target_DEC -  np.arcsin(DegToRad*((np.sin(DegToRad*target_DEC) + Y_ob*np.cos(DegToRad*target_DEC))/(np.sqrt(1 + X_ob**2 + Y_ob**2))))
    print('RA in rad', a, 'DEC in rad', d)
    print('RA',a,target_RA, 'DEC',d, target_DEC)

    return(a,d)

输入例如是一个数组，其中星星的位置以图像的像素和真实世界的程度为单位

star pos = [[1948.2, 1205.8, 132.34058333333334, -3.4429722222222225], [153.90000000000001, 1892.5, 131.08620833333333, -5.0947499999999994]
# star_pos [x_pos in pix, y_pos in pix, RA, DEC]

(x_pix, y_pix) = (0.0135, 0.0135)
# pixel size

focallength = 0.7168

(target_RA, target_DEC) = (131.683014444 -3.91890194444)
# But I am not sure how accurate that is, it is more of an assumption. I would say if I look at the star map manually, it looks quite a bit off...

我希望看到 (0,0) 点 RA 约为 133° 和 DEC -5.75°

4

0 回答 0