fix PA bug

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])