Commit before breaking everything

[matches/honours.git] / research / analysis and stuff / theory.py

diff --git a/research/analysis and stuff/theory.py b/research/analysis and stuff/theory.py
new file mode 100755 (executable)
index 0000000..8367f1d
--- /dev/null
+++ b/research/analysis and stuff/theory.py
@@ -0,0 +1,111 @@
+#!/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.plot(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())