test

getInputPointing.py

@@ -4,6 +4,7 @@ from astropy import units as u
 from pylab import *
 import numpy as np
 import galsim
+import math
 # from numba import jit
 
 class Chip(object):
@@ -111,20 +112,49 @@ def getSelectPointingList(center = [60,-40], radius = 2):
     
     center = center#ra dec
     radius = radius # degree
+
+    radii_dec = 1
+    radii_ra = 1/math.cos(math.pi*center[1]/180)
+
+    if radii_ra > 180:
+        radii_ra = 180
+
     c_eclip = coord.SkyCoord(points[:,2]*u.degree, points[:,1]*u.degree,frame='barycentrictrueecliptic')
     c_equtor = c_eclip.transform_to('icrs')
 
-    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
+    # print(np.min((c_equtor.ra*u.degree).value), np.max((c_equtor.ra*u.degree).value))
+
+    # c_equtor_sel = c_equtor
+    # points_sel = points
+
+    ra_range_lo = center[0]-radii_ra
+    ra_range_hi = center[0]+radii_ra
+
+    if ra_range_lo< 0:
+        ids1 = ((c_equtor.ra*u.degree).value<ra_range_hi) | ((c_equtor.ra*u.degree).value>360+ra_range_lo)
+    elif ra_range_hi > 360:
+        ids1 = ((c_equtor.ra*u.degree).value>ra_range_lo) | ((c_equtor.ra*u.degree).value<ra_range_hi-360)
+    else:
+        ids1 = ((c_equtor.ra*u.degree).value > ra_range_lo) & ((c_equtor.ra*u.degree).value < ra_range_hi)
+
+
+    dec_range_lo = center[1]-radii_dec
+    if center[1]-radii_dec < -90:
+        dec_range_lo = -90
+
+    dec_range_hi = center[1]+radii_dec
+    if center[1]+radii_dec > 90:
+        dec_range_hi = 90
+    
+    ids3 = (c_equtor[ids1].dec*u.degree).value > dec_range_lo
+    ids4 = (c_equtor[ids1][ids3].dec*u.degree).value < dec_range_hi
 
-    num = points[ids1][ids2][ids3][ids4].shape[0]
+    num = points[ids1][ids3][ids4].shape[0]
 
     p_result = np.zeros([num, 5])
     i = 0
 
-    for p,p_ in zip(points[ids1][ids2][ids3][ids4],c_equtor[ids1][ids2][ids3][ids4]):
+    for p,p_ in zip(points[ids1][ids3][ids4],c_equtor[ids1][ids3][ids4]):
         ra = (p_.ra*u.degree).value
         dec = (p_.dec*u.degree).value
         # print(ra, dec)