diff --git a/bin/test_aperture.py b/bin/test_aperture.py
index 65efee1613c6246ee60ebab5b68de16d77bb520d..03bcee70685592c1f7527009bd1ef6e46f779ee1 100644
--- a/bin/test_aperture.py
+++ b/bin/test_aperture.py
@@ -1,3 +1,4 @@
+from pykat.utilities.optics.gaussian_beams import gauss_param
 from pykat import finesse
 from pykat.commands import xaxis
 import pylab as pl
@@ -19,7 +20,7 @@ kat.parseCommands(code)
 
 maxtem = np.arange(0, 2, 2)
 
-kat.m1.n2.q = 1j*(math.pi * 1e-3)**2/1064e-9
+kat.m1.n2.q = gauss_param(w0=1e-3, z=0)
 
 for tem in maxtem:
     print "Calculating maxtem ", tem, "..."
diff --git a/pykat/commands.py b/pykat/commands.py
index 45ab2ca9a5f6307e5d04745170266b44f0866406..483fa0a67560b0626258ab487081221457b8b2e4 100644
--- a/pykat/commands.py
+++ b/pykat/commands.py
@@ -38,8 +38,17 @@ class cavity(Command):
         
     def getFinesseText(self):
         return 'cav {0} {1} {2} {3} {4}'.format(self.__name, self.__c1, self.__n1, self.__c2, self.__n2);
-    
 
+class gauss(object):
+    @staticmethod
+    def parseFinesseText(text, kat):  
+        values = text.split(" ")
+        
+        if not values[0].startswith("gauss"):
+            raise exceptions.RuntimeError("'{0}' not a valid Finesse gauss command".format(text))
+        
+        
+        
 class attr():
     @staticmethod
     def parseFinesseText(text, kat):  
diff --git a/pykat/finesse.py b/pykat/finesse.py
index 0582e5f97fce9eb62cda34847b24f5d8f4b26def..c55fbfb9b4774112d957073eaf84da2bf1f84945 100644
--- a/pykat/finesse.py
+++ b/pykat/finesse.py
@@ -149,6 +149,9 @@ class kat(object):
         
         commands=self.remove_comments(commands)
         
+        after_process = [] # list of commands that should be processed after 
+                           # objects have been set and created
+        
         for line in commands.split("\n"):
             #for line in commands:
             if len(line.strip()) >= 2:
@@ -200,15 +203,25 @@ class kat(object):
                     obj = pykat.detectors.photodiode.parseFinesseText(line)
                 elif(first == "xaxis" or first == "x2axis" or first == "xaxis*" or first == "x2axis*"):
                     obj = pykat.commands.xaxis.parseFinesseText(line)
+                elif(first == "gauss" or first == "gauss*" or first == "gauss**"):
+                    after_process.append(line)
                 else:
                     print "Parsing `{0}` into pykat object not implemented yet, added as extra line.".format(line)
                     obj = line
                     # manually add the line to the block contents
                     self.__blocks[self.__currentTag].contents.append(line) 
                 
-                if not isinstance(obj, str):
+                if obj != None and not isinstance(obj, str):
                     self.add(obj)
                     
+        # now process all the varous gauss/attr etc. commands which require
+        # components to exist first before they can be processed
+        for line in after_process:
+            first = line.split(" ",1)[0]
+            
+            if first == "gauss" or first == "gauss*" or first == "gauss**":
+                pykat.commands.gauss.parseFinesseText(line)
+            
         self.__currentTag = NO_BLOCK 
             
     def run(self, printout=1, printerr=1, save_output=False, save_kat=False,kat_name=None) :
diff --git a/pykat/utilities/__init__.py b/pykat/utilities/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/pykat/utilities/optics/ABCD.py b/pykat/utilities/optics/ABCD.py
new file mode 100644
index 0000000000000000000000000000000000000000..58ecc2e61951a564ec9f25a4d6f8a905d09f7203
--- /dev/null
+++ b/pykat/utilities/optics/ABCD.py
@@ -0,0 +1,16 @@
+import numpy as np
+
+def apply(ABCD, q1, n1, n2):
+    return n2 * (ABCD[0,0] * q1/n1 + ABCD[0,1] / (ABCD[1,0] * q1/n1 + ABCD[1,1]))
+
+def mirror_trans(n1, n2, Rc):
+    return np.matrix([[1,0],[(n2-n1)/Rc,1]])
+    
+def mirror_refl(n, Rc):
+    return np.matrix([[1,0],[-2*n/Rc,1]])
+    
+def space(n,L):
+    return np.matrix([[1, L/n],[0,1]])
+    
+def lens(f):
+    return np.matrix([[1, 0],[-1/f,1]])
\ No newline at end of file
diff --git a/pykat/utilities/optics/__init__.py b/pykat/utilities/optics/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/pykat/utilities/optics/gaussian_beams.py b/pykat/utilities/optics/gaussian_beams.py
new file mode 100644
index 0000000000000000000000000000000000000000..568dbd5f29c2ac18faaafc371ae22acad4940104
--- /dev/null
+++ b/pykat/utilities/optics/gaussian_beams.py
@@ -0,0 +1,126 @@
+import pykat.exceptions as pkex
+import numpy
+import math
+
+class gauss_param(object):
+    """
+    Gaussian beam complex parameter
+    
+    gauss_param is effectively a complex number with extra
+    functionality to determine beam parameters.
+    
+    defaults to 1064e-9m for wavelength and refractive index 1
+    usage:
+        q = gauss_param(w0=w0, z=z)
+        q = gauss_param(z=z, zr=zr)
+        q = gauss_param(wz=wz, rc=rc)
+        
+        or change default wavelength and refractive index with:
+        
+        q = gauss_param(wavelength, nr, w0=w0, zr=zr)
+    """
+    
+    def __init__(self, wavelength=1064e-9, nr=1, *args, **kwargs):
+        self.__q = None
+        self.__lambda = wavelength
+        self.__nr = nr
+        
+        if len(args) == 1:
+            self.__q = args[0]
+        elif len(kwargs) == 2:        
+            
+            if "w0" in kwargs and "z" in kwargs:
+                q = float(kwargs["z"]) + 1j *float(math.pi*kwargs["w0"]**2/(self.__lambda/self.__nr) )
+            elif "z" in kwargs and "zr" in kwargs:
+                q = float(kwargs["z"]) + 1j *float(kwargs["zr"]) 
+            elif "rc" in kwargs and "wz" in kwargs:
+                one_q = 1 / kwargs["rc"] - 1j * self.__lamda / (math.pi * self.__nr * kwargs["wz"]**2)
+                q = 1/one_q
+            else:
+                raise pkex.BasePyKatException("Must specify: z and w0 or z and zr or rc and wz, to define the beam parameter")
+                
+            self.__q = q
+        else:
+            raise pkex.BasePyKatException("Incorrect usage for gauss_param constructor")
+    
+    @property
+    def wavelength(self): return self.__lambda
+    
+    @property
+    def nr(self): return self.__nr
+    
+    @property
+    def q(self): return self.__q
+    
+    @property
+    def z(self): return self.__q.real
+    
+    @property
+    def zr(self): return self.__q.imag
+    
+    @property
+    def wz(self):
+        return math.sqrt(self.__lambda /(self.__nr * math.pi) * abs(self.__q) / self.__q.imag)
+    
+    @property
+    def w0(self):
+        return math.sqrt(self.__q.imag * self.__lambda / (self.__nr * math.pi))    
+
+    @property
+    def Rc(self):
+        if self.__q.real != 0:
+            return abs(self.__q) / self.__q.real
+        else:
+            return float("inf")
+    
+    def conjugate(self):
+        return gauss_param(self.__lambda, self.__nr, self.__q.conjugate())
+    
+    def __complex__(self):
+        return self.__q
+    
+    def __str__(self):
+        return str(self.__q)
+    
+    def __mul__(self, a):
+        return gauss_param(self.__lambda, self.__nr, self.__q * complex(a))
+    
+    def __imul__(self, a):
+        self.__q += complex(a)
+        return self
+        
+    __rmul__ = __mul__
+    
+    def __add__(self, a):
+        return gauss_param(self.__lambda, self.__nr, self.__q + complex(a))
+    
+    def __iadd__(self, a):
+        self.__q += complex(a)
+        return self
+        
+    __radd__ = __add__
+    
+    def __sub__(self, a):
+        return gauss_param(self.__lambda, self.__nr, self.__q - complex(a))
+    
+    def __isub__(self, a):
+        self.__q -= complex(a)
+        return self
+        
+    __rsub__ = __sub__
+    
+    def __div__(self, a):
+        return gauss_param(self.__lambda, self.__nr, self.__q / complex(a))
+    
+    def __idiv__(self, a):
+        self.__q /= complex(a)
+        return self
+    
+    def __pow__(self, q):
+        return gauss_param(self.__lambda, self.__nr, self.__q**q)
+
+    def __neg__(self, q):
+        return gauss_param(self.__lambda, self.__nr, -self.__q)
+        
+    def __eq__(self, q):
+        return complex(q) == self.__q
\ No newline at end of file