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("")
```