In [292]:
%pylab
import scipy.interpolate
import desimodel.io
def get_radius_mm(theta):
"""
Returns the radius in mm given a radius in degrees using the platescale data
relative to the center of the focal plane as (0,0)
Parameters
----------
theta: A float that represents the angle from the center of the focal plane
"""
platescale = desimodel.io.load_platescale()
# Uses a quadratic one-dimensional interpolation to approximate the radius in degrees versus radius in mm
fn = scipy.interpolate.interp1d(platescale['theta'], platescale['radius'], kind = 'quadratic')
radius = float(fn(theta))
return radius
# Implements the haversine formula
def radec2xy(telra, teldec, ra, dec):
"""
Returns arrays of the x, y positions of given celestial objects
on the focal plane given an arbitrary telescope pointing in RA and Dec and
arrays of the RA and Dec of celestial objects in the sky
Parameters
----------
telra: a scalar float signifying the telescope's RA pointing in degrees
teldec: a scalar float signifying the telescope's Dec pointing in degrees
ra: An array of RA values for locations in the sky
dec: An array of declination values for locations in the sky
"""
import numpy as np
import math
# Inclination is 90 degrees minus the declination in degrees
inc = 90 - dec
x0 = np.sin(math.radians(inc)) * np.cos(math.radians(ra))
y0 = np.sin(math.radians(inc)) * np.sin(math.radians(ra))
z0 = np.cos(math.radians(inc))
coord = [x0, y0, z0]
print(coord)
# Clockwise rotation around y axis by declination of the tile center
decrotate = np.zeros(shape=(3,3))
teldec_rad = math.radians(teldec)
decrotate[0] = [np.cos(teldec_rad), 0, np.sin(teldec_rad)]
decrotate[1] = [0, 1, 0]
decrotate[2] = [-np.sin(teldec_rad), 0, np.cos(teldec_rad)]
# Clockwise rotation around the z-axis by the right ascension of the tile center
rarotate = np.zeros(shape=(3,3))
telra_rad = math.radians(telra)
rarotate[0] = [np.cos(telra_rad), np.sin(telra_rad), 0]
rarotate[1] = [-np.sin(telra_rad), np.cos(telra_rad), 0]
rarotate[2] = [0, 0, 1]
coord1 = np.matmul(rarotate, coord)
coord2 = np.matmul(decrotate, coord1)
x = coord2[0]
y = coord2[1]
z = coord2[2]
newteldec = 0
newtelra = 0
ra_rad = np.arctan2(y, x)
dec_rad = (np.pi / 2) - np.arccos(z / np.sqrt((x**2) + (y**2) + (z**2)))
radius_rad = 2 * np.arcsin(np.sqrt((np.sin((dec_rad - newteldec) / 2)**2) + ((np.cos(newteldec)) * np.cos(dec_rad) * (np.sin((ra_rad - newtelra) / 2)**2))))
radius_deg = math.degrees(radius_rad)
q_rad = np.arctan2(-z, -y)
q_deg = math.degrees(q_rad) % 360
if(q_deg == 360):
q_deg = 0
radius_mm = get_radius_mm(radius_deg)
# CAN DELETE Q_DEG: IT WAS USED FOR DEBUGGING AND ERROR CHECKING
x_focalplane = radius_mm * np.cos(q_rad)
y_focalplane = radius_mm * np.sin(q_rad)
return x_focalplane, y_focalplane
ra = [8.40634632111, 8.927313423598427]
dec = [-9.93649291992, -9.324956250231294]
x, y = radec2xy(8.37, -10.65, 8.927313423598427, -9.324956250231294)
#x, y = radec2xy(8.37, -10.65, ra, dec)
print(x)
print(y)
#8.40634632111 -9.93649291992
"""newra, newdec = xy2radec(8.37, -10.65, -138.345, -333.179)
[0.97483131854101557, 0.1531306530391553, -0.16203364925652533]
8.927313423598427
-9.324956250231294"""
Out[292]:
In [293]:
def is_on_gfa(tileid, ra, dec):
"""
Checks if a target is on any of the 10 GFAs given a tileid and an array of RA and Dec pointings.
Parameters
----------
tileid: An integer representing the tile's unique ID which can retrieve its RA and Dec pointing
ra: An array of RA values for locations in the sky
dec: An array of declination values for locations in the sky
"""
import desimodel.footprint
import desimodel.io
# If any calculated area is under the threshold area, it is mathematically impossible
THRESHOLD_AREA = 469.7
TOLERANCE = 100
# The area boundary's value is the area of the gfa plus some tolerance.
AREA_BOUNDARY = 469.8037221 + TOLERANCE
targetindices = []
gfatable = desimodel.io.load_gfa()
telra, teldec = desimodel.footprint.get_tile_radec(tileid)
#print(telra, teldec)
targetx, targety = desimodel.focalplane.radec2xy(telra, teldec, ra, dec)
targetx = np.asarray(targetx)
targety = np.asarray(targety)
# Method to check if point is inside the rectangle
for gfaid in range(0, 40, 4):
# a1 through a4 are edge lengths of the rectangle formed by corners of the GFAs
a1 = np.sqrt((gfatable['X'][gfaid] - gfatable['X'][gfaid + 1])**2
+ (gfatable['Y'][gfaid] - gfatable['Y'][gfaid + 1])**2)
a2 = np.sqrt((gfatable['X'][gfaid + 1] - gfatable['X'][gfaid + 2])**2
+ (gfatable['Y'][gfaid + 1] - gfatable['Y'][gfaid + 2])**2)
a3 = np.sqrt((gfatable['X'][gfaid + 2] - gfatable['X'][gfaid + 3])**2
+ (gfatable['Y'][gfaid + 2] - gfatable['Y'][gfaid + 3])**2)
a4 = np.sqrt((gfatable['X'][gfaid + 3] - gfatable['X'][gfaid])**2
+ (gfatable['Y'][gfaid + 3] - gfatable['Y'][gfaid])**2)
# b1 through b4 are the line segments from each corner to the target location
b1 = np.sqrt((gfatable['X'][gfaid] - targetx)**2
+ (gfatable['Y'][gfaid] - targety)**2)
b2 = np.sqrt((gfatable['X'][gfaid + 1] - targetx)**2
+ (gfatable['Y'][gfaid + 1] - targety)**2)
b3 = np.sqrt((gfatable['X'][gfaid + 2] - targetx)**2
+ (gfatable['Y'][gfaid + 2] - targety)**2)
b4 = np.sqrt((gfatable['X'][gfaid + 3] - targetx)**2
+ (gfatable['Y'][gfaid + 3] - targety)**2)
# Calculating areas of triangles using Heron's Formula
#gfaarea = a1 * a2
#gfaarea2 = a2 * a3
u1 = (a1 + b1 + b2) / 2.0
u2 = (a2 + b2 + b3) / 2.0
u3 = (a3 + b3 + b4) / 2.0
u4 = (a4 + b4 + b1) / 2.0
area1 = np.sqrt(u1 * (u1 - a1) * (u1 - b1) * (u1 - b2))
area2 = np.sqrt(u2 * (u2 - a2) * (u2 - b2) * (u2 - b3))
area3 = np.sqrt(u3 * (u3 - a3) * (u3 - b3) * (u3 - b4))
area4 = np.sqrt(u4 * (u4 - a4) * (u4 - b4) * (u4 - b1))
targetarea = area1 + area2 + area3 + area4
assert np.all(targetarea > THRESHOLD_AREA)
#use np.where!!!!!!!!!
# This if statement may be redundant and thus unnecessary
if(any(targetarea < AREA_BOUNDARY) and any(targetarea > THRESHOLD_AREA)):
print("The target is on GFA", int(gfaid / 4))
newtargetindices = np.where(targetarea < AREA_BOUNDARY)
targetindices.append(newtargetindices[0])
"""targetindices is an array of targets with their respective indices in the
RA and Dec list passed in that fall on certain GFAs denoted by the index
in the targetindices array."""
return targetindices
In [294]:
def is_on_gfa_test(targetx, targety, TOLERANCE):
import desimodel.footprint
import desimodel.io
# If any calculated area is under the threshold area, it is mathematically impossible
THRESHOLD_AREA = 469.7
#TOLERANCE = 400
# The area boundary's value is the area of the gfa plus some tolerance.
AREA_BOUNDARY = 469.8037221 + TOLERANCE
targetindices = []
gfatable = desimodel.io.load_gfa()
targetx = np.asarray(targetx)
targety = np.asarray(targety)
# Method to check if point is inside the rectangle
for gfaid in range(0, 40, 4):
# a1 through a4 are edge lengths of the rectangle formed by corners of the GFAs
a1 = np.sqrt((gfatable['X'][gfaid] - gfatable['X'][gfaid + 1])**2
+ (gfatable['Y'][gfaid] - gfatable['Y'][gfaid + 1])**2)
a2 = np.sqrt((gfatable['X'][gfaid + 1] - gfatable['X'][gfaid + 2])**2
+ (gfatable['Y'][gfaid + 1] - gfatable['Y'][gfaid + 2])**2)
a3 = np.sqrt((gfatable['X'][gfaid + 2] - gfatable['X'][gfaid + 3])**2
+ (gfatable['Y'][gfaid + 2] - gfatable['Y'][gfaid + 3])**2)
a4 = np.sqrt((gfatable['X'][gfaid + 3] - gfatable['X'][gfaid])**2
+ (gfatable['Y'][gfaid + 3] - gfatable['Y'][gfaid])**2)
# b1 through b4 are the line segments from each corner to the target location
b1 = np.sqrt((gfatable['X'][gfaid] - targetx)**2
+ (gfatable['Y'][gfaid] - targety)**2)
b2 = np.sqrt((gfatable['X'][gfaid + 1] - targetx)**2
+ (gfatable['Y'][gfaid + 1] - targety)**2)
b3 = np.sqrt((gfatable['X'][gfaid + 2] - targetx)**2
+ (gfatable['Y'][gfaid + 2] - targety)**2)
b4 = np.sqrt((gfatable['X'][gfaid + 3] - targetx)**2
+ (gfatable['Y'][gfaid + 3] - targety)**2)
# Calculating areas of triangles using Heron's Formula
#gfaarea = a1 * a2
#gfaarea2 = a2 * a3
u1 = (a1 + b1 + b2) / 2.0
u2 = (a2 + b2 + b3) / 2.0
u3 = (a3 + b3 + b4) / 2.0
u4 = (a4 + b4 + b1) / 2.0
area1 = np.sqrt(u1 * (u1 - a1) * (u1 - b1) * (u1 - b2))
area2 = np.sqrt(u2 * (u2 - a2) * (u2 - b2) * (u2 - b3))
area3 = np.sqrt(u3 * (u3 - a3) * (u3 - b3) * (u3 - b4))
area4 = np.sqrt(u4 * (u4 - a4) * (u4 - b4) * (u4 - b1))
targetarea = area1 + area2 + area3 + area4
assert np.all(targetarea > THRESHOLD_AREA)
#use np.where!!!!!!!!!
# This if statement may be redundant and thus unnecessary
if(any(targetarea < AREA_BOUNDARY) and any(targetarea > THRESHOLD_AREA)):
print("The target is on GFA", int(gfaid / 4))
newtargetindices = np.where(targetarea < AREA_BOUNDARY)
targetindices.append(newtargetindices[0])
"""targetindices is an array of targets with their respective indices in the
RA and Dec list passed in that fall on certain GFAs denoted by the index
in the targetindices array."""
#print(gfatable['X'], gfatable['Y'])
return targetindices
#targetindices = is_on_gfa_test(400, 20)
#print(targetindices)
In [295]:
import matplotlib.pyplot as plt
#telra and teldec = 7.11 -10.53
miny = -575
maxy = -150
minx = -100
maxx = 275
x = []
y = []
import desimodel.focalplane
import desimodel.footprint
while miny < maxy:
startx = minx
while startx < maxx:
x.append(startx)
y.append(miny)
startx += 5
miny += 5
#import desimodel.io
#tiles = desimodel.io.load_tiles()
x = np.asarray(x)
y = np.asarray(y)
targetindices = is_on_gfa_test(x, y, 2000)
plt.plot(x, y, 'b.')
for gfaid in targetindices:
for indices in gfaid:
plt.plot(x[indices], y[indices], 'r.')
targetindices = is_on_gfa_test(x, y, 0.001)
for gfaid in targetindices:
for indices in gfaid:
plt.plot(x[indices], y[indices], 'g.')
"""ra.append(8.452)
dec.append(-9.6977)
ra.append(6.644525362152656)
dec.append(-9.055)"""
#savefig('gfa2_tolerance_2000.jpg')
Out[295]:
In [7]:
def find_points_in_tel_range(telra, teldec, ra, dec, radius = None):
"""Return a list of indices of points that are within a radius of an arbitrary telra, teldec.
This function is optimized to query a lot of points with a single telra and teldec.
radius is in units of degrees. The return value is a list
that contains the index of points that are in each tile.
The indices are not sorted in any particular order.
if tiles is a scalar, a single list is returned.
default radius is from desimodel.focalplane.get_tile_radius_deg()
Note: This is simply a modified version of find_points_in_tiles, but this function does not know about tiles.
"""
from scipy.spatial import cKDTree as KDTree
import desimodel.focalplane
import numpy as np
if radius is None:
radius = desimodel.focalplane.get_tile_radius_deg()
# check for malformed input shapes. Sorry we currently only
# deal with vector inputs. (for a sensible definition of indices)
assert ra.ndim == 1
assert dec.ndim == 1
points = _embed_sphere(ra, dec)
tree = KDTree(points)
# radius to 3d distance
threshold = 2.0 * np.sin(np.radians(radius) * 0.5)
xyz = _embed_sphere(telra, teldec)
indices = tree.query_ball_point(xyz, threshold)
return indices
def _embed_sphere(ra, dec):
""" embed RA DEC to a uniform sphere in three-d """
import numpy as np
phi = np.radians(np.asarray(ra))
theta = np.radians(90.0 - np.asarray(dec))
r = np.sin(theta)
x = r * np.cos(phi)
y = r * np.sin(phi)
z = np.cos(theta)
return np.array((x, y, z)).T
import numpy as np
answer = find_points_in_tel_range(0, 0, np.array([1.5, 0, 1.9, 0]), np.array([0, 1.5, 0, 1.3]))
print(answer)