In [245]:
# Working on identifying which celestial objects land on the GFAs
from astropy.io import fits
from astropy.table import Table
import desimodel.io
import desimodel.footprint
%pylab
sweep = Table.read('sweep-050p030-060p035.fits')
sweep1 = Table.read('sweep-060p065-070p070.fits')
plot(sweep['RA'], sweep['DEC'], 'b.')
tiles = desimodel.io.load_tiles()
#print(tiles.dtype.names)
#print(tiles['TILEID'])
radec = desimodel.footprint.get_tile_radec(4)
#print(radec)
#print(tiles['RA'])
#print(tiles['DEC'])
for targets in sweep:
if(targets['FLUX_R'] > 0):
print(targets['FLUX_R'])
#sweep
sweep1
Out[245]:
In [237]:
"""newra, newdec = xy2radec(8.37, -10.65, -138.345, -333.179)
8.927313423598427
-9.324956250231294"""
import desimodel.io
gfatable = desimodel.io.load_gfa()
gfatable
import desimodel.focalplane
tiles = desimodel.io.load_tiles()
#print(indices)
#indices = desimodel.footprint.find_tiles_over_point(tiles, 8.37, -10.65)
#[8254, 15072, 198, 2213, 4225, 6242, 8255, 10265, 12279]
#print(tiles[8254])
#tile index 8254 has tileid 23658
telra, teldec = desimodel.footprint.get_tile_radec(23658)
#print(telra, teldec)
ra, dec = desimodel.focalplane.xy2radec(7.11, -10.53, 116.279135121, -372.885546514)
#print(ra, dec)
#Certain gfa location is in ra, dec respectively 6.644525362152656 -9.055425745149217
In [238]:
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
MIN_TOLERANCE = .001
# The area boundary's value is the area of the gfa plus some tolerance.
AREA_BOUNDARY = 469.8037221 + MIN_TOLERANCE + 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
#is_on_gfa(23658, 6.644525362152656, -9.055425745149217)
import desimodel.focalplane
#test point: (x,y) = (119, -380)
#test tileid = 23658
telra1, teldec1 = desimodel.footprint.get_tile_radec(23658)
testra, testdec = desimodel.focalplane.xy2radec(telra1, teldec1, -335, -210)
#is_on_gfa(23658, testra, testdec)
#print(testra, testdec)
testra2 = [testra, 6.644525362152656, 7]
testdec2 = [testdec, -9.055425745149217, -10]
is_on_gfa(23658, testra2, testdec2)
gfatable
Out[238]:
In [239]:
import matplotlib.pyplot as plt
#telra and teldec = 7.11 -10.53
"""mindec = -9.8
maxdec = -9.6
minra = 6.5
maxra = 8.47"""
mindec = -12.5
maxdec = -7
minra = 5.4
maxra = 8.9
ra = []
dec = []
import desimodel.focalplane
import desimodel.footprint
while mindec < maxdec:
startra = minra
while startra < maxra:
ra.append(startra)
dec.append(mindec)
startra += .05
mindec += .05
import desimodel.io
tiles = desimodel.io.load_tiles()
ra = np.asarray(ra)
dec = np.asarray(dec)
"""index = 0
for specific_tile in tiles:
if(tiles[index]['TILEID'] == 23658):
print(index)
index += 1
#INDEX IS TILES[8254]"""
points = desimodel.footprint.find_points_in_tiles(tiles[8254], ra, dec)
"""for lists in points:
if lists:
print(lists)"""
#print(points)
new_ra = []
new_dec = []
for point_indices in points:
plt.plot(ra[point_indices], dec[point_indices], 'b.')
new_ra.append(ra[point_indices])
new_dec.append(dec[point_indices])
targetindices = is_on_gfa(23658, new_ra, new_dec)
for gfaid in targetindices:
for indices in gfaid:
plt.plot(new_ra[indices], new_dec[indices], 'g.')
"""ra.append(8.452)
dec.append(-9.6977)
ra.append(6.644525362152656)
dec.append(-9.055)"""
"""targetindices = is_on_gfa(23658, ra, dec)
print(targetindices)
#print(tiles.dtype.names)
plt.plot(ra, dec, 'b.')
for gfaid in targetindices:
for indices in gfaid:
plt.plot(ra[indices], dec[indices], 'g.')
#ra.remove(ra[indices])
#dec.remove(dec[indices])"""
#savefig('entire_gfa_tolerance_100.jpg')
Out[239]: