From a626ceb8bd731649d2c317964fd945ea2dea7490 Mon Sep 17 00:00:00 2001 From: xin Date: Tue, 13 Sep 2022 13:21:09 +0800 Subject: [PATCH] fix PA bug --- getPointingList.py | 36 +++++++++++++++++++++++++++++++++--- 1 file changed, 33 insertions(+), 3 deletions(-) diff --git a/getPointingList.py b/getPointingList.py index 979e793..ed27eee 100755 --- a/getPointingList.py +++ b/getPointingList.py @@ -4,6 +4,33 @@ from pylab import * import numpy as np import math as m # test + +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 getobsPA(ra, dec): + l1 = np.array([0,0,1]) + l2 = transRaDec2D(ra, dec) + polar_ec = coord.SkyCoord(0*u.degree, 90*u.degree,frame='barycentrictrueecliptic') + polar_eq = polar_ec.transform_to('icrs') + + # print(polar_eq.ra.value,polar_eq.dec.value) + polar_d = transRaDec2D(polar_eq.ra.value, polar_eq.dec.value) + l1l2cross = np.cross(l2,l1) + pdl2cross = np.cross(l2,polar_d) + 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 + return angle + + + def loadSatOrbitDat(datDir='',fileNum=50): oData = loadtxt(datDir+'1.txt') @@ -81,10 +108,11 @@ def producePointingList(out = '' , center = [60,-40], radius = 5): for p,p_ in zip(points[ids1][ids2][ids3][ids4],c_equtor[ids1][ids2][ids3][ids4]): ra = (p_.ra*u.degree).value dec = (p_.dec*u.degree).value + pa = getobsPA(ra, dec) + 90 # print(ra, dec) lon = p[2] lat = p[1] - elip_txt = elip_txt + str(round(ra,6)) + ' ' + str(round(dec,6)) + ' ' + str(round(lon,6)) + ' ' + str(round(lat,6)) + ' -113.4333' + elip_txt = elip_txt + str(round(ra,6)) + ' ' + str(round(dec,6)) + ' ' + str(round(lon,6)) + ' ' + str(round(lat,6)) + ' ' + str(pa) elip_txt = elip_txt + '\n' # for i in np.arange(0,num,1): @@ -124,10 +152,11 @@ def producePointingList2(out = '' , center = [60,-40], radius = 5, survey_file = for p,p_ in zip(points[ids1][ids2][ids3][ids4],c_equtor[ids1][ids2][ids3][ids4]): ra = (p_.ra*u.degree).value dec = (p_.dec*u.degree).value + pa = getobsPA(ra, dec) + 90 # print(ra, dec) lon = p[2] lat = p[1] - elip_txt = elip_txt + str(round(ra,6)) + ' ' + str(round(dec,6)) + ' ' + str(round(lon,6)) + ' ' + str(round(lat,6)) + ' -113.4333' + elip_txt = elip_txt + str(round(ra,6)) + ' ' + str(round(dec,6)) + ' ' + str(round(lon,6)) + ' ' + str(round(lat,6)) + ' ' + str(pa) oTime = p[0] satPos, sOrbitId =locateSat(time=oTime, OrbitData = orbitDat, startId = sOrbitId, orbDataLen = orbitDat.shape[0]) @@ -206,9 +235,10 @@ def producePointingList3(out = '' , center = [60,-40], radius = 5, survey_file = ra = (p_.ra*u.degree).value dec = (p_.dec*u.degree).value # print(ra, dec) + pa = getobsPA(ra, dec) + 90 lon = p[2] lat = p[1] - elip_txt = elip_txt + str(round(ra,6)) + ' ' + str(round(dec,6)) + ' ' + str(round(lon,6)) + ' ' + str(round(lat,6)) + ' -113.4333' + elip_txt = elip_txt + str(round(ra,6)) + ' ' + str(round(dec,6)) + ' ' + str(round(lon,6)) + ' ' + str(round(lat,6)) + ' ' + str(pa) oTime = p[0] satPos, sOrbitId =locateSat(time=oTime, OrbitData = orbitDat, startId = sOrbitId, orbDataLen = orbitDat.shape[0]) -- GitLab