From 8709ad0dec8dd575173d46f35712ada9deb1425b Mon Sep 17 00:00:00 2001
From: xin <zhangxinbjfu@gmail.com>
Date: Mon, 18 Apr 2022 09:52:52 +0800
Subject: [PATCH] circle range

---
 getPointingList.py | 85 +++++++++++++++++++++++++++++++++++++++++++---
 1 file changed, 80 insertions(+), 5 deletions(-)

diff --git a/getPointingList.py b/getPointingList.py
index 49cdd4f..979e793 100755
--- a/getPointingList.py
+++ b/getPointingList.py
@@ -2,6 +2,7 @@ import astropy.coordinates as coord
 from astropy import units as u
 from pylab import *
 import numpy as np
+import math as m
 # test
 def loadSatOrbitDat(datDir='',fileNum=50):
     oData = loadtxt(datDir+'1.txt')
@@ -12,6 +13,11 @@ def loadSatOrbitDat(datDir='',fileNum=50):
     
     return oData
 
+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 locateSat(time=2459766., OrbitData = None, startId = 0, orbDataLen = 10000):
@@ -169,15 +175,84 @@ def producePointingList2(out = '' , center = [60,-40], radius = 5, survey_file =
     pointfn.close()
 
 
+def producePointingList3(out = '' , center = [60,-40], radius = 5, survey_file = 'skyMapOrSurveyList/E17.5_b17.5_beta_11.6_opt_transtime_1_CMG_1_dp_2_0.25_da_10_Texp_1.5_DEC60_500_0.1_800_1000_+5deg.dat'):
+    points = np.loadtxt(survey_file)
+
+    num = points.shape[0]
+    center = center#ra dec
+    radius = radius # degree
+    c_eclip = coord.SkyCoord(points[:,2]*u.degree, points[:,1]*u.degree,frame='barycentrictrueecliptic')
+    c_equtor = c_eclip.transform_to('icrs')
+
+    cent_xyz = transRaDec2D(center[0],center[1])
+    xyz = transRaDec2D((c_equtor.ra*u.degree).value,(c_equtor.dec*u.degree).value)
+
+    xyz_dp = xyz[0]*cent_xyz[0] + xyz[1]*cent_xyz[1] + xyz[2]*cent_xyz[2]
+
+    ids = xyz_dp>=m.cos(radius*m.pi/180.)
+
+    # ids1 = (c_equtor.ra*u.degree).value > center[0]-radius
+    # ids2 = (c_equtor[ids1].ra*u.degree).value < center[0]+radius
+    # ids3 = (c_equtor[ids1][ids2].dec*u.degree).value > center[1]-radius
+    # ids4 = (c_equtor[ids1][ids2][ids3].dec*u.degree).value < center[1]+radius
+
+    obDataDir = 'orbit20160925/'
+    orbitDat = loadSatOrbitDat(obDataDir,50)
+    sOrbitId = 0
+
+    elip_txt = '# ra    dec    lon(ecliptic)    lat(ecliptic)    pos_angle    time(julian)    sat_x    sat_y    sat_z    sun_x    sun_y    sun_z    moon_x    moon_y    moon_z    sat_vx    sat_vy    sat_vz    exp_time    isDeep\n'
+
+    for p,p_ in zip(points[ids],c_equtor[ids]):
+        ra = (p_.ra*u.degree).value
+        dec = (p_.dec*u.degree).value
+        # 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'
+
+        oTime = p[0]
+        satPos, sOrbitId =locateSat(time=oTime, OrbitData = orbitDat, startId = sOrbitId, orbDataLen = orbitDat.shape[0])
+        tempOrbitId = sOrbitId
+        deltT = 0.1 # unit s
+        nTime = p[0]+deltT/86400.
+        satPosN, _ =locateSat(time=nTime, OrbitData = orbitDat, startId = sOrbitId, orbDataLen = orbitDat.shape[0])
+        sat_v = (satPosN - satPos)/deltT
+        mm1 = np.sqrt(np.sum(sat_v*sat_v))
+        mm2 = np.sqrt(np.sum(satPos*satPos))
+        # print(np.sum(sat_v*satPos)/(mm1*mm2))
+        # print(sat_v,p[6:9])
+        if(np.abs(np.sum(sat_v*satPos)/(mm1*mm2))>0.01):
+            print(oTime)
+            break
+
+        
+        dat_col = [0,3,4,5,6,7,8,9,10,11,16,14]
+        for col in dat_col[:-2]:
+            elip_txt = elip_txt + '    ' + str(p[col])
+        elip_txt = elip_txt + '    ' + str(sat_v[0]) + '    ' + str(sat_v[1]) + '    ' + str(sat_v[2])
+
+        for col in dat_col[-2:]:
+            elip_txt = elip_txt + '    ' + str(p[col])
+
+
+        elip_txt = elip_txt + '\n'
+        
+
+    pointfn = open(out,'w')
+    pointfn.write(elip_txt)
+    pointfn.flush()
+    pointfn.close()
+
+
 if __name__ == "__main__":
-    isRealSurvey = False
+    isRealSurvey = True
     survey_file = 'skyMapOrSurveyList/E17.5_b17.5_beta_11.6_opt_transtime_1_CMG_1_dp_2_0.25_da_10_Texp_1.5_DEC60_500_0.1_800_1000_+5deg.dat'
-    outFileName = 'pointing_test_false.dat'
-    radi = 1.3
+    outFileName = 'pointing_test_NGP_1..dat'
+    radi = 1.2
 
-    center_pos = [60, -40]
+    center_pos = [192.8595, 27.1283]
     if isRealSurvey:
-        producePointingList2(out = outFileName,center = center_pos, radius = radi, survey_file = survey_file)
+        producePointingList3(out = outFileName,center = center_pos, radius = radi, survey_file = survey_file)
     else:
         producePointingList(out = outFileName,center = center_pos, radius = radi)
 
-- 
GitLab