In [1]:
#Example of an Interactive program.
#I believe this is useful as a reminder computers are a great companion to research because they do as they are told.
#They are an aid not the answer.
#This program is used for the normalization and measurement of high resolution spectra,
#which is probably not important to all who will try this program
#but is a useful example in how to create a GUI to analyse data should you find the tools of your trade clunky or inadequate.
#The program takes two files. One fits file (it does not handle all wavelength fits headers, yet) and one MOOG formated line list.
#An example of both files are attached to this program.
#To use the program I recommend copying and pasting the program text into a file.
#Then you can run the program from there using the command "python file.py" where file.py is whatever you called the file. ENJOY!
#N.B. This is not the final version of this program so please do not use it as such. There are numerous errors in both calculations and cosmetics
In [1]:
import matplotlib
#opens a window on your computer
matplotlib.use('TkAgg')
from matplotlib.backends.backend_tkagg import FigureCanvasTkAgg, NavigationToolbar2TkAgg
#implement the deault mpl key bindings
from matplotlib.backend_bases import key_press_handler,MouseEvent
import tkMessageBox as box
import tkFileDialog as Tkf
import numpy as np
import pyfits
#Scipy for spline
from scipy.interpolate import interp1d
#scipy for integration
from scipy.integrate import simps
#Astropymodeling for gauss fit
from astropy.modeling import models, fitting
import matplotlib.pyplot as plt
import sys,os
if sys.version_info[0] < 3:
import Tkinter as Tk
else:
import tkinter as Tk
#dictionary to convert proton number into element symbol
edict = {}
edict['1'] = 'H'
edict['2'] = 'Hi'
edict['3'] = 'Li'
edict['4'] = 'Be'
edict['5'] = 'B'
edict['6'] = 'C'
edict['7'] = 'N'
edict['8'] = 'O'
edict['9'] = 'F'
edict['10'] = 'Ne'
edict['11'] = 'Na'
edict['12'] = 'Mg'
edict['13'] = 'Al'
edict['14'] = 'Si'
edict['15'] = 'P'
edict['16'] = 'S'
edict['17'] = 'Cl'
edict['18'] = 'Ar'
edict['19'] = 'K'
edict['20'] = 'Ca'
edict['21'] = 'Sc'
edict['22'] = 'Ti'
edict['23'] = 'V'
edict['24'] = 'Cr'
edict['25'] = 'Mn'
edict['26'] = 'Fe'
edict['27'] = 'Co'
edict['28'] = 'Ni'
edict['29'] = 'Cu'
edict['30'] = 'Zn'
edict['31'] = 'Ga'
edict['32'] = 'Ge'
edict['33'] = 'As'
edict['34'] = 'Se'
edict['35'] = 'Br'
edict['36'] = 'Kr'
edict['37'] = 'Rb'
edict['38'] = 'Sr'
edict['39'] = 'Y'
edict['40'] = 'Zr'
edict['41'] = 'Nb'
edict['42'] = 'Mo'
edict['43'] = 'Tc'
edict['44'] = 'Ru'
edict['45'] = 'Rh'
edict['46'] = 'Pd'
edict['47'] = 'Ag'
edict['48'] = 'Cd'
edict['49'] = 'In'
edict['50'] = 'Sn'
edict['51'] = 'Sb'
edict['52'] = 'Te'
edict['53'] = 'I'
edict['54'] = 'Xe'
edict['55'] = 'Cs'
edict['56'] = 'Ba'
edict['57'] = 'La'
edict['58'] = 'Ce'
edict['59'] = 'Pr'
edict['60'] = 'Nd'
edict['61'] = 'Pm'
edict['62'] = 'Sm'
edict['63'] = 'Eu'
edict['64'] = 'Gd'
edict['65'] = 'Tb'
edict['66'] = 'Dy'
edict['67'] = 'Ho'
edict['68'] = 'Er'
edict['69'] = 'Tm'
edict['70'] = 'Yb'
edict['71'] = 'Lu'
edict['72'] = 'Hf'
edict['73'] = 'Ta'
edict['74'] = 'W'
edict['75'] = 'Re'
edict['76'] = 'Os'
edict['77'] = 'Ir'
edict['78'] = 'Pt'
edict['79'] = 'Au'
edict['80'] = 'Hg'
edict['81'] = 'Tl'
edict['82'] = 'Pb'
edict['83'] = 'Bi'
edict['84'] = 'Bo'
edict['85'] = 'At'
edict['86'] = 'Rn'
edict['87'] = 'Tr'
edict['88'] = 'Ra'
edict['89'] = 'Ac'
edict['90'] = 'Th'
edict['91'] = 'Pa'
edict['92'] = 'U'
edict['93'] = 'Np'
edict['94'] = 'Pu'
edict['95'] = 'Am'
edict['96'] = 'Cm'
#dictionary for ionization level
idict = {}
idict['0'] = 'I'
idict['1'] = 'II'
idict['2'] = 'III'
idict['3'] = 'IV'
idict['4'] = 'V'
idict['5'] = 'VI'
idict['6'] = 'VII'
#Convert from air to vacuum wavelengths
def air_vac(wav0):
wav0 = wav0/(1.0 + 2.735182*10**(-4.0) + (131.4182/wav0**2.) + (2.76249*10.**8./wav0**4.))
return wav0
#correct wavelengths for rv
def rv_correct(lam0,rv):
cs = 2.99792458e5 #km/s
lam = lam0*(np.sqrt((1.+(rv/cs))/(1.-(rv/cs))))
return lam
#Line EW measurement edit popup
class ew_popup:
global root
def __init__(self,parent):
top = self.top = Tk.Toplevel(parent)
# self.myLabel = Tk.Label(top,text='EW Edit')
# self.myLabel.pack()
self.gButton = Tk.Button(top,text='Gauss Fit',command=self.sendg)
self.gButton.pack(side=Tk.LEFT)
self.gButton = Tk.Button(top,text='Simpson Rule Integration',command=self.sends)
self.gButton.pack(side=Tk.LEFT)
self.gButton = Tk.Button(top,text='Voigt Fit',command=self.sendv)
self.gButton.pack(side=Tk.LEFT)
self.gButton = Tk.Button(top,text='Delete',command=self.sendd)
self.gButton.pack(side=Tk.LEFT)
self.gButton = Tk.Button(top,text='Cancel',command=self.sendc)
self.gButton.pack(side=Tk.LEFT)
def sendd(self):
self.editparm = 'd'
self.top.destroy()
def sendg(self):
self.editparm = 'g'
self.top.destroy()
def sendc(self):
self.editparm = 'c'
self.top.destroy()
def sendv(self):
self.editparm = 'v'
self.top.destroy()
def sends(self):
self.editparm = 's'
self.top.destroy()
#Class for asking for order number to plot first
#No longer used switched to text box
class order_popup:
global root
def __init__(self,parent):
top = self.top = Tk.Toplevel(parent)
self.myLabel = Tk.Label(top,text='Order')
self.myLabel.pack()
self.myEntryBox = Tk.Entry(top)
self.myEntryBox.pack()
self.mySubmitButton = Tk.Button(top,text='Submit',command=self.send)
self.mySubmitButton.pack()
def send(self):
self.order = self.myEntryBox.get()
self.top.destroy()
#Class for asking for fits dimension with echelle data
#Only does this if echelle data not in zeroth dimension
class dimension_popup:
global root
def __init__(self,parent):
top = self.top = Tk.Toplevel(parent)
self.myLabel = Tk.Label(top,text='Fits Dimension')
self.myLabel.pack()
self.myEntryBox = Tk.Entry(top)
self.myEntryBox.pack()
self.mySubmitButton = Tk.Button(top,text='Submit',command=self.send)
self.mySubmitButton.pack()
def send(self):
self.order = self.myEntryBox.get()
self.top.destroy()
#Main gui create and modification loop
class gui_c(Tk.Frame):
def __init__(self,parent):
Tk.Frame.__init__(self,parent,background='white')
self.parent = parent
#Default setting for finding the fits dimension
self.fitsd = 0
#Intialize clicker
self.click = 0
#Initial gauss find tolerance
self.findtol = .1
#set default order
self.order = 1
#radial velocity corrections to line list
# self.rv = -269.16
# self.hrv = -2.14
#radial velocity corrections to line list
self.rv = 0.0
self.hrv = 0.0
#looking fore continuum
self.seeking = False
#line list input
self.linelistin = False
#EW have been measured
self.ewmeas = False
#Order has been measured
self.oewmeas = False
#Edit EW measurement
self.ewedit = False
#maximum line width in one direction for gauss fit
self.maxlinewidth = 2.
#Make sure file is saved before you destroy data
self.saved = True
#Variable for simpsons integration
self.sint = False
#Stellar parameters for Moog and Marcs someday
self.teff = 0.
self.feh = 0.
self.logg = 0.
self.vturb = 0.
#Start the creation of the window and GUI
self.centerWindow()
self.FigureWindow()
self.initUI()
#Create area and window for figure
def FigureWindow(self):
#Create the figure
self.f,self.a = plt.subplots(figsize=(8,10))
#Create window for the plot
self.canvas = FigureCanvasTkAgg(self.f,master=self)
#Draw the plot
self.canvas.draw()
#Turn on matplotlib widgets
self.canvas.get_tk_widget().pack(side=Tk.TOP,fill=Tk.BOTH,expand=1)
#Display matplotlib widgets
self.toolbar = NavigationToolbar2TkAgg(self.canvas,self)
self.toolbar.update()
self.canvas._tkcanvas.pack(side=Tk.TOP,fill=Tk.BOTH,expand=1)
#Connect mpl to mouse clicking
self.f.canvas.mpl_connect('button_press_event',self.on_click_event)
#Connect mpl to mouse clicking
self.f.canvas.mpl_connect('key_press_event',self.on_key_event)
#create button for initialization of Continuum normalization
specframe = Tk.Frame(self)
#Check box for vaccuum correction
self.checkvar = Tk.IntVar()
self.C1 = Tk.Checkbutton(master=specframe,text="Vaccum",variable=self.checkvar,onvalue=1,offvalue=0,command=self.vac_mod)
self.C1.pack(side=Tk.LEFT)
#create button to go up an order
upbutton = Tk.Button(master=specframe,text='Increase Order',command=self.increaseorder)
upbutton.pack(side=Tk.LEFT)
#create button to go down an order
downbutton = Tk.Button(master=specframe,text='Decrease Order',command=self.decreaseorder)
downbutton.pack(side=Tk.LEFT)
#Put in same frame as other spec information
downbutton = Tk.Button(master=specframe,text='Normalize Continuum',command=self.click_cont)
downbutton.pack(side=Tk.LEFT)
#Create button to ask for Moog formated line list
mooggrab = Tk.Button(master=specframe,text='Linelist',command=self.line_grab)
mooggrab.pack(side=Tk.LEFT)
specframe.pack(side=Tk.TOP)
#Create button to ask for Moog formated line list
moogsave = Tk.Button(master=specframe,text='Save EW',command=self.line_save)
moogsave.pack(side=Tk.RIGHT)
#Crate button to reset plot limits without clearing data
plotrbutton = Tk.Button(master=specframe,text='Reset Limits',command=self.reset_limits)
plotrbutton.pack(side=Tk.RIGHT)
#Crate button to reset plot limits without clearing data
openbutton = Tk.Button(master=specframe,text='Open File',command=self.open_file)
openbutton.pack(side=Tk.RIGHT)
#Pack the spec frame together
# frame2 = Tk.Frame(self)
specframe.pack(side=Tk.TOP)
#Resets the limits of the plot
def reset_limits(self):
self.a.set_xlim([self.wav0.min(),self.wav0.max()])
maxer = np.max(self.spec0[self.spec0.size/2-.1*self.spec0.size:self.spec0.size/2+.1*self.spec0.size])
self.a.set_ylim([0.,maxer*.1+maxer])
self.canvas.draw()
#Save Updated EW
def line_save(self):
self.linelistout = Tkf.asksaveasfile(mode='w')
if self.linelistout is None:
return
if self.checkvar.get() == 0:
#Write out line ew in MOOG format
for i,j in enumerate(self.line):
self.linelistout.write(' {0:10} {1:10} {2:10} {3:10} {5} {4:10}\n'.format(str(round(self.line[i],2)),str(self.lineel[i]),str(round(self.lineep[i],2)),str(round(self.linegf[i],2)),str(round(self.lineew[i],2)),' '*20))
#Always write out in air wavelengths
if self.checkvar.get() == 1:
#Write out line ew in MOOG format
for i,j in enumerate(self.aline):
self.linelistout.write(' {0:10} {1:10} {2:10} {3:10} {5} {4:10}\n'.format(str(round(self.aline[i],2)),str(self.lineel[i]),str(round(self.lineep[i],2)),str(round(self.linegf[i],2)),str(round(self.lineew[i],2)),' '*20))
#close the file
self.linelistout.close()
#Removes flag to warn user of unsave data
self.saved = True
#Grab Moog line list
def line_grab(self):
m = 0
#If the values already exist set the line ew measure to True
self.ewmeas = True
while m == 0:
self.linelist = self.get_file()
try:
fil = open(self.linelist,'r')
fils = fil.readlines()[1:]
line = []
lineel = []
lineels = []
lineep = []
linegf = []
lineew = []
for i in fils:
line.append(float(i[:10]))
lineel.append(i[10:20].replace(' ',''))
eltext = i[10:20].replace(' ','').split('.')
ion = idict[eltext[1]]
el = edict[eltext[0]]
lineels.append(el+ion)
lineep.append(float(i[20:30]))
linegf.append(float(i[30:40]))
#Try to grab EW from line list file
try:
lineew.append(round(float(i[50:70]),1))
#Append value -9999.9 if EW have not been measured
except IndexError:
#If the any values don't exist set the line ew measure to False
self.ewmeas = False
lineew.append(round(-9999.9,1))
m =1
except ValueError:
self.error = 7
self.onError()
#Convert the lists to numpy arrays and read them into a variable in the object
self.line = np.array(line)
self.lineel = np.array(lineel)
self.lineels = np.array(lineels)
self.lineep = np.array(lineep)
self.linegf = np.array(linegf)
self.lineew = np.array(lineew)
self.linelistin = True
self.text_loc()
self.plot_lines()
#Find location for labeling lines
def text_loc(self):
#If the spectrum is normalized just point linelist at normal y values
if np.median(self.spec0) < 10:
self.elloc = 1.1+np.zeros(self.line.size)
self.ewloc = 1.2+np.zeros(self.line.size)
#Else try half heartly to fit the continuum and scale
else:
forfit, = np.where(np.abs(self.spec0-np.median(self.spec0)) <= np.std(self.spec0)*3.)
input = rv_correct(self.line,self.rv-self.hrv)
specloc = np.interp(input,self.wav0[forfit],self.spec0[forfit])
self.elloc = specloc*.1+specloc
self.ewloc = specloc*.15+specloc
#Plot lines from the Moog line list
def plot_lines(self):
maxw = max(self.wav0)
minw = min(self.wav0)
findlines = rv_correct(self.line,self.rv-self.hrv)
use, = np.where((findlines >= minw) & (findlines <= maxw))
#locations of lines in order
#location on spectrum
self.olines = findlines[use]
#rest wavelengths of lines in order
self.orestlines = self.line[use]
#An array containing the text values to be updated when a manual measurement is applied
self.ewtext = []
for i in use:
self.a.text(rv_correct(self.line[i],self.rv-self.hrv),self.elloc[i],str(self.lineels[i]),rotation='vertical')
self.ewtext.append(self.a.text(rv_correct(self.line[i],self.rv-self.hrv),self.ewloc[i],str(self.lineew[i]),rotation='vertical'))
self.canvas.draw()
#Do wavelength conversion of linelist to vaccum
def vac_mod(self):
#release cursor from entry box and back to the figure
#needs to be done otherwise key strokes will not work
self.f.canvas._tkcanvas.focus_set()
self.aline = self.line.copy()
if ((self.checkvar.get() == 1) & (self.linelistin)):
self.line = air_vac(self.line)
self.spec_plot()
else:
#Uncheck box if linelist is not read in
if (self.linelistin) == False:
self.error = 8
self.onError()
#Uncorrect if mistakenly pressed
if (self.checkvar.get() == 0):
self.line = self.aline
self.spec_plot()
self.C1.deselect()
#Basic click event
def on_click_event(self,click):
#Click envents for continuum selection
if self.seeking:
self.contx.append(click.xdata)
self.conty.append(click.ydata)
self.pcontx.append(click.x)
self.pconty.append(click.y)
self.a.scatter(click.xdata,click.ydata,marker='x',color='red',zorder=8,s=35)
self.canvas.draw()
if self.ewedit :
self.contx.append(click.xdata)
self.conty.append(click.ydata)
self.mangaussfitlist.append(self.a.scatter(click.xdata,click.ydata,marker='x',color='red',zorder=8,s=35))
#Commands to inform the user what is available during line fitting
self.clickcommands.remove()
if len(self.contx) == 1:
self.clickcommands = self.a.text(self.olines[self.editewlineindex],self.ymarker,'Click Line Center')
if len(self.contx) == 2:
self.clickcommands = self.a.text(self.olines[self.editewlineindex],self.ymarker,'Left FWHM')
if len(self.contx) == 3:
self.clickcommands = self.a.text(self.olines[self.editewlineindex],self.ymarker,'Right FWHM')
if len(self.contx) == 4:
self.clickcommands = self.a.text(self.olines[self.editewlineindex],self.ymarker,'Right Continumm')
self.canvas.draw()
#Once all data points have been ented
if len(self.contx) == 5:
self.calc_gauss()
self.check_ew()
#Remove the X markers from the gaussian fit
for i in self.mangaussfitlist:
i.set_visible(False)
self.canvas.draw()
#Do simpsons rule integration
if self.sint:
self.contx.append(click.xdata)
self.conty.append(click.ydata)
self.intlist.append(self.a.scatter(click.xdata,click.ydata,marker='x',color='red',zorder=8,s=35))
self.clickcommands.remove()
if len(self.contx) == 1:
self.clickcommands = self.a.text(self.olines[self.editewlineindex],self.ymarker,'Right Side')
#Do the integration
if len(self.contx) == 2:
range, = np.where((self.wav0 >= np.min(self.contx)) & (self.wav0 <= np.max(self.contx)))
self.ewupdate = 1000.*simps(1.-self.spec0[range],x=self.wav0[range],even='avg')
self.canvas.draw()
#check EW
self.check_ew()
try:
self.fitlist[self.editewlineindex][0].remove()
except ValueError:
print 'line already removed'
self.fitlist[self.editewlineindex] = self.a.plot(self.wav0[range],self.spec0[range],color='red',zorder=8)
#Remove the X markers from the gaussian fit
for i in self.intlist:
i.set_visible(False)
self.canvas.draw()
#Function to check if the EW should be saved
def check_ew(self):
if box.askyesno('Save','Should the New EW be added? \n '+str(round(self.ewupdate,2))):
close = np.abs(self.orestlines[self.editewlineindex]-self.line)
update, = np.where(close == close.min())
self.lineew[update] = round(float(self.ewupdate),2)
#Plot the new EW
self.ewtext[self.editewlineindex].remove()
self.ewtext[self.editewlineindex] = self.a.text(rv_correct(self.line[update],self.rv-self.hrv),self.ewloc[update],str(self.lineew[update][0]),color='red',rotation='vertical')
self.canvas.draw()
#After the ew has been calculated remind use to save if they quit
self.saved = False
else:
self.fitlist[self.editewlineindex][0].remove()
#After the all the points have been clicked drop editing and calculate gaussian
self.ewedit = False
# self.spec_plot()
#Select points for continuum normalization
def click_cont(self):
if self.seeking == False:
self.seeking = True
self.contx = []
self.conty = []
#Pixel value so points can be deleted
self.pcontx = []
self.pconty = []
self.clickcommands = self.a.text(np.median(self.wav0),0.0,'Click on the Continuum then Click q to quit and fit \n d deletes bad points')
self.canvas.draw()
#Actually normalize the continuum
def cont_norm(self):
#First a second order legendre polynomail
#Changed to spline based on Gray Analysis of stellar photospheres
# lcoef = np.polynomial.legendre.legfit(self.contx,self.conty,len(self.contx)-1)
if len(self.contx) > 3:
#sort the continuum in x
sorted = np.argsort(np.array(self.contx))
#Spline interpolation of continuum points
fspline = interp1d(np.array(self.contx)[sorted],np.array(self.conty)[sorted],kind='cubic',bounds_error=False)
#Calc the continuum using spline
# self.cont = np.polynomial.legendre.legval(self.wav0,lcoef)
self.cont = fspline(self.wav0)
contplot = self.a.plot(self.wav0,self.cont,color='red',zorder=16)
self.canvas.draw()
#If the continuum normalization is good apply it
if box.askyesno("Continuum",'Is the Continuum Normalization Good?'):
self.spec.data[self.order-1] = self.spec0/self.cont
self.spec_set()
self.spec_plot()
else:
#Ask to keep points in array if continuum normalization is improper
if box.askyesno('Continuum','Keep Trying?'):
self.seeking = True
#remove pervious continuum
contplot[0].remove()
else:
self.spec_plot()
else:
#Error if not enough points for spline fit
self.error = 15
self.seeking = True
self.onError()
#Defind key strokes
def on_key_event(self,event):
#Exits continuum normalization and fits
if event.key == 'q':
self.seeking = False
self.cont_norm()
#Deletes continuum normalization point
if event.key == 'd':
if ((self.seeking) & (self.ewedit == False)) :
subx = np.array(self.pcontx)-event.x
suby = np.array(self.pconty)-event.y
close = np.sqrt(subx**2.+suby**2.)
delcol, = np.where(close == np.min(close))
delcol = delcol[0]
#Delete index with closest pixel value to where d was pressed
pointx = self.contx[delcol]
pointy = self.conty[delcol]
del self.conty[delcol]
del self.contx[delcol]
del self.pconty[delcol]
del self.pcontx[delcol]
self.a.scatter(pointx,pointy,marker='x',color='blue',zorder=15,s=35)
self.canvas.draw()
#key event for editing EW Measurement
if event.key == 'e':
#If the EW for a given measure have been measured
# if (((self.oewmeas[self.order-1]) | (self.ewmeas)) & (self.seeking == False)):
if (self.seeking == False):
subx = np.array(self.olines)-event.xdata
#Find the nearest point in wavelength space
self.editewlineindex, = np.where(np.abs(subx) == np.min(np.abs(subx)))
#If two points are the same distance away error
if self.editewlineindex.size > 1:
self.error = 12
else:
self.edit_ew()
else:
self.error = 11
self.onError()
#Another way to trigger continuum normalization
if event.key == 'c':
self.click_cont()
#Keystroke to reset limits
if event.key == 'r':
self.reset_limits()
#Another way to add and subtract orders to be similar to IRAF
if event.key == 'j':
self.decreaseorder()
#Another way to add and subtract orders to be similar to IRAF
if event.key == 'k':
self.increaseorder()
#automattically tries to fit gaussian to lines in line list
if event.key == 'a':
self.auto_gauss_fit()
#Edit EW measurement
def edit_ew(self):
#set EW range to focus on line
line = self.olines[self.editewlineindex]
self.a.set_xlim([line-5.,line+5.])
findmin, = np.where((self.wav0 >= line-5.) & (self.wav0 <= line+5.))
ymin = np.min(self.spec0[findmin])
self.a.set_ylim([ymin,1.2])
self.ymarker = np.abs(np.diff(self.a.get_ylim()))*.1+np.min(self.a.get_ylim())
self.canvas.draw()
choice = ew_popup(root)
root.wait_window(choice.top)
edit = choice.editparm
#manual gauss fit
if edit == 'g':
#points for the gaussian fit
self.mangaussfitlist = []
self.man_gauss_fit()
#replace line EW with -9999.9
if edit == 'd':
minfind = np.abs(self.orestlines[self.editewlineindex]-self.line)
remove, = np.where(minfind == minfind.min())
self.lineew[remove] = -9999.9
#Plot the new EW
self.ewtext[self.editewlineindex].remove()
self.ewtext[self.editewlineindex] = self.a.text(rv_correct(self.line[remove],self.rv-self.hrv),self.ewloc[remove],str(self.lineew[remove][0]),color='red',rotation='vertical')
#Remove poor fits
try:
self.fitlist[self.editewlineindex][0].remove()
except ValueError:
print 'line already removed'
self.canvas.draw()
#Return to normal part of spectrum
self.reset_limits()
#Do straight integration
if edit == 's':
self.intlist = []
self.man_int()
#Do not edit if c go straight to spec_plot
if edit == 'c':
self.reset_limits()
#Do a Voigt fit Not yet implimented
if edit == 'v':
self.error = 9
self.onError()
#Return to normal part of spectrum
self.reset_limits()
#Calculate a Gaussian from manual fit
def calc_gauss(self):
#First element is the continuum
cont = (self.conty[0]+self.conty[-1])/2.
#Second element is line center
depth = cont-self.conty[1]
#Width of feature (estimate of 2simga
width = np.abs(self.contx[2]-self.contx[3])
g_model = models.Gaussian1D(amplitude=depth,mean=self.contx[1],stddev=width/2.)
fit_g = fitting.LevMarLSQFitter()
#Use input parameters to find LSQ fit to 3 sigma of line
use, = np.where((self.wav0 >= g_model.mean-3*(width/2.)) & (self.wav0 <= g_model.mean+3.*(width/2.)))
gauss_out = fit_g(g_model,self.wav0[use],1.-self.spec0[use])
#Read EW in mA
ew = 1000.*gauss_out.amplitude.value*1.06447*2.*np.abs(gauss_out.stddev.value)
#Plot found EW
sig3 = np.linspace(self.contx[1]-4.*gauss_out.stddev.value,self.contx[1]+4.*gauss_out.stddev.value,100)
#Put gaussian plots in a list so they can be deleted if unsatisfactory
#First remove old plot
#need the zero since for whatever reason plots are different than text
try:
self.fitlist[self.editewlineindex][0].remove()
except ValueError:
print 'No previous measure'
#Now plot new gaussian in place
self.fitlist[self.editewlineindex] = self.a.plot(np.array(sig3),1.-gauss_out(sig3),color='red',zorder=3)
self.a.text(self.contx[1],.5,str(round(ew,1)),color='red',rotation='vertical')
self.canvas.draw()
# the EW to pass to the plotting routine
self.ewupdate = round(ew,1)
#Simpson's Rule integration
def man_int(self):
self.sint = True
self.contx = []
self.conty = []
self.clickcommands = self.a.text(self.olines[self.editewlineindex],self.ymarker,'Click Left intregation point')
self.canvas.draw()
#Command to manually fit Gaussian
def man_gauss_fit(self):
self.ewedit = True
self.contx = []
self.conty = []
self.clickcommands = self.a.text(self.olines[self.editewlineindex],self.ymarker,'Click Left Continuum')
self.canvas.draw()
#Command to fit gaussians to lines in given order
def auto_gauss_fit(self):
if self.linelistin:
#set EW measured flag for order to be true
self.oewmeas[self.order-1] = True
#List for gaussian plots
self.fitlist = []
#Use the min and maximum wavlength values in order to search for lines
maxw = max(self.wav0)
minw = min(self.wav0)
findlines = rv_correct(self.line,self.rv-self.hrv)
use, = np.where((findlines >= minw) & (findlines <= maxw))
#array for newly calculated ew
newew = []
for i in use:
# self.a.set_xlim([findlines[i]-5.,findlines[i]+5.])
#look for minimum value with some tolerance value of line center (initally setting to .1)
lookbot, = np.where(np.abs(self.wav0-findlines[i]) <= self.findtol)
minspec = self.spec0[lookbot].min()
lookbot, = np.where(self.spec0 == minspec)
#Only use line if minimum is found
if lookbot.size > 0:
minwav = self.wav0[np.where(self.spec0 == minspec)][0]
#Find the line depth and half max
depthspec = np.float(1.-minspec)
halfspec = 1.-(depthspec*.5)
#find the blue line side information
#If more than one minimum pick the first
bluefind = lookbot[0]
listwav = []
listspec = []
distance = 0.
#look for the half max unless it is larger than 2 Angstroms in width
while ((self.spec0[bluefind] <= halfspec) & (distance <= self.maxlinewidth) & (bluefind >= 0)):
bluefind = bluefind-1
listwav.append(np.float(self.wav0[bluefind]))
listspec.append(1.-self.spec0[bluefind])
distance = minwav-self.wav0[bluefind]
#reverse the blue wavelengths so they are in the correct order
listwav.reverse()
listspec.reverse()
#append the center of the line to the list
listwav.append(np.float(minwav))
listspec.append(depthspec)
#find the red line side information
redfind = lookbot[0]
distance = 0.
#look for the half max unless it is larger than 2 Angstroms in width
while ((self.spec0[redfind] <= halfspec) & (distance <= self.maxlinewidth) & (redfind < self.spec0.size-2)):
redfind = redfind+1
listwav.append(np.float(self.wav0[redfind]))
listspec.append(1.-self.spec0[redfind])
distance = self.wav0[redfind]-minwav
#Only include lines with greater than 3 pixels in the line
#and the sigma is greater than 1 pixel seperation in wavelength
redhm = np.mean(listwav[-2:])
bluehm = np.mean(listwav[:2])
sigmaest = (redhm-bluehm)/2.
if ((len(listspec) > 2) & (sigmaest >= np.diff(self.wav0).min())):
#Fit a gaussian model to EW
g_model = models.Gaussian1D(amplitude=depthspec,mean=minwav,stddev=sigmaest)
fit_g = fitting.LevMarLSQFitter()
#Use input parameters to find LSQ fit to 3 sigma of line
useg, = np.where((self.wav0 >= g_model.mean-3*(sigmaest)) & (self.wav0 <= g_model.mean+3.*(sigmaest)))
gauss_out = fit_g(g_model,self.wav0[useg],1.-self.spec0[useg])
# gauss_out = fit_g(g_model,np.array(listwav),np.array(listspec))
#Read EW in mA
ew = 1000.*gauss_out.amplitude.value*1.06447*2.*np.abs(gauss_out.stddev.value)
if ew < 200.:
#Append new EW to list
newew.append(round(ew,1))
#Plot found EW
sig3 = np.linspace(minwav-4.*gauss_out.stddev.value,minwav+4.*gauss_out.stddev.value,100)
self.fitlist.append(self.a.plot(np.array(sig3),1.-gauss_out(sig3),color='red',zorder=3))
#Up
self.a.text(findlines[i],.6,str(round(ew,1)),color='red',rotation='vertical')
else:
newew.append(round(-9999.9,1))
self.fitlist.append(self.a.plot(np.arange(2),np.arange(2),color='red',zorder=3))
#need the zero since for whatever reason plots are different than text
self.fitlist[-1][0].remove()
self.a.text(findlines[i],.5,str(-9999.9),color='red',rotation='vertical')
else:
newew.append(-9999.9)
#Create the plot so everything remains in the same order but immediately drop the plot from the window
self.fitlist.append(self.a.plot(np.arange(2),np.arange(2),color='red',zorder=3))
#need the zero since for whatever reason plots are different than text
self.fitlist[-1][0].remove()
self.a.text(findlines[i],.2,str(-9999.9),color='red',rotation='vertical')
# print listspec
# print listwav
# print 'No line found'
#draw found EW
else:
newew.append(-9999.9)
#Create the plot so everything remains in the same order but immediately drop the plot from the window
self.fitlist.append(self.a.plot(np.arange(2),np.arange(2),color='red',zorder=3))
#need the zero since for whatever reason plots are different than text
self.fitlist[-1][0].remove()
self.a.text(findlines[i],.4,str(-9999.9),color='red',rotation='vertical')
self.canvas.draw()
#Ask if Line EW should be updated
if box.askyesno('Update','Should EW be Updated?'):
#Update lines in order
for j,i in enumerate(use):
self.lineew[i] = newew[j]
#Plot the new EW
self.ewtext[j].remove()
self.ewtext[j] = self.a.text(findlines[i],self.ewloc[i],str(self.lineew[i]),color='red',rotation='vertical')
self.canvas.draw()
#Remind the used the changes have not been saved if they try to exit
self.saved = False
# self.spec_plot()
else:
self.error = 8
self.onError()
self.canvas.draw()
#Command to increase the order to plot new spectrum
def increaseorder(self):
self.order = self.order+1
if self.order > self.ordermax:
self.order = 1
self.sorder.set(str(int(self.order)))
self.a.clear()
self.spec_set()
self.spec_plot()
#Command to decrease order to plot new spectrum
def decreaseorder(self):
self.order = self.order-1
if self.order < 1:
self.order = self.ordermax
self.sorder.set(str(int(self.order)))
self.a.clear()
self.spec_set()
self.spec_plot()
#Create window in center of screen
def centerWindow(self):
w = 2000
h = 1200
sw = self.parent.winfo_screenwidth()
sh = self.parent.winfo_screenheight()
x = (sw-w)/2
y = (sh-h)/2
self.parent.geometry('%dx%d+%d+%d' % (w,h,x,y))
#Intialize the GUI
def initUI(self):
self.parent.title("Echelle Spectrum Analysis")
#create frame for plotting
frame = Tk.Frame(self,relief=Tk.RAISED,borderwidth=1)
frame.pack(fill=Tk.BOTH,expand=1)
self.pack(fill=Tk.BOTH,expand=1)
#set up okay and quit buttons
quitButton = Tk.Button(self,text="Quit",command=self.onExit)
quitButton.pack(side=Tk.RIGHT,padx=5,pady=5)
saveButton = Tk.Button(self,text="Save Fits",command=self.save_fits)
saveButton.pack(side=Tk.RIGHT)
#set up velocity information to be added to line list
rvText = Tk.StringVar()
rvText.set("Radial Velocity (km/s)")
rvDir = Tk.Label(self,textvariable=rvText,height=4)
rvDir.pack(side=Tk.LEFT)
#Add so Velocity can be updated
rv = Tk.StringVar()
rv.set(str(round(self.rv,2)))
self.rvval = Tk.Entry(self,textvariable=rv,width=10)
self.rvval.bind("<Return>",self.rv_callback)
self.rvval.pack(side=Tk.LEFT,padx=5,pady=5)
#set up heliocentric velocity information to be added to line list
hrvText = Tk.StringVar()
hrvText.set("Heliocentric Correction (km/s)")
hrvDir = Tk.Label(self,textvariable=hrvText,height=4)
hrvDir.pack(side=Tk.LEFT)
#Add so Velocity can be updated
hrv = Tk.StringVar()
hrv.set(str(round(self.hrv,2)))
self.hrvval = Tk.Entry(self,textvariable=hrv,width=10)
self.hrvval.bind("<Return>",self.rv_callback)
self.hrvval.pack(side=Tk.LEFT,padx=5,pady=5)
#set up effective temperature
teffText = Tk.StringVar()
teffText.set("Teff (K)")
teffDir = Tk.Label(self,textvariable=teffText,height=4)
teffDir.pack(side=Tk.LEFT)
#Add so effective temperature can be updated
teff = Tk.StringVar()
teff.set(str(round(self.teff,2)))
self.teffval = Tk.Entry(self,textvariable=teff,width=10)
self.teffval.bind("<Return>",self.stellar_param)
self.teffval.pack(side=Tk.LEFT,padx=5,pady=5)
#set up [Fe/H]
fehText = Tk.StringVar()
fehText.set("[Fe/H]")
fehDir = Tk.Label(self,textvariable=fehText,height=4)
fehDir.pack(side=Tk.LEFT)
#Add so [Fe/H] can be updated
feh = Tk.StringVar()
feh.set(str(round(self.feh,2)))
self.fehval = Tk.Entry(self,textvariable=feh,width=10)
self.fehval.bind("<Return>",self.stellar_param)
self.fehval.pack(side=Tk.LEFT,padx=5,pady=5)
#set up log(g)
loggText = Tk.StringVar()
loggText.set("log(g)")
loggDir = Tk.Label(self,textvariable=loggText,height=4)
loggDir.pack(side=Tk.LEFT)
#Add so log(g) can be updated
logg = Tk.StringVar()
logg.set(str(round(self.logg,2)))
self.loggval = Tk.Entry(self,textvariable=logg,width=10)
self.loggval.bind("<Return>",self.stellar_param)
self.loggval.pack(side=Tk.LEFT,padx=5,pady=5)
#set up microturbulent velocity
vturbText = Tk.StringVar()
vturbText.set("Vturb (km/m)")
vturbDir = Tk.Label(self,textvariable=vturbText,height=4)
vturbDir.pack(side=Tk.LEFT)
#Add so microturbulent velocity can be updated
vturb = Tk.StringVar()
vturb.set(str(round(self.vturb,2)))
self.vturbval = Tk.Entry(self,textvariable=vturb,width=10)
self.vturbval.bind("<Return>",self.stellar_param)
self.vturbval.pack(side=Tk.LEFT,padx=5,pady=5)
#set up order number
orderText = Tk.StringVar()
orderText.set("Order")
orderDir = Tk.Label(self,textvariable=orderText,height=4)
orderDir.pack(side=Tk.LEFT)
#Add so order number can be updated
self.sorder = Tk.StringVar()
self.sorder.set(str(int(self.order)))
self.orderval = Tk.Entry(self,textvariable=self.sorder,width=10)
self.orderval.bind("<Return>",self.on_order_box)
self.orderval.pack(side=Tk.LEFT,padx=5,pady=5)
#set up popup menu
self.menu = Tk.Menu(self.parent,tearoff=0)
self.menu.add_command(label="Beep",command=self.bell())
self.menu.add_command(label="Exit",command=self.onExit)
self.parent.bind("<Button-3>",self.showMenu)
self.pack()
#set up Submenu
menubar = Tk.Menu(self.parent)
self.parent.config(menu=menubar)
fileMenu = Tk.Menu(menubar)
subMenu = Tk.Menu(fileMenu)
subMenu.add_command(label="New File",command=self.open_file)
fileMenu.add_cascade(label="Import",menu=subMenu,underline=0)
#create another item in menu
fileMenu.add_separator()
fileMenu.add_command(label='Exit',underline=0,command=self.onExit)
menubar.add_cascade(label="File",underline=0,menu=fileMenu)
self.open_file()
def open_file(self):
#ask to open fits file
self.filename = self.get_file()
#open the fits file
self.fitfile = self.check_spec()
#number of dimension of fits file
self.dmax = len(self.fitfile)
#use the fits dimension with the spectrum in it (default is 0)
m = 0
while m == 0:
try:
self.ordermax = self.fitfile[self.fitsd].shape[0]
m = 1
except IndexError:
#ask for fits dimension with ordered spectrum inside
self.error = 2
self.onError()
self.fitsd = self.get_fitsd()
#set up boolean array with whether order has been measured or not
self.oewmeas = np.ones(self.ordermax,dtype=bool)
#Set the default pyfits information to be extracted
self.spec = self.fitfile[self.fitsd]
#starts plotting the fits file
#First ask for order
#Changed to text box
# self.order = self.on_order()
#set standard header
self.standard_header()
#Set up the infromation need to plot the spectrum
self.spec_set()
#Draw said spectrum
self.spec_plot()
#RV Callback corrects line list for radial velocities
def rv_callback(self,event):
#release cursor from entry box and back to the figure
#needs to be done otherwise key strokes will not work
self.f.canvas._tkcanvas.focus_set()
if self.linelistin:
self.rv = float(self.rvval.get())
self.hrv = float(self.hrvval.get())
self.spec_plot()
self.text_loc()
self.plot_lines()
else:
self.error = 8
self.onError()
def stellar_param(self,event):
#release cursor from entry box and back to the figure
#needs to be done otherwise key strokes will not work
self.f.canvas._tkcanvas.focus_set()
try:
self.teff = float(self.teffval.get())
self.feh = float(self.fehval.get())
self.logg = float(self.loggval.get())
self.vturb = float(self.vturbval.get())
#error if not floats
except ValueError:
self.error = 10
self.onError()
#save fits file
def save_fits(self):
self.fileout = Tkf.asksaveasfilename(filetypes=(('FITS File','*.fits'),('FITS File','*.FITS'),('FITS File','*.fit'),('FITS File','*.FIT')))
#cannot overwrite current file or core will dump
if self.fileout == self.fitfile.filename():
self.error = 14
self.onError()
else:
try:
self.spec.writeto(self.fileout,output_verify='ignore')
except IOError:
if box.askyesno("Over write file?","Do you want to over write "+self.fileout+"?"):
self.spec.writeto(self.fileout,clobber=True)
#Function to get the dimension of
def get_fitsd(self):
m = 0
while m == 0:
try:
inputd = dimension_popup(root)
root.wait_window(inputd.top)
dim = int(inputd.order)
if ((dim > 0) & (dim <= self.dmax)):
m = 1
else:
#Error order is out of range
self.error = 5
self.onError()
except ValueError:
#Error order is not an integer
self.error = 4
self.onError()
return dim
#Function for retrieving order from popup
def on_order(self):
m = 0
while m == 0:
try:
inputO = order_popup(root)
root.wait_window(inputO.top)
order = int(inputO.order)
if ((order > 0) & (order <= self.ordermax)):
m = 1
else:
#Error order is out of range
self.error = 3
self.onError()
except ValueError:
#Error order is not an integer
self.error = 4
self.onError()
return order
#Function for retrieving order from text box
def on_order_box(self,event):
#release cursor from entry box and back to the figure
#needs to be done otherwise key strokes will not work
self.f.canvas._tkcanvas.focus_set()
m = 0
while m == 0:
try:
order = self.orderval.get()
order = int(order)
if ((order > 0) & (order <= self.ordermax)):
m = 1
self.order = order
self.spec_set()
self.spec_plot()
else:
#Error order is out of range
self.error = 3
self.onError()
except ValueError:
#Error order is not an integer
self.error = 4
self.onError()
#Shows a menu with right click
def showMenu(self,e):
self.menu.post(e.x_root,e.y_root)
#Exits the program
def onExit(self):
if self.saved == False:
if box.askyesno('Warning','You have unsaved EW Data!\n Are you sure you want to Exit?'):
self.quit()
else:
self.line_save()
else:
self.quit()
#Tells Why Order information is incorrect
def onError(self):
if self.error == 1:
box.showerror("Error","File Not Found")
if self.error == 2:
box.showerror("Error","Fits Dimension Does not contain proper ordered Spectrum")
if self.error == 5:
box.showerror("Error","Fits Dimension is Out of range (0-"+str(self.dmax)+")")
if self.error == 3:
box.showerror("Error","Order Value is out of Range (1-"+str(self.ordermax)+')')
if self.error == 4:
box.showerror("Error","Value Must be an Integer")
if self.error == 6:
box.showerror("Error","File is not in Fits Format")
if self.error == 7:
box.showerror("Error","File is not formated for Moog")
if self.error == 8:
box.showerror("Error","No line list has been read")
if self.error == 9:
box.showerror("Error","Voigt Function is Currently in Developement")
if self.error == 10:
box.showerror("Error","Value Must be Float")
if self.error == 11:
box.showerror("Error","EW Have not been Measured (press a)")
if self.error == 12:
box.showerror("Error","Could not determine line desired for editing \n Please try again")
if self.error == 14:
box.showerror("Error","Cannot overwrite currently opened fits file \n Please choose a different file name")
if self.error == 15:
box.showerror("Error","Not enough points in continuum \n Please try again")
#routine for opening file
#check to make sure file is found
def get_file(self):
filename = Tkf.askopenfilename()
while os.path.isfile(filename) == False:
self.error = 1
self.onError()
filename = Tkf.askopenfilename()
return filename
#make sure file is a fits file
def check_spec(self):
ffile = 0
while ffile == 0:
try:
fitfile = pyfits.open(self.filename)
ffile = 1
except IOError:
self.error = 6
self.onError()
self.filename = self.get_file()
return fitfile
#Uses IRAF WAT Headers
def standard_header(self):
text = ''
for i in self.spec.header["WAT2*"].itervalues():
add = i
#Correct for any WAT Values less than 68 Characters
while len(add) < 68:
add = add+' '
text = text+add
#Create an array containing one order of information per array
# text = text.replace('spe c','spec').replace('sp ec','spec').replace('s pec','spec')
breakt = text.replace('multipspec ','').split('spec')[2:]
#loop over all orders creating a standard header which will be able to be read
for i,j in enumerate(breakt):
split = j.split(' = ')
astr = split[0]
while len(astr) < 2:
astr = '0'+astr
wavd = split[1].split(' ')
try:
lam0 = np.float(wavd[3])
dlam = np.float(wavd[4])
self.spec.header['CRVL1_'+astr] = lam0
self.spec.header['CDLT1_'+astr] = dlam
except ValueError:
print 'Trying any way'
#set variables spectrum of a given order
def spec_set(self):
ostring = str(self.order)
while len(ostring) < 2:
ostring = '0'+ostring
self.dlam = self.spec.header['CDLT1_'+ostring]
self.lam0 = self.spec.header['CRVL1_'+ostring]
self.wav0 = np.arange(self.spec.data[self.order-1].size).astype(float)*self.dlam+self.lam0
self.spec0 = self.spec.data[self.order-1]
#actually plot the spectrum
def spec_plot(self):
self.a.clear()
self.a.plot(self.wav0,self.spec0,color='black',zorder=1)
self.a.set_xlabel(r'Wavelength ($\AA$)')
self.a.set_ylabel(r'Counts')
self.a.set_xlim([min(self.wav0),max(self.wav0)])
self.canvas.draw()
if self.linelistin:
self.text_loc()
self.plot_lines()
#Starts the main loop of the program
def main():
#Main loop
global root
root = Tk.Tk()
app = gui_c(root)
root.mainloop()
if __name__=="__main__":
#create root frame
main()
In []: