Commit 4e8cf2dc authored by Yan Zhaojun's avatar Yan Zhaojun
Browse files

debug

parent 5bb06e78
Loading
Loading
Loading
Loading
+19 −19
Original line number Diff line number Diff line
@@ -24,31 +24,31 @@ def transRaDec2D(ra, dec):
    return np.array([x1, y1, z1])


def earth_angle(time_jd, x_sat, y_sat, z_sat, ra_obj, dec_obj):
# def earth_angle(time_jd, x_sat, y_sat, z_sat, ra_obj, dec_obj):

    ra_sat = np.arctan2(y_sat, x_sat) / np.pi * 180
    dec_sat = np.arctan2(z_sat, np.sqrt(x_sat**2+y_sat**2)) / np.pi * 180
    radec_sat = SkyCoord(ra=ra_sat*u.degree,
                         dec=dec_sat*u.degree, frame='gcrs')
    lb_sat = radec_sat.transform_to('geocentrictrueecliptic')
#     ra_sat = np.arctan2(y_sat, x_sat) / np.pi * 180
#     dec_sat = np.arctan2(z_sat, np.sqrt(x_sat**2+y_sat**2)) / np.pi * 180
#     radec_sat = SkyCoord(ra=ra_sat*u.degree,
#                          dec=dec_sat*u.degree, frame='gcrs')
#     lb_sat = radec_sat.transform_to('geocentrictrueecliptic')

    # get the obj location
    radec_obj = SkyCoord(ra=ra_obj*u.degree,
                         dec=dec_obj*u.degree, frame='gcrs')
    lb_obj = radec_obj.transform_to('geocentrictrueecliptic')
#     # get the obj location
#     radec_obj = SkyCoord(ra=ra_obj*u.degree,
#                          dec=dec_obj*u.degree, frame='gcrs')
#     lb_obj = radec_obj.transform_to('geocentrictrueecliptic')

    # calculate the angle between sub-satellite point and the earth side
    earth_radius = 6371     # km
    sat_height = np.sqrt(x_sat**2 + y_sat**2 + z_sat**2)
    angle_a = np.arcsin(earth_radius/sat_height) / np.pi * 180
#     # calculate the angle between sub-satellite point and the earth side
#     earth_radius = 6371     # km
#     sat_height = np.sqrt(x_sat**2 + y_sat**2 + z_sat**2)
#     angle_a = np.arcsin(earth_radius/sat_height) / np.pi * 180

    # calculate the angle between satellite position and the target position
    angle_b = lb_sat.separation(lb_obj)
#     # calculate the angle between satellite position and the target position
#     angle_b = lb_sat.separation(lb_obj)

    # calculat the earth angle
    angle = 180 - angle_a - angle_b.degree
#     # calculat the earth angle
#     angle = 180 - angle_a - angle_b.degree

    return angle
#     return angle

###############################################################################