#!/usr/bin/python -u

#
# @file theory.py
# @purpose Do theory stuff
#		Because sadly, Mathematica does not work on my new laptop
# @author Sam Moore
# @date August 2012
#

import sys
import os
import math
import re # Regular expressions - for removing comments
#import odict #ordered dictionary

import Gnuplot, Gnuplot.funcutils

gnuplot = Gnuplot.Gnuplot()


#
# @function se_dist
# @purpose Model the secondary electron distribution using M. Furman's formula
#
def se_dist(s, Ep, sigma_p, dE, n_max=0.75):
	results = []
	E = 0.0
	while (E < 8.0):
		results.append([E, n_max *s * E / (s - 1.0 + E**s) , math.exp(-(E - Ep)**2.0 / (2.0*sigma_p**2.0))])
		E += dE
	return results

#
# @function integrate
# @purpose integrate function numerically over a range of arguments
# 
def integrate(f, xmin, xmax, dx):
	x = xmin
	result = 0.00
	while (x <= xmax):
		result += f(x) * dx
		x += dx
	return result

# @function derivative
# @purpose get derivative of a function at a value
def der(f, x, dx):
	return (f(x + dx) - f(x)) / dx

# 
# @function tcs
# @purpose Model S(E)
def tcs(f, sigma, Emin, Emax, dE):
	results = []
	E = Emin
	while (E < Emax):
		results.append([E, (1.0 - sigma(0)) * f(-E) - integrate(lambda e : f(e - E) * sigma(E), Emin, Emax, dE)])
		E += dE
	return results

def delta(x):
	if (x == 0):
		return 1.0
	else:
		return 0.0

def table(f, xmin, xmax, dx):
	result = []
	x = xmin
	while (x <= xmax):
		result.append([x, f(x)])
		x += dx
	return result

def gaussian(x, sigma):
	return math.exp(- (x**2.0)/(2.0 * sigma**2.0)) / (sigma * (2.0 * math.pi)**0.50)

def step(x, sigma, T):
	return 1.0 / (math.exp((x - sigma)/T) + 1.0)

# 
# @function write_data
# @purpose Write a list of data to a file suitable for gnuplotting
#
def write_data(data, fileName):
	out = open(fileName, "w")
	for line in data:
		for column in line:
			out.write(str(column) + "\t")
		out.write("\n")
	out.close()

def Plot(data):
	gnuplot.replot(Gnuplot.Data(data, with_="lp"))


def main():	
	
	#test = se_dist(4.0, 7.0, 0.2, 0.01)
	#write_data(test, "se_dist.dat")	
	
	
	#print("Press enter to exit")
	#sys.stdin.readline()
		
	return 0


if __name__ == "__main__":
	sys.exit(main())