Skip to content

Instantly share code, notes, and snippets.

@thriveth
Last active September 9, 2016 18:04
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save thriveth/6323a65869f3a1911013b88454d00a46 to your computer and use it in GitHub Desktop.
Save thriveth/6323a65869f3a1911013b88454d00a46 to your computer and use it in GitHub Desktop.
Python, Chaco, TRaits(UI) -based interactive Voigt profile editor. Good to play around with to understand how b and N affect the profile.
#!/usr/bin/env python
# encoding: utf-8
# Copyright 2013, 2016 T. Emil Rivera-Thorsen
# trive@astro.su.se
''' To run: python vp-interactive-py
'''
import scipy as sp
import scipy.constants as con
from traits.api import *
from traitsui.api import View, RevertButton, OKButton, CancelButton, Item
from chaco.chaco_plot_editor import ChacoPlotItem
# cgs units (required in T-G voigt implementation):
cgsc = 2.998e10
cgsme = 9.1095e-28
cgse = 4.8032e-10
# Voigt-Hjerting function, as approximate by Tepper-Garcia 2006:
def H(a, x):
""" The H(a, u) function of Tepper Garcia 2006
"""
P = x ** 2
H0 = sp.e ** (-(x ** 2))
Q = 1.5 / x ** 2
H = H0 - a / sp.sqrt(sp.pi) / P * \
(H0 * H0 * (4 * P * P + 7 * P + 4 + Q) - Q - 1)
return H
# Class that handles user interacting and calculations that are best done
# inside said class. Plus shows the outpot.
# Hooray for the Enthought Tool Suite!
class VoigtEdit(HasTraits):
"""The interface to edit the thingie.
"""
pos = Range(800., 25000.)
z = Range(0., 10.)
N = Range(1.e6, 5e22)
B = Range(1., 1000) # In cgs. Must be conv to km/s before passed on.
b = Property(Float, depends_on=['B']) # Holds the cm/s version of B.
f = Float
g = Float
# The rest are properties of the above Traits.
vp = Property(Array, depends_on=['C', 'a', 'N', 'x'])
# or Ca instead of c and a
wl = Property(Array, depends_on=['pos'])
x = Property(Array, depends_on=['pos', 'dwld'])
C = Property(Float, depends_on=['f', 'g'])
wavlen = Property(Array, depends_on=['wl', 'z'])
dwld = Property(Float, depends_on=['b', 'wl0'])
a = Property(Float, depends_on=['wl0', 'g', 'b'])
wl0 = Property(Float, depends_on=['pos'])
Ca = Property(Float, depends_on=['f', 'b', 'wl0'])
fake_x = Property(Array, depends_on=['x'])
def _get_b(self):
return self.B * 10000. # km/s -> cm/s
def _pos_default(self):
return 1215.67 # rest wl of line center in Å
def _N_default(self):
return 14.5e20
def _b_default(self):
return 20.
def _z_default(self):
return 0.
def _f_default(self):
f = 0.416 # Write different getter if other than Lyman Alpha is wanted.
return f
def _g_default(self):
g = 6.265e8 # Write different getter if other than Lyman Alpha is wanted.
return g
def _get_wl0(self):
return self.pos * 1.e-10
def _get_wl(self):
return sp.linspace(self.pos - 100, self.pos + 100, 50000) * 1.e-10
def _get_dwld(self):
return self.b * self.wl0 / con.c
def _get_x(self):
return (self.wl - self.wl0) / self.dwld
def _get_C(self):
return 4 * sp.pi * cgse * self.f / (cgsme * cgsc * self.g)
def _get_Ca(self):
Ca = cgse ** 2 / (cgsme * cgsc) * sp.sqrt(sp.pi ** 3) / sp.pi * self.f \
* self.wl0 / self.b
return Ca
def _get_a(self):
return self.wl0 * self.g / (4 * sp.pi * self.b)
def _get_wavlen(self):
return self.wl * 1.e10 * (1. + self.z)
def _get_fake_x(self):
return self.x * 1.e10
def _get_vp(self):
return sp.exp(-self.Ca * self.N * H(self.a, self.x))
view = View(
ChacoPlotItem(
'wavlen', 'vp', y_bounds=(-.2, 1.2),
y_auto=False,
show_label=False, x_label='wavelength',
width=750, height=500,
y_label='relative flux'
),
Item('z'),
Item('N', label='N, cm^-2'),
Item('B', label='b, km/s'),
resizable=True,
title='Voigt Profile editor',
buttons=[RevertButton, CancelButton, OKButton]
)
# Just a little demo.
if __name__ == '__main__':
v = VoigtEdit()
v.configure_traits()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment