#/usr/bin/env python
"""
Ein kleines Skript zum Berechnen der Saha Gleichung
"""

from math import pi, exp, sqrt
import sys
from numpy import logspace,loadtxt,ones
from Gnuplot import Gnuplot, GnuplotOpts, Data, funcutils
import Gnuplot

kB = 1.3806488e-23 # Boltzmannkonstante
h  = h= 6.62606957e-34 #Plancksches Wirkungsquantum
m_e = 9.10938291e-31  #Elektronenmase

def saha(n,T,E):
    return (2*pi*m_e*kB*T)**(1.5)/(n*h**3)*exp(-E/kB/T)

def solve_saha(n,T,E):
    a = saha(n,T,E)
    return a/(1+a)

if __name__ == '__main__':
    if len(sys.argv) != 4:
        print 'saha.py erwartet drei Argumente: n, T, E'
        print 'E in eV eingeben!'
    else:
        n = float(sys.argv[1])
        T = float(sys.argv[2])
        E = float(sys.argv[3])
        E = E*1.6e-19
        print 'Der Ionisationsgrad betraegt: ', solve_saha(n,T,E)
        saha_file = open('saha_New.dat','w')
        for T in logspace(2,5,1000):
            string = str(T) + ' ' + str(solve_saha(n,T,E)) + '\n'
            saha_file.write(string)
        saha_file.close()
            
        g = Gnuplot.Gnuplot()
        x,y = loadtxt('saha_New.dat',unpack=True)
        pldat = Data(x,y,with_='line lw 2 lc rgb "black"',title='Ionisationsgrad')
        plline = Data(x,ones(len(x)), with_='line lw 1 lc rgb "black"',title='')
        g('set zero 1.e-20')
        g('set size square')
        g('set xlabel "Temperatur [K]"')
        g('set ylabel "Ionisationsgrad"')
        g('set logscale xy')
        g('set yrange [1.e-150:1.e20]')
        g('set term postscript eps enhanced color "Arial" 24')
        g('set output "saha_new.eps"')
        g.plot(pldat,plline)