Galactic Redshift Script

On my spectral analysis of a quasar project, I had cheated by already knowing the redshift of the quasar from professional research, and was able to reverse-discover what elements should be in my spectrum. It led me to wonder how I might discover the elements in my spectra naturally. The ratios of the observed lines is really all there is to go off of, but there are just too many combinations of possible elements observable in a galaxy to guess this by hand. So, I wrote a Python script to do this for me.

I used a table of emission lines observed in galaxies, converted to a spreadsheet. The script takes every combination of elements in this chart, computes the ratio of the wavelengths of the longer line to the shorter line, and compares it to the observed wavelengths entered in the command line. For example, using the data I had from my quasar project, I would use the script thusly,

% python3.8 redshift.py 445.333 654.016 730.710

You simply just run the script with every emission line mean wavelength from the command line. It will first show you the ratios of all the combinations of your wavelengths,

Observed wavelengths and ratios
1        4453.33000
2        6540.16000
   6540.16000 / 4453.33000 = 1.468600
3        7307.10000
   7307.10000 / 4453.33000 = 1.640817
   7307.10000 / 6540.16000 = 1.117266

After it chugs through a large amount of possibilities from the reference data, it will show you the top ten matches based on the RMS of the errors between the observed and rest ratios, e.g., here are the top three from my data,

RMS: 0.000230
Rest Ratios:
[[1.4684041497472178], [1.6410161306706743, 1.1175507308074357]]
Elements:
[Fe V]: 3839.270000 -> 4453.330000
[Fe VI]: 5637.600000 -> 6540.160000
[O I]: 6300.304000 -> 7307.100000
Redshift:
0.1599465122528052

RMS: 0.000426
Rest Ratios:
[[1.4683580844685535], [1.6412941771876017, 1.1177751493646317]]
Elements:
C III]: 1908.734000 -> 4453.330000
Mg II]: 2802.705000 -> 6540.160000
O III: 3132.794000 -> 7307.100000
Redshift:
1.3330348690831522

RMS: 0.000835
Rest Ratios:
[[1.467395237750597], [1.6406031906864216, 1.1180376959661793]]
Elements:
[Ne III]: 3868.760000 -> 4453.330000
[Fe VI]: 5677.000000 -> 6540.160000
Si II: 6347.100000 -> 7307.100000
Redshift:
0.15146512185773703

The Gaussian fits of my data were not great looking back at what I did, a noobie mistake, but the real emission lines come up as the second best fit! So if you emission lines are measured accurately enough, this script should give you the elements in your spectrum as well as the calculated mean redshift.

The script and data file can be downloaded here. You will need the pandas library installed for the LibreOffice file operations. The full script is shown below,

# -*-coding: utf-8 -*-
import sys
import math
import statistics

from itertools import combinations
import pandas as pd

def rms(array):
    n = len(array)
    square = 0.0
    root = 0.0
    mean = 0.0
    
    #calculating Squre
    for i in range(0, n):
        square += (array[i] ** 2)
    #Calculating Mean
    mean = (square/ (float)(n))
    #Calculating Root
    root = math.sqrt(mean)
    return root

def flatten_list(_2d_list):
    flat_list = []
    # Iterate through the outer list
    for element in _2d_list:
        if type(element) is list:
            # If the element is of type list, iterate through the sublist
            for item in element:
                flat_list.append(item)
        else:
            flat_list.append(element)
    return flat_list


os = [float(a) for a in sys.argv[1:]]

os.sort()

ratios = []

if os[0] < 1000:
   os = [10 * o for o in os]

print("Observed wavelengths and ratios")

for i, w in enumerate(os):   
   print("%i        %0.5f" % (i + 1, w))
   
   rs = []
   
   for j in range(0, i):
      rs.append(os[i] / os[j])
      print("   %0.5f / %0.5f = %0.6f" % (os[i], os[j], os[i] / os[j]))
   
   if len(rs) > 0:
      ratios.append(rs)

data = pd.read_excel('data.ods', engine='odf')

print("")
print("")

errors = []

max_row = 0

for i in range(0, len(data.index)):
   if os[-1] >= data.at[i, "λ (Å)"]:
      max_row = i

combs = list(combinations(range(0, max_row), len(os)))

for comb in combs:
   ratios_e = []
   
   for i in range(0, len(comb)):
      r_es = []
      
      for j in range(0, i):
         r_es.append(data.at[comb[i], "λ (Å)"] / data.at[comb[j], "λ (Å)"])
         #print("   %0.5f / %0.5f = %0.6f" % (data.at[i, "λ (Å)"], data.at[j, "λ (Å)"], data.at[i, "λ (Å)"] / data.at[j, "λ (Å)"]))

      if len(r_es) > 0:
         ratios_e.append(r_es)

   # Use RMS of ratio errors.
   score = rms([o - e for o, e in zip(flatten_list(ratios), flatten_list(ratios_e))])
   
   errors.append( (score, ratios_e, comb) )

errors = sorted(errors, key=lambda x: x[0])

for i in range(0, 10):
   e = errors[i]
   
   print("RMS: %f" % (e[0]))

   print("Rest Ratios:")
   print(e[1])
   
   print("Elements:")
   
   comb = e[2]
   
   lambda_es = []
   
   for j in range(0, len(comb)):
      lambda_es.append(data.at[comb[j], "λ (Å)"])
      print("%s: %f -> %f" % (data.at[comb[j], "Ion"], data.at[comb[j], "λ (Å)"], os[j]))
      
   print("Redshift:")   
   print(statistics.mean([(o - e) / e for o, e in zip(os, lambda_es)]))
   
   print("")