import astropy as ap
import numpy as np
from astropy.io import fits
import os
import math as m
import time as t
# ra = y, dec = x
# I mess around with the ordering of x and y coordinates in pairs, so beware!
def pixels_to_sky(x_coord, y_coord, head, convert=360000):
"""
Purpose: Convert pixel location in a .fits image to right ascension and declination in 0.01 arcseconds
:param x_coord: The x-coordinate to be converted
:param y_coord: The y-coordinate to be converted
:param head: The FITS header of the image. The centre point of the .fits image must be
defined; specifically the 'CRPIX1' and 'CRPIX2' in the FITS header (these entries should be in degrees)
:param convert: Defines a unit conversion away from degrees to the unit of each pixel in the .fits image
:return: The right ascension and declination of the given point in the image in 0.01 arcseconds as floats
"""
# right ascension is oriented on y-axis, declination is oriented on x-axis
y_height = head['CRPIX1']
x_width = head['CRPIX2']
ra_ref = head['CRVAL1'] * convert # 360000 puts them into 0.01 arcseconds, this was convenient for unit pixels
dec_ref = head['CRVAL2'] * convert
ra_coord = y_coord - y_height + ra_ref # distance between RA&Dec in sky is equal to pixel distance in Y&X
dec_coord = x_coord - x_width + dec_ref
return dec_coord, ra_coord
def sky_to_pixels(dec_coord, ra_coord, head, convert=360000):
"""
Purpose: Convert right ascension and declination to pixel location in a .fits image
:param dec_coord: The declination to be converted in 0.01 arcseconds
:param ra_coord: The right ascension to be converted in 0.01 arcseconds
:param head: The FITS header of the image. The centre point of the .fits image must be
defined; specifically the 'CRPIX1', 'CRPIX2','CRVAL1' and 'CRVAL2' in the FITS header (these entries should be in
degrees)
:param convert: Defines a unit conversion away from degrees to the unit of each pixel in the .fits image
:return: The x- and y-coordinate of the given point in the image
"""
y_height = head['CRPIX1']
x_width = head['CRPIX2']
ra_ref = head['CRVAL1'] * convert
dec_ref = head['CRVAL2'] * convert
x_coord = int(dec_coord - dec_ref + x_width)
y_coord = int(ra_coord - ra_ref + y_height)
return x_coord, y_coord
def write_point(x_coord, y_coord, data, mag=10):
"""
Purpose: Write a point onto a .fits file image array. The image must already be opened.
:param x_coord: x-coordinate in pixels of the point to be written
:param y_coord: y-coordinate in pixels of the point to be written
:param data: The data of the .fits image stored as an array
:param mag: Magnitude of the image as a float. Default is 10
:return: (none)
"""
# for our procedure we had normalized magnitudes, so we set pixel values arbitrarily (but similar)
data[y_coord, x_coord] = 2.512 ** mag
return
def get_info():
"""
Purpose: Prompt user input of the name of a file with each image RA, Dec (both in degrees), and magnitude
comma-separated, with each image on a new line, and return a list lists of these values for each image
:return: A list of lists, with each list corresponding to a single lensed image with the RA, Dec, and magnitude, as
well as the name of the input file
"""
chars = ["0", "1", "2", "3", "4", "5", "6", "7", "8", "9", ".", "-", "\n", ","]
image_info = input('Please input the name of a file with the appropriate information and format: ')
raw = open(image_info+'.csv', 'r') # pretty basic file reading stuff
boss = []
for line in raw:
for item in line:
if item not in chars: # check to ensure characters are only coordinates
print('Invalid character in .csv file found. Check for whitespaces or letters '
'(recommend opening in Notepad rather than a spreadsheet program). Program closing now.')
quit()
slave = line.rstrip().split(',')
boss.append(slave)
return boss, image_info
def write_all(length, convert=360000):
"""
Purpose: Makes a .fits image and inputs points representing quasar images, optionally with relative magnitudes
:param length: The length/width of the fits image to be created
:param convert: Defines a unit conversion away from degrees to the unit of each pixel in the .fits image
:return: The file name of a new spoofed .fits file, as well as a list-of-lists containing the pixel locations of
each image in (x,y) format
"""
books = [] # need this for 'bookkeeping' later.
info, quasar_name = get_info()
data = np.zeros((int(length), int(length))) # create array of appropriate size
head = fits.Header() # create the header
height, width = data.shape
y_height = height // 2 # find middle of image
x_width = width // 2
# if any other header parts need to be written, this is the place to do it
head['CRPIX1'] = y_height # arbitrarily puts the first image at the centre
head['CRPIX2'] = x_width
image_prime = info[0]
# ra_prime = int(float(image_prime[0]) * convert)
# dec_prime = int(float(image_prime[1]) * convert)
books.append([y_height, x_width])
head['CRVAL1'] = float(image_prime[0]) # in degrees
head['CRVAL2'] = float(image_prime[1])
head['CUNIT1'] = '0.01 arcseconds'
head['CUNIT2'] = '0.01 arcseconds'
write_point(x_width, y_height, data) # prime image is written and coordinate system defined
for image in info[1:]: # no need to take first image here; it is already written!
ra_image = int(float(image[0]) * convert)
dec_image = int(float(image[1]) * convert)
# mag = float(image[2]) # not using this right now, just normalizing image
x_image, y_image = sky_to_pixels(dec_image, ra_image, head)
write_point(x_image, y_image, data)
books.append([y_image, x_image])
fits.writeto("spoofed_" + quasar_name + '.fits', data, head)
print("A new spoofed file should be available in the root folder.")
return "spoofed_" + quasar_name + '.fits', books
def semi_minor(a, x, y, x_c, y_c, angle):
"""
Calculates the sum of the errors of putting the points to be fitted into the ellipse equation with given
coefficients
:param a: The semi-major axis of the ellipse
:param x: The x-coordinate of the point being calculated
:param y: The y-coordinate of the point being calculated
:param x_c: The x-coordinate of the center point
:param y_c: The y-coordinate of the center point
:param angle: The angle between the positive x-axis and the center point, rotated counter-clockwise around the tip
:return: The semi-minor axis as calculated at this point, None if a mathematical singularity is found
"""
b = None
numerator = (- (a ** 2) * ((y_c * m.cos(angle) - y * m.cos(angle) - x_c * m.sin(angle) + x * m.sin(angle)) ** 2))
denominator = (-(a ** 2)
+ (x_c ** 2) * (m.cos(angle) ** 2) - 2 * x_c * x * (m.cos(angle) ** 2)
+ (x ** 2) * (m.cos(angle) ** 2)
+ 2 * x_c * y_c * m.cos(angle) * m.sin(angle)
- 2 * x * y_c * m.cos(angle) * m.sin(angle)
- 2 * x_c * y * m.cos(angle) * m.sin(angle)
+ 2 * x * y * m.cos(angle) * m.sin(angle)
+ (y_c ** 2) * (m.sin(angle) ** 2)
- 2 * y_c * y * (m.sin(angle) ** 2)
+ (y ** 2) * (m.sin(angle) ** 2))
if a != 0 and denominator != 0: # check for mathematical singularities
argument = numerator / denominator
if argument > 0:
b = m.sqrt(argument) # from equation for an ellipse, derived using Mathematica
return b
def ellipse_insanity(books, error=10):
"""
Fits an ellipse to a lensed quasar image via a brute-force algorithm
:param books: A list-of-lists containing the y,x coordinates of each quasar image
:param error: The error range in pixels while comparing semi-minor axis calculations for each point
:return: The a list of lists of semi-major, semi-minor, center-x, center-y and theta values for fitting ellipses
"""
# figure out how large the major axis should be, and decide which of the two extremum points is the tip
x_t = int(input("Please input the FITS X of the first opposite image: "))
y_t = int(input("Please input the FITS Y of the first opposite image: "))
x_tw = int(input("Please input the FITS X of the second opposite image: "))
y_tw = int(input("Please input the FITS Y of the second opposite image: "))
opposites = m.sqrt((x_t - x_tw) ** 2 + (y_t - y_tw) ** 2)
start_time = t.time() # just for fun
a_distance = 0
b_distance = 0
for item in books:
a_distance += m.sqrt((x_t - item[1]) ** 2 + (y_t - item[0]) ** 2)
b_distance += m.sqrt((x_tw - item[1]) ** 2 + (y_tw - item[0]) ** 2)
if a_distance > b_distance:
tip = y_t, x_t
else:
tip = y_tw, x_tw
# determine max and min values of theta
omegas = []
for item in books:
shorta = item[1]-tip[1] # dec difference
shortb = item[0]-tip[0] # ra difference
if shorta != 0 and shortb != 0:
longc = m.sqrt(shorta ** 2 + shortb ** 2) # point distance
bigsee = m.acos((longc ** 2 - shorta ** 2 - shortb ** 2)/(-2 * shorta * shortb)) # law of cosines
omegas.append(bigsee)
min_omega = min(omegas)
max_omega = max(omegas)
# bounds for searching for the centre (it must be inside the image shape)
y_vals = [item[0] for item in books]
x_vals = [item[1] for item in books]
x_max = max(x_vals)
x_min = min(x_vals)
y_max = max(y_vals)
y_min = min(y_vals)
positives = []
calculations = 0
# brute-force the centre point of the ellipse, increasing allowed error in minor axis calc by 5 each time
# zero results are found
while len(positives) == 0 and error < 100:
for y_c in range(y_min, y_max + 1):
for x_c in range(x_min, x_max + 1):
# calculate theta with law of cosines
part1 = x_c - tip[1] # dec difference
part2 = y_c - tip[0] # ra difference
part3 = m.sqrt(part1 ** 2 + part2 ** 2)
if part1 != 0 and part2 != 0:
theta = m.acos((part3 ** 2 - part1 ** 2 - part2 ** 2)/(-2 * part1 * part2))
elif part1 == 0 and tip[0] < y_c:
theta = m.pi/2 # special case where theta is 90 degrees (tip is below center)
elif part1 == 0 and tip[0] > y_c:
theta = m.pi / 2 # special case where theta is 270 degrees (tip is above center)
elif part2 == 0 and tip[1] > x_c:
theta = m.pi # special case where theta is 180 degrees (tip is right of center)
else:
theta = 0 # special case where theta is 0 degrees (tip is left of center)
a = m.sqrt((tip[1] - x_c) ** 2 + (tip[0] - y_c) ** 2)
b_main = semi_minor(a, tip[1], tip[0], x_c, y_c, theta)
mini_flag = []
for item in books: # calculates the theoretical minor axis for all points and compares them
b_item = semi_minor(a, item[1], item[0], x_c, y_c, theta)
if b_item is None or b_main is None:
mini_flag.append(False)
else: # checks if b is within error, both b's are not None, theta is within bounds, and a is
# sufficiently large to make physical sense (larger than half the opposite points distance)
mini_flag.append(abs(b_main-b_item) <= error and a >= 0.5 * opposites
and max_omega >= theta >= min_omega)
if False not in mini_flag:
positives.append([a, b_main, x_c, y_c, theta, a - 0.5 * opposites])
calculations += 1
error += 5
# just some statistics for testing
print("\nFound " + str(len(positives)) + " ellipses in " + str((t.time() - start_time)) +
" seconds, with semi-minor axis error of " + str(error - 5))
if error >= 100:
print("Minor-axis error has gotten too large. System is unable to be fitted by this algorithm. Program closing")
quit()
# pick average result for each parameter
ave_result = []
for i in range(0, 5): # average all values
mini_result = []
for j in positives:
mini_result.append(j[i])
ave_result.append(sum(mini_result) / len(mini_result)) # probably the reason results are mainly circular
return ave_result
def lens_inversion(image, x_c, y_c, convert=360000):
"""
Purpose: Finds the coordinates of the lens in degrees by reflecting the luminosity centroid location through the
center point of the fitted ellipse
:param image: The image of the lensed quasar
:param x_c: The x-coordinate of the centre of the drawn ellipse
:param y_c: The y-coordinate of the centre of the drawn ellipse
:param convert: Defines a unit conversion away from the units in the .fits file, which should have unitary pixels
:return: The right ascension and declination of the lens in degrees, as well as the pixel location of the lens
"""
x_centroid = int(input("Please input the FITS X of the centroid: "))
y_centroid = int(input("Please input the FITS Y of the centroid: "))
x_lens = 2 * x_c - x_centroid # flip centroid over major and minor axes to get lens location
y_lens = 2 * y_c - y_centroid
hdul = fits.open(image) # convert to astronomy coordinates
head = hdul[0].header
dec_lens, ra_lens = pixels_to_sky(x_lens, y_lens, head)
hdul.close()
return ra_lens / convert, dec_lens / convert, x_lens, y_lens
def ellipse_draw(ellipse, spoofed, books, length=1200, error=0.01):
"""
Trace an ellipse onto an array
:param ellipse: A list of semi-major axis, semi-minor axis, centre-x, centre-y, and tilt angle of an ellipse
:param spoofed: The filename of the spoofed ellipse
:param books: The points onto which the ellipse fits
:param length: The size of the image
:param error: The error range in which to draw
:return: (none)
"""
new_data = np.zeros((int(length), int(length)))
a = ellipse[0]
b = ellipse[1]
x_c = ellipse[2]
y_c = ellipse[3]
theta = ellipse[4]
phase = 0
# Find absolute bounds for ellipse
y_vals = [item[0] for item in books]
x_vals = [item[1] for item in books]
x_max = max(x_vals)
x_min = min(x_vals)
y_max = max(y_vals)
y_min = min(y_vals)
for y in range(y_min-80, y_max+81): # need large bounds to get a really visible image
for x in range(x_min-80, x_max+81):
if abs(((((x - x_c) * m.cos(theta+phase)+(y-y_c) * m.sin(theta+phase)) ** 2) / a ** 2
+ (((x - x_c) * m.sin(theta+phase)-(y-y_c) * m.cos(theta+phase)) ** 2) / b ** 2)-1) <= error:
new_data[y, x] = 100
fits.writeto(str(spoofed) + "_drawn.fits", new_data)
return
def master_function(draw=False):
"""
Purpose: Encompass all functions for lens inversion and allow user to perform the operation for multiple lensed
quasars or files.
:param draw: Toggles the drawing function
:return: (none)
"""
# 1200 chosen to represent 0.2' (600 pixels = 0.1', 1 pixel = 0.01")
terminate = 'n'
while terminate != 'y':
spoofed, images = write_all(1200)
print("Please open the new spoofed file in AstroImageJ.")
ellipse = ellipse_insanity(images)
x_c = ellipse[2]
y_c = ellipse[3]
ra_lens, dec_lens, x_lens, y_lens = lens_inversion(spoofed, x_c, y_c)
print("\nLens information: ")
print("Right Ascension in Degrees: " + str(ra_lens))
print("Declination in Degrees: " + str(dec_lens))
print("X-coordinate on image: " + str(int(x_lens)))
print("Y-coordinate on image: " + str(int(y_lens)))
if draw:
ellipse_draw(ellipse, spoofed, images)
terminate = input("If you wish to exit this session, type 'y'. Other wise, type 'n' or another letter to solve "
"a new lensed system. ")
return
master_function(draw=True)
print("Program ran to completion.")