#/usr/bin/env python
"""
Ein kleines Skript zum Berechnen der Saha Gleichung
"""
from pylab import *

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!'
        print 'n in [1/m^3], T in [K]'
    else:
        n = float(sys.argv[1])
        T = float(sys.argv[2])
        E = float(sys.argv[3])
        E = E*1.6e-19
        print 'Sahagleichung ergibt:', saha(n,T,E)
        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()

        x,y = loadtxt('saha_New.dat',unpack=True)
        fig, ax = subplots()
        fig.set_size_inches(8, 8)
        ax.tick_params(axis='x',labelsize=16)
        ax.tick_params(axis='y',labelsize=16)
        ax.set_xlabel('Temperatur [K]',fontsize=16)
        ax.set_ylabel('Ionisationsgrad',fontsize=16)
        ax.tick_params(which='major',length=8,width=1)
        ax.tick_params(which='minor',length=4,width=1)
        ax.set_xscale("log",nonposy='clip')
        ax.set_yscale("log",nonposy='clip')
        ax.plot(x,y,'k-',linewidth=2.5)
        ax.set_ylim([1.e-300,1.e20])
        a=ax.get_yticks()
        print a
        ax.yaxis.set_ticks(logspace(-300,20,17),('a'))
        ax.set_yticks(logspace(-300,20,17))
        print a
        #x0,x1 = ax.get_xlim()
        #y0,y1 = ax.get_ylim()
        #ax.set_aspect('equal')
        #ax.set(adjustable='box-forced', aspect='equal')
        savefig('saha_new.eps')
        savefig('saha_new.png')
        savefig('saha_new.pdf')
        show()

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