diff --git a/pykat/maths/__init__.py b/pykat/math/__init__.py
similarity index 100%
rename from pykat/maths/__init__.py
rename to pykat/math/__init__.py
diff --git a/pykat/maths/hermite.py b/pykat/math/hermite.py
similarity index 100%
rename from pykat/maths/hermite.py
rename to pykat/math/hermite.py
diff --git a/pykat/maths/zernike.py b/pykat/math/zernike.py
similarity index 100%
rename from pykat/maths/zernike.py
rename to pykat/math/zernike.py
diff --git a/pykat/optics/gaussian_beams.py b/pykat/optics/gaussian_beams.py
index ce156326b8a1d7bdd09c66223aab9b866bdaba6c..43d047008c7aa2c89e57a1aa26172d9e7e967d91 100644
--- a/pykat/optics/gaussian_beams.py
+++ b/pykat/optics/gaussian_beams.py
@@ -10,7 +10,8 @@ import copy
 import warnings
 import cmath
 from math import factorial
-from scipy.special import hermite, jacobi
+from scipy.special import hermite, eval_jacobi
+from pykat.math.jacobi import jacobi
 from pykat.SIfloat import SIfloat
 
 class gauss_param(object):
@@ -271,11 +272,11 @@ class beam_param(gauss_param):
     pass
 
 # Should be renamed to HG_mode?    
-class HG_beam(object):
+class HG_mode(object):
     """ Hermite-Gauss beam profile. Example usage:
     import pykat.optics.gaussian_beams as gb
     qx=gb.beam_param(w0=1e-3,z=0)
-    beam=gb.HG_beam(qx,n=2,m=0)
+    beam=gb.HG_mode(qx,n=2,m=0)
     beam.plot()
     """    
     def __init__(self, qx, qy=None, n=0, m=0):
@@ -386,7 +387,7 @@ class HG_beam(object):
         return np.outer(_un, _um)
         
     def plot(self, ndx=100, ndy=100, xscale=4, yscale=4):
-        """ Make a simple plot the HG_beam """
+        """ Make a simple plot the HG_mode """
         import pykat.plotting 
         import matplotlib.pyplot as plt
         
@@ -405,15 +406,11 @@ class HG_beam(object):
         cbar = fig.colorbar(axes)
         plt.show()
         
-        
-
 
-        
-
-def hg2lg(n,m):
+def HG2LG(n,m):
     """A function for Matlab which returns the coefficients and mode indices of
     the LG modes required to create a particular HG mode.
-    Usage: coefficients,ps,ls = hg2lg(n,m)
+    Usage: coefficients,ps,ls = HG2LG(n,m)
     
     n,m:          Indces of the HG mode to re-create.
     coeffcients:  Complex coefficients for each order=n+m LG mode required to
@@ -444,10 +441,61 @@ def hg2lg(n,m):
 
         # Coefficient
         c = (signl*1j)**m * np.sqrt(factorial(N-m)*factorial(m)/(2**N * factorial(np.abs(l)+p)*factorial(p)))
-        c = c * (-1.0)**p * (-2.0)**m * scipy.special.eval_jacobi(m,np.abs(l)+p-m,p-m,0.0)
+        c = c * (-1.0)**p * (-2.0)**m * eval_jacobi(m,np.abs(l)+p-m,p-m,0.0)
 
         coefficients[j] = c
         
     return coefficients, ps, ls 
         
 
+    
+def LG2HG(p,l):
+    """ Function to compute the amplitude coefficients
+    of Hermite-Gauss modes whose sum yields a Laguerre Gauss mode
+    of order n,m.
+    Usage: coefficients, ns, ms = LG2HG(p,l)
+    p:    Radial LG index
+    l:    Azimuthal LG index
+    The LG mode is written as u_pl with 0<=|l|<=p.
+    The output is a series of coefficients for u_nm modes,
+    given as complex numbers and respective n,m indices
+    coefficients (complex array): field amplitude for mode u_nm
+    ns (int array): n-index of u_nm
+    ms (int array): m-index of u_nm
+
+    
+    The formula is adpated from M.W. Beijersbergen et al 'Astigmatic
+    laser mode converters and transfer of orbital angular momentum',
+    Opt. Comm. 96 123-132 (1993)
+    We adapted coefficients to be compatible with our
+    definition of an LG mode, it differs from
+    Beijersbergen by a (-1)^p factor and has exp(il\phi) rather
+    than exp(-il\phi).  Also adapted for allowing -l.
+    Andreas Freise, Charlotte Bond    25.03.2007"""
+
+    # Mode order
+    N=2*p+np.abs(l)
+
+    # Indices for coefficients
+    n = np.abs(l)+p
+    m = p
+
+    # create empty vectors
+    coefficients = np.linspace(0,0,N+1,dtype=np.complex_)
+    ns = np.linspace(0,0,N+1)
+    ms = np.linspace(0,0,N+1)
+
+    # l positive or negative
+    signl = np.sign(l)
+    if (l==0):
+        signl = 1.0
+    
+    # Beijersbergen coefficients
+    for j in np.arange(0,N+1):
+        ns[j]=N-j
+        ms[j]=j
+
+        c=(-signl*1j)**j * math.sqrt(factorial(N-j)*factorial(j)/(2**N * factorial(n)*factorial(m)))
+        coefficients[j] = c * (-1.0)**p * (-2)**j * jacobi(j,n-j,m-j,0.0)
+    
+    return coefficients, ns, ms
diff --git a/pykat/optics/knm.py b/pykat/optics/knm.py
index 9a412baa4ba4b6417e9a9a1f6b21bf55e3161ce6..c36aaf5b44979f1496990799a1b29885567a223f 100644
--- a/pykat/optics/knm.py
+++ b/pykat/optics/knm.py
@@ -1,5 +1,5 @@
 from itertools import combinations_with_replacement as combinations
-from pykat.optics.gaussian_beams import beam_param, HG_beam
+from pykat.optics.gaussian_beams import beam_param, HG_mode
 from pykat.exceptions import BasePyKatException
 from pykat.optics.romhom import u_star_u
 from pykat.external.progressbar import ProgressBar, ETA, Percentage, Bar
@@ -7,10 +7,10 @@ from scipy.interpolate import interp2d
 from scipy.integrate import dblquad
 from pykat.optics.romhom import ROMWeights
 from math import factorial
-from pykat.maths.hermite import hermite
+from pykat.math.hermite import hermite
 from scipy.misc import comb
 from scipy.integrate import newton_cotes
-from pykat.maths import newton_weights
+from pykat.math import newton_weights
 
 import time
 import pykat.optics.maps
@@ -65,8 +65,8 @@ def adaptive_knm(mode_in, mode_out, q1, q2, q1y=None, q2y=None, smap=None, delta
     if len(mode_in) != 2 or len(mode_out) != 2:
         raise BasePyKatException("Both mode in and out should be a container with modes [n m]")
     
-    Hg_in  = HG_beam(qx=q1, qy=q1y, n=mode_in[0], m=mode_in[1])
-    Hg_out = HG_beam(qx=q2, qy=q2y, n=mode_out[0], m=mode_out[1])
+    Hg_in  = HG_mode(qx=q1, qy=q1y, n=mode_in[0], m=mode_in[1])
+    Hg_out = HG_mode(qx=q2, qy=q2y, n=mode_out[0], m=mode_out[1])
     
     Nfuncs = []
     Nfuncs.append(0)
@@ -168,8 +168,8 @@ def riemann_HG_knm(x, y, mode_in, mode_out, q1, q2, q1y=None, q2y=None,
     dy = abs(y[1] - y[0])    
         
     if cache is None:
-        Hg_in  = HG_beam(qx=q1, qy=q1y, n=mode_in[0], m=mode_in[1])
-        Hg_out = HG_beam(qx=q2, qy=q2y, n=mode_out[0], m=mode_out[1])
+        Hg_in  = HG_mode(qx=q1, qy=q1y, n=mode_in[0], m=mode_in[1])
+        Hg_out = HG_mode(qx=q2, qy=q2y, n=mode_out[0], m=mode_out[1])
         
         U1 = Hg_in.Unm(x+delta[0], y+delta[1])
         U2 = Hg_out.Unm(x,y).conjugate()
@@ -505,8 +505,8 @@ def square_aperture_HG_knm(mode_in, mode_out, q, R):
     m = mode_in[1]
     _m = mode_out[1]
     
-    hg1 = HG_beam(q, n=n, m=m)
-    hg2 = HG_beam(q, n=_n, m=_m)
+    hg1 = HG_mode(q, n=n, m=m)
+    hg2 = HG_mode(q, n=_n, m=_m)
         
     kx = hg1.constant_x * hg2.constant_x.conjugate()
     ky = hg1.constant_y * hg2.constant_y.conjugate()
diff --git a/pykat/optics/romhom.py b/pykat/optics/romhom.py
index 238616e6748aa49add0028f8be58be0ae3dc6ef6..392523581af8502b4999b75e23ab326d86efb004 100644
--- a/pykat/optics/romhom.py
+++ b/pykat/optics/romhom.py
@@ -13,11 +13,11 @@ import itertools
 from copy import copy
 from pykat.external.progressbar import ProgressBar, ETA, Percentage, Bar
 from itertools import combinations_with_replacement as combinations
-from pykat.optics.gaussian_beams import beam_param, HG_beam
+from pykat.optics.gaussian_beams import beam_param
 from scipy.linalg import inv
 from math import factorial
-from pykat.maths.hermite import *
-from pykat.maths import newton_weights
+from pykat.math.hermite import *
+from pykat.math import newton_weights
 from scipy.integrate import newton_cotes
 from multiprocessing import Process, Queue, Array, Value, Event
 from pykat.exceptions import BasePyKatException
@@ -658,4 +658,4 @@ def makeWeightsNew(smap, EIxFilename, EIyFilename=None, verbose=True, newtonCote
     if verbose:
         p.finish()
     
-    return ROMWeights(w_ij_Q1=w_ij_Q1, w_ij_Q2=w_ij_Q2, w_ij_Q3=w_ij_Q3, w_ij_Q4=w_ij_Q4, EIx=EIx, EIy=EIy, direction=direction)
\ No newline at end of file
+    return ROMWeights(w_ij_Q1=w_ij_Q1, w_ij_Q2=w_ij_Q2, w_ij_Q3=w_ij_Q3, w_ij_Q4=w_ij_Q4, EIx=EIx, EIy=EIy, direction=direction)