Creating and manipulating FITS files

I wrote a program that manipulated data from VLA observations and created FITS files to be used in creating spectral energy distributions of high redshift radio galaxies. I really tried to be as efficient as possible. I made lists of the variables and lists I had to use. I looped everything where I could, but it still came out horribly convoluted.

```
#!/usr/bin/python

import math
from astropy.io import ascii
import astropy.cosmology as cosmo
import pyfits as fits
import matplotlib.pyplot as plt
from matplotlib import lines
import numpy as np
import os

# removing any duplicates from the current directory
try:
os.system("rm -f sed_and_alpha.png")
except:
pass
try:
os.system("rm -f sed_and_alpha.fits")
except:
pass
#constants
CBAND = 4.86*10**9
XBAND = 8.46*10**9
CF = (10**-26)

#functions
def get_luminosity_distance(redshift):
'''finds luminosity distance'''
lum_dist = cosmo.luminosity_distance(redshift)
tmp = str(lum_dist)
tmp = tmp.replace("Mpc","")
tmp = float(tmp)
return tmp3.0910**22
print lum_dist

def get_alpha(s1,s2,v1,v2):
'''gets alpha from two flux densities'''
top = math.log10(s1/s2)
bottom = math.log10(v2/v1)
return top/bottom

def get_luminosity(D_l,z,a,nu1,nu2,Sv2):
'''gets luminosity from alpha'''
Sv2 = (Sv2/1000)*CF
lum = (4math.pi((D_l2)))(1/((1+z)(1+a)))((nu1/nu2)**a)*(Sv2)
return lum

#lists
available_redshifts = [] # For making sure there are no blank spaces
luminosity_distribution = [] # which calculated luminosities are available
alpha_list = [] # available alphas
new_luminosity_distribution = [] # luminosity distr. after str[-2:]
default_data = {}
list22 = []
list23 = []
list24 = []
list25 = []
list26 = []
list27 = []
list28 = []
exponent_list = [
list22,
list23,
list24,
list25,
list26,
list27,
list28,

The first, really simple thing I would do is tidy up the way you're using scientific notation. Rather than writing 4.86*10**9, you can write 4.86e9, and so on. It's easier to read and less error-prone.

Second, some of the things you're trying to do are duplicated in library functions. For example, your get_luminosity_distance function can be replaced with cosmo.luminosity_distance(redshift).si.value. (at least, in the current version of astropy. I don't know which version you're using).

I would store your exponent_lists in a dictionary; you're trying to do that already with the exponent_list_names list. To check the range of exponent values, you could then do something like:

for idx, luminosity in enumerate(new_luminosity_distribution):
    exponent_dict[int(np.log10(luminosity))].append(idx)

which is much more compact than

# checking the range of exponent values
for i in range(len(new_luminosity_distribution)-1):
    tmp = str(new_luminosity_distribution[i])
    tmp = tmp[-2:]
    if int(tmp) == 22:
        list22.append(i)
    elif int(tmp) == 23:
        list23.append(i)
    elif int(tmp) == 24:
        list24.append(i)
    elif int(tmp) == 25:
        list25.append(i)
    elif int(tmp) == 26:
        list26.append(i)
    elif int(tmp) == 27:
        list27.append(i)
    elif int(tmp) == 28:
        list28.append(i)

That actually brings up another issue: when you're looping over the elements of a list, rather than doing

for i in range(len(thelist)):
    ... etc. ...

it's much more pythonic to do:

for item in thelist:
    ... etc. ...

if you just need the items, or use enumerate (as I did above) if you need the index, as well.

StackExchange Code Review Q#119401, answer score: 3