Merge branch 'master' of gitmaster.atlas.aei.uni-hannover.de:pykat/pykat

bin/test_2d_plot.py

+from pykat import finesse
+from pykat.utilities.optics.gaussian_beams import gauss_param
+from pykat.detectors import *
+from pykat.components import *
+from pykat.commands import *
+from pykat.structs import *
+import matplotlib.image as mpimg
+from mpl_toolkits.mplot3d.axes3d import Axes3D
+
+import numpy as np
+import pylab as pl
+
+code = """
+l l1 1 0 0 n1
+s s1 10 1 n1 n2
+
+ad ad1 0 n2
+beam b1 0 n2
+maxtem 0 
+
+xaxis b1 x lin -10.0 10 100
+x2axis b1 y lin -6 6 100
+
+"""
+
+kat = finesse.kat()
+
+kat.parseCommands(code)
+kat.s1.n1.q = gauss_param(w0=1e-3, z=0)
+
+out = kat.run(printout=0,printerr=0)
+
+fig = pl.figure()
+ax = fig.add_subplot(1, 1, 1, projection='3d')
+x, y = np.meshgrid(out.x, out.y)
+p = ax.plot_wireframe(x, y, out["b1"])   
+pl.xlabel(out.xlabel)
+pl.ylabel(out.ylabel)
+pl.show()
+
+pl.figure()
+extent = [np.min(out.x),np.max(out.x),np.min(out.y),np.max(out.y)]
+imgplot = pl.imshow(out["b1"], extent=extent)
+#imgplot.set_interpolation('bicubic')
+imgplot.set_interpolation('nearest')
+pl.colorbar()
+pl.xlabel(out.xlabel)
+pl.ylabel(out.ylabel)
+pl.show()

bin/test_plot.py

@@ -39,7 +39,7 @@ kat.maxtem = 0
 out = kat.run(printout=0,printerr=0)
 
 pl.figure()
-pl.plot(out.x,out.y)
+pl.plot(out.x, out["pd_cav"])
 pl.xlabel(out.xlabel)
 pl.ylabel("Intensity [W]")
 pl.legend(out.ylabels)

pykat/commands.py

@@ -77,7 +77,8 @@ class attr():
        
 class xaxis(Command):
     
-    def __init__(self, scale, limits, comp, param, steps):
+    def __init__(self, scale, limits, comp, param, steps, axis_type="xaxis"):
+        self._axis_type = axis_type
         
         if scale == "lin":
             scale = Scale.linear
@@ -121,23 +122,43 @@ class xaxis(Command):
     def parseFinesseText(text):
         values = text.split(" ")
         
-        if values[0] != "xaxis" and values[0] != "xaxis*" and values[0] != "x2axis" and values[0] != "x2axis*":
+        if values[0] != "xaxis" and values[0] != "xaxis*":
             raise exceptions.RuntimeError("'{0}' not a valid Finesse xaxis command".format(text))
         
+        axis_type = values[0]
+        
         values.pop(0) # remove initial value
         
         if len(values) != 6:
             raise exceptions.RuntimeError("xaxis Finesse code format incorrect '{0}'".format(text))
 
-        return xaxis(values[2], [values[3], values[4]], values[0], values[1], values[5])
+        return xaxis(values[2], [values[3], values[4]], values[0], values[1], values[5], axis_type=axis_type)
         
     def getFinesseText(self):
         # store either the component name of the string provided
         comp_name = self.__comp.name if isinstance(self.__comp, Component) else self.__comp
         param_name = self.__param.name if isinstance(self.__param, Param) else self.__param
         
-        return 'xaxis {0} {1} {2} {3} {4} {5}'.format(
+        return '{axis_type} {0} {1} {2} {3} {4} {5}'.format(
                 comp_name, param_name, self.scale,
-                min(self.limits), max(self.limits), self.steps);
+                min(self.limits), max(self.limits), self.steps, axis_type=self._axis_type);
+                
+class x2axis(xaxis):
+    def __init__(self, scale, limits, comp, param, steps):
+        xaxis.__init__(self, scale, limits, comp, param, steps, axis_type="x2axis")        
+
+    @staticmethod
+    def parseFinesseText(text):
+        values = text.split(" ")
         
+        if values[0] != "x2axis" and values[0] != "x2axis*":
+            raise exceptions.RuntimeError("'{0}' not a valid Finesse xaxis command".format(text))
+        
+        axis_type = values[0]
+        
+        values.pop(0) # remove initial value
+        
+        if len(values) != 6:
+            raise exceptions.RuntimeError("xaxis Finesse code format incorrect '{0}'".format(text))
 
+        return x2axis(values[2], [values[3], values[4]], values[0], values[1], values[5])
\ No newline at end of file

pykat/components.py

@@ -143,6 +143,8 @@ class Component(object):
     @property
     def id(self): return self.__id
     
+    def __str__(self): return self.name
+    
 class Param(float):
     def __new__(self,name,value):
         return float.__new__(self,SIfloat(value))
@@ -345,7 +347,7 @@ class space(Component):
         
         self._requested_node_names.append(node1)
         self._requested_node_names.append(node2)
-        
+        self._QItem = None
         self.__L = SIfloat(L)
         self.__n = SIfloat(n)
         

pykat/detectors.py

@@ -44,10 +44,10 @@ class Detector(object) :
     @property 
     def node(self): return self.__node
     
-    def __getname(self):
-        return self.__name        
+    @property
+    def name(self): return self.__name        
 
-    name = property(__getname)
+    def __str__(self): return self.name
     
 class photodiode(Detector):
 

pykat/finesse.py

@@ -34,7 +34,7 @@ import pykat
 import warnings
 import re
 
-from pykat.exceptions import *
+import pykat.exceptions as pkex
 
 from pykat.node_network import NodeNetwork
 from pykat.detectors import Detector
@@ -67,6 +67,48 @@ class katRun(object):
         with open(filename,'r') as infile:
             return pickle.load(infile)
     
+    def get(self, value): return self[value]
+    
+    def __getitem__(self, value):
+        if value in self.ylabels:
+            idx = self.ylabels.index(value)
+            if len(self.y.shape) == 1:
+                return self.y
+            else:
+                return self.y[:, idx]
+        else:
+            raise  pkex.BasePyKatException("No output by the name {0} found", value)
+      
+class katRun2D(object):
+    def __init__(self):
+        self.runDateTime = datetime.datetime.now()
+        self.x = None
+        self.y = None
+        self.z = None
+        self.xlabel = None
+        self.ylabel = None
+        self.zlabels = None
+        self.katScript = None
+        self.katVersion = None
+        
+    def saveKatRun(self, filename):
+        with open(filename,'w') as outfile:
+            pickle.dump(self, outfile)
+    
+    @staticmethod
+    def loadKatRun(filename):
+        with open(filename,'r') as infile:
+            return pickle.load(infile)
+    
+    def get(self, value): return self[value]
+    
+    def __getitem__(self, value):
+        idx = [i for i in range(len(self.zlabels)) if self.zlabels[i].split(" ")[0] == str(value)]
+        
+        if len(idx) > 0:
+            return self.z[idx].squeeze()
+        else:
+            raise  pkex.BasePyKatException("No output by the name {0} found".format(str(value)))
       
 class Block:
     def __init__(self, name):
@@ -225,8 +267,10 @@ class kat(object):
                     obj = pykat.components.lens.parseFinesseText(line)
                 elif(first[0:2] == "pd"):
                     obj = pykat.detectors.photodiode.parseFinesseText(line)
-                elif(first == "xaxis" or first == "x2axis" or first == "xaxis*" or first == "x2axis*"):
+                elif(first == "xaxis" or first == "xaxis*"):
                     obj = pykat.commands.xaxis.parseFinesseText(line)
+                elif(first == "x2axis" or first == "x2axis*"):
+                    obj = pykat.commands.x2axis.parseFinesseText(line)
                 elif(first == "gauss" or first == "gauss*" or first == "gauss**"):
                     after_process.append(line)
                 else:
@@ -263,15 +307,16 @@ class kat(object):
         except BasePyKatException as ex:
             print ex
             
-    def run(self, printout=1, printerr=1, save_output=False, save_kat=False,kat_name=None) :
+    def run(self, printout=0, printerr=0, save_output=False, save_kat=False,kat_name=None) :
         """ 
         Runs the current simulation setup that has been built thus far.
         It returns a katRun object which is populated with the various
         data from the simulation run.
         """
+        
         try:        
-            r = katRun()
-            r.katScript = "".join(self.generateKatScript())       
+            if not hasattr(self, "xaxis"):
+                raise pkex.BasePyKatException("No xaxis was defined")
             
             if len(self.__katdir) == 0:
                 # Get the environment variable for where Finesse is stored
@@ -296,6 +341,17 @@ class kat(object):
             if not (os.path.isfile(kat_exec) and os.access(kat_exec, os.X_OK)):
                 raise MissingFinesse()
                 
+            print "--------------------------------------------------------------"
+            start = datetime.datetime.now()
+            print "Running kat - Started at " + str(start)
+            
+            if hasattr(self, "x2axis"):
+                r = katRun2D()
+            else:
+                r = katRun()
+                
+            r.katScript = "".join(self.generateKatScript())   
+            
             # create a kat file which we will write the script into
             if self.__tempname == None:
                 katfile = tempfile.NamedTemporaryFile(suffix=".kat", dir=self.__tempdir)
@@ -307,27 +363,36 @@ class kat(object):
             katfile.flush()
             
             cmd=[kat_exec, '--perl1']
+            
             if self.__time_code:
                 cmd.append('--perf-timing')
+
             cmd.append('--no-backspace')
+            # set default format so that less repeated numbers are printed to the
+            # output file, should speed up running and parsing of output files
+            cmd.append('-format=%.15g')
 
             cmd.append(katfile.name)
-            if self.verbose:
-                print cmd
+            
+            #if self.verbose:
+                #print cmd
+                
             p=subprocess.Popen(cmd, shell=False, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
             err = ""
             
+            print "Finesse binary output:"
+            
             for line in iter(p.stderr.readline, ""):
-                err += line
-                vals = line.split("-")
+                #err += line 
                 
-                if len(vals) == 2:
+                if len(line) > 0:
+                    if line.rstrip().endswith('%'):
+                        vals = line.split("-")
                         action = vals[0].strip()
                         prc = vals[1].strip()[:-1]
-                    
-                    #sys.stdout.write("\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b")
                         sys.stdout.write("\r{0} {1}%".format(action, prc))
-                    sys.stdout.flush()
+                    elif line[0] == '*':
+                        sys.stdout.write(line)
                         
             [out,errpipe] = p.communicate()
             
@@ -348,11 +413,6 @@ class kat(object):
             base = os.path.basename(root[0])            
             outfile = root[0] + ".out"
             
-            [r.x,r.y,hdr] = self.readOutFile(outfile)
-            
-            r.xlabel = hdr[0]
-            r.ylabels = hdr[1:]
-            
             if save_output:        
                 newoutfile = "{0}.out".format(base)
                 
@@ -364,7 +424,19 @@ class kat(object):
                     
                 os.rename(outfile, newoutfile)
 
-                print "Output data saved to '{0}'".format(newoutfile)
+                print "\nOutput data saved to '{0}'".format(newoutfile)
+            
+            if hasattr(self, "x2axis"):
+                [r.x,r.y,r.z,hdr] = self.readOutFile(outfile)
+                
+                r.xlabel = hdr[0]
+                r.ylabel = hdr[1]
+                r.zlabels = map(str.strip, hdr[2:])
+            else:
+                [r.x,r.y,hdr] = self.readOutFile(outfile)
+            
+                r.xlabel = hdr[0]
+                r.ylabels = map(str.strip, hdr[1:])
                             
             if save_kat:
                 if kat_name == None:
@@ -397,6 +469,9 @@ class kat(object):
             
         except FinesseRunError as fe:
             print fe
+        finally:
+            print ""
+            print "Finished in " + str(datetime.datetime.now()-start)
             
         
     def add(self, obj):
@@ -435,15 +510,27 @@ class kat(object):
 
     def readOutFile(self, filename):
         
-        outfile = open(filename,'r')
-        
+        with open(filename,'r') as outfile:
             # read first to lines to get to header line
             outfile.readline()
             outfile.readline()
             
             hdr = outfile.readline().replace('%','').replace('\n','').split(',')
         
-        data = np.loadtxt(filename,comments='%')
+        data = np.loadtxt(filename,comments='%',skiprows=4)
+        
+        if hasattr(self, "x2axis"):
+            # need to parse 2D outputs slightly different as they are effectively 2D matrices
+            # written in linear form
+            x = data[0::(1+self.x2axis.steps),0]
+            y = data[0:(1+self.x2axis.steps),1]
+            # get rows and columns lined up so that we can reshape a single column of all x/y data
+            # into a matrix
+            z = data[:,2:].transpose().reshape(data.shape[1]-2, 1+self.xaxis.steps, 1+self.x2axis.steps)
+            # once you do this the data for y and x axes need swapping
+            z = z.swapaxes(1,2)
+            return [x, y, z, hdr]
+        else:
             shape_len = len(data.shape)
             
             if shape_len > 1:

pykat/utilities/optics/gaussian_beams.py

 import pykat.exceptions as pkex
 import numpy
 import math
+import copy
 
 class gauss_param(object):
     """
@@ -142,3 +143,48 @@ class gauss_param(object):
     def imag(self): return self.__q.imag
     @imag.setter
     def imag(self, value): self.__q.imag = SIfloat(value)
+    
+class HG_gauss_beam(object):
+    
+    def __init__(self, qx, qy=None, n=0, m=0):
+        self._qx = copy.deepcopy(qx)
+        self._2pi_qrt = math.pow(2.0/math.pi, 0.25)
+        
+        if qy == None:
+            self._q0 = copy.deepcopy(qx)
+        else:
+            self._q0 = copy.deepcopy(qy)
+    
+        self.n = n
+        self.m = m
+        self._calc_constants()
+        
+    @property
+    def n(self): return self._n
+    @n.setter
+    def n(self,value): 
+        self._n = float(value)
+        self._calc_constants()
+
+    @property
+    def m(self): return self._m
+    @m.setter
+    def m(self,value): 
+        self._m = float(value)
+        self._calc_constants()
+        
+    def _calc_constants(self):
+        self.__xpre_const = math.pow(2.0/math.pi, 0.25)
+        self.__xpre_const *= math.sqrt(1/(2**self._n * math.factorial(self._n)))
+        self.__xpre_const *= math.sqrt(1j*self._qx.imag / self._qx)
+        self.__xpre_const *= math.pow((1j*self._qx.imag * self._qx.conjugate)/(-1j*self._qx.imag * self._qx), self._n/2.0)
+        
+        self.__ypre_const = math.pow(2.0/math.pi, 0.25)
+        self.__ypre_const *= math.sqrt(1/(2**self._m * math.factorial(self._m)))
+        self.__ypre_const *= math.sqrt(1j*self._qy.imag / self._qy)
+        self.__ypre_const *= math.pow((1j*self._qy.imag * self._qy.conjugate)/(-1j*self._qy.imag * self._qy), self._m/2.0)
+                
+    def Unm(self, n, m, x, y):    
+        return self.__xpre_const * special
+        
+        
\ No newline at end of file