thesis/theory.py

   1 #!/usr/bin/python -u
   2
   3 #
   4 # @file theory.py
   5 # @purpose Do theory stuff
   6 #               Because sadly, Mathematica does not work on my new laptop
   7 # @author Sam Moore
   8 # @date August 2012
   9 #
  10
  11 import sys
  12 import os
  13 import math
  14 import re # Regular expressions - for removing comments
  15 #import odict #ordered dictionary
  16
  17 import Gnuplot, Gnuplot.funcutils
  18
  19 gnuplot = Gnuplot.Gnuplot()
  20
  21
  22 #
  23 # @function se_dist
  24 # @purpose Model the secondary electron distribution using M. Furman's formula
  25 #
  26 def se_dist(s, Ep, sigma_p, dE, n_max=0.75):
  27         results = []
  28         E = 0.0
  29         while (E < 8.0):
  30                 results.append([E, n_max *s * E / (s - 1.0 + E**s) , math.exp(-(E - Ep)**2.0 / (2.0*sigma_p**2.0))])
  31                 E += dE
  32         return results
  33
  34 #
  35 # @function integrate
  36 # @purpose integrate function numerically over a range of arguments
  37 #
  38 def integrate(f, xmin, xmax, dx):
  39         x = xmin
  40         result = 0.00
  41         while (x <= xmax):
  42                 result += f(x) * dx
  43                 x += dx
  44         return result
  45
  46 # @function derivative
  47 # @purpose get derivative of a function at a value
  48 def der(f, x, dx):
  49         return (f(x + dx) - f(x)) / dx
  50
  51 #
  52 # @function tcs
  53 # @purpose Model S(E)
  54 def tcs(f, sigma, Emin, Emax, dE):
  55         results = []
  56         E = Emin
  57         while (E < Emax):
  58                 results.append([E, (1.0 - sigma(0)) * f(-E) - integrate(lambda e : f(e - E) * der(sigma, E, dE), Emin, Emax, dE)])
  59                 E += dE
  60         return results
  61
  62 def delta(x):
  63         if (x == 0):
  64                 return 1.0
  65         else:
  66                 return 0.0
  67
  68 def table(f, xmin, xmax, dx):
  69         result = []
  70         x = xmin
  71         while (x <= xmax):
  72                 result.append([x, f(x)])
  73                 x += dx
  74         return result
  75
  76 def gaussian(x, sigma):
  77         return math.exp(- (x**2.0)/(2.0 * sigma**2.0)) / (sigma * (2.0 * math.pi)**0.50)
  78
  79 def step(x, sigma, T):
  80         return 1.0 / (math.exp((x - sigma)/T) + 1.0)
  81
  82 #
  83 # @function write_data
  84 # @purpose Write a list of data to a file suitable for gnuplotting
  85 #
  86 def write_data(data, fileName):
  87         out = open(fileName, "w")
  88         for line in data:
  89                 for column in line:
  90                         out.write(str(column) + "\t")
  91                 out.write("\n")
  92         out.close()
  93
  94 def Plot(data):
  95         gnuplot.replot(Gnuplot.Data(data, with_="lp"))
  96
  97
  98 def main():
  99
 100         #test = se_dist(4.0, 7.0, 0.2, 0.01)
 101         #write_data(test, "se_dist.dat")
 102
 103
 104         #print("Press enter to exit")
 105         #sys.stdin.readline()
 106
 107         return 0
 108
 109
 110 if __name__ == "__main__":
 111         sys.exit(main())