diff --git a/bin/test_hg_beam.py b/bin/test_hg_beam.py
new file mode 100644
index 0000000000000000000000000000000000000000..370407fdbaff41b0cb74ba768068d722b25c43fa
--- /dev/null
+++ b/bin/test_hg_beam.py
@@ -0,0 +1,11 @@
+import pykat
+from pykat.utilities.optics.gaussian_beams import *
+
+gx = beam_param(w0=2e-3, z=0)
+gy = beam_param(w0=1e-3, z=0)
+
+beam = HG_beam(gx,gy,0,0)
+beam.n = 5
+beam.m = 6
+
+beam.plot()
\ No newline at end of file
diff --git a/bin/test_remove.py b/bin/test_remove.py
new file mode 100644
index 0000000000000000000000000000000000000000..ce5222441618400c4584edcd31f4b6952529e45a
--- /dev/null
+++ b/bin/test_remove.py
@@ -0,0 +1,40 @@
+from pykat import finesse
+from pykat.detectors import *
+from pykat.components import *
+from pykat.commands import *
+from pykat.structs import *
+
+import numpy as np
+import pylab as pl
+
+code = """
+l l1 1 0 0 n1 ### test
+s s1 10 1 n1 n2
+m m1 0.5 0.5 0 n2 n3
+s s2 10 1 n3 n4
+m m2 0.5 0.5 0 n4 n5
+s s3 10 1 n5 n6
+
+yaxis abs:deg
+
+pd pd_cav n3
+
+cav c1 m1 n3 m2 n4
+
+attr m1 Rc 1
+"""
+
+kat = finesse.kat()
+
+kat.parseCommands(code)
+
+kat.add(xaxis("lin", [0, 360], kat.m2.phi, 100))
+
+kat.m1.Rcx = -1000.0
+kat.m1.Rcy = -1000.0
+kat.m2.Rcx = 1000.0
+kat.m2.Rcy = 1000.0
+
+kat.maxtem = 0
+
+kat.parseCommands("xaxis m1 phi lin -100 100 10")
\ No newline at end of file
diff --git a/pykat/commands.py b/pykat/commands.py
index d8f80cd5e51068108594c61bcbf4cf56ce2cfcfa..dfbe6f7e922319582c5c046119c139ab42735ce3 100644
--- a/pykat/commands.py
+++ b/pykat/commands.py
@@ -12,7 +12,7 @@ from structs import *
from pykat.param import Param, putter
import pykat.exceptions as pkex
from collections import namedtuple
-from pykat.utilities.optics.gaussian_beams import gauss_param
+from pykat.utilities.optics.gaussian_beams import beam_param
class Command(object):
__metaclass__ = abc.ABCMeta
@@ -73,22 +73,22 @@ class gauss(object):
if not values[0].endswith("*"):
if len(values) == 6:
- gp = gauss_param(kat.lambda0, w0=values[-2], z=values[-1])
+ gp = beam_param(kat.lambda0, w0=values[-2], z=values[-1])
elif len(values) == 8:
- gpx = gauss_param(kat.lambda0, w0=values[-4], z=values[-3])
- gpy = gauss_param(kat.lambda0, w0=values[-2], z=values[-1])
+ gpx = beam_param(kat.lambda0, w0=values[-4], z=values[-3])
+ gpy = beam_param(kat.lambda0, w0=values[-2], z=values[-1])
elif values[0].endswith("*"):
if len(values) == 6:
- gp = gauss_param(kat.lambda0, z=values[-2], zr=values[-1])
+ gp = beam_param(kat.lambda0, z=values[-2], zr=values[-1])
elif len(values) == 8:
- gpx = gauss_param(kat.lambda0, z=values[-4], zr=values[-3])
- gpy = gauss_param(kat.lambda0, z=values[-2], zr=values[-1])
+ gpx = beam_param(kat.lambda0, z=values[-4], zr=values[-3])
+ gpy = beam_param(kat.lambda0, z=values[-2], zr=values[-1])
elif values[0].endswith("**"):
if len(values) == 6:
- gp = gauss_param(kat.lambda0, w=values[-2], rc=values[-1])
+ gp = beam_param(kat.lambda0, w=values[-2], rc=values[-1])
elif len(values) == 8:
- gpx = gauss_param(kat.lambda0, w=values[-4], rc=values[-3])
- gpy = gauss_param(kat.lambda0, w=values[-2], rc=values[-1])
+ gpx = beam_param(kat.lambda0, w=values[-4], rc=values[-3])
+ gpy = beam_param(kat.lambda0, w=values[-2], rc=values[-1])
else:
raise pkex.BasePyKatException("Unexpected ending to gauss command '{0}'".format(text))
diff --git a/pykat/node_network.py b/pykat/node_network.py
index a0047929ee2d8e5ec2448855f5493ac3191ddb39..f2684c5a7471aa0b0f0314f2ec7624a483cc6592 100644
--- a/pykat/node_network.py
+++ b/pykat/node_network.py
@@ -9,7 +9,7 @@ import pykat.gui.graphics
import pykat.exceptions as pkex
from pykat.components import Component
from pykat.detectors import Detector
-from pykat.utilities.optics.gaussian_beams import gauss_param
+from pykat.utilities.optics.gaussian_beams import beam_param
class NodeNetwork(object):
def __init__(self, kat):
@@ -276,11 +276,11 @@ class Node(object):
self.__q_comp = component
if len(args) == 1:
- self.__q_x = gauss_param(self._network.kat.lambda0, q=args[0])
- self.__q_y = gauss_param(self._network.kat.lambda0, q=args[0])
+ self.__q_x = beam_param(self._network.kat.lambda0, q=args[0])
+ self.__q_y = beam_param(self._network.kat.lambda0, q=args[0])
elif len(args) == 2:
- self.__q_x = gauss_param(self._network.kat.lambda0, q=args[0])
- self.__q_y = gauss_param(self._network.kat.lambda0, q=args[1])
+ self.__q_x = beam_param(self._network.kat.lambda0, q=args[0])
+ self.__q_y = beam_param(self._network.kat.lambda0, q=args[1])
else:
raise pkex.BasePyKatException("Must specify either 1 Gaussian beam parameter or 2 for astigmatic beams")
diff --git a/pykat/utilities/optics/gaussian_beams.py b/pykat/utilities/optics/gaussian_beams.py
index 20e8fcb834755261eadc71ff2ee2cce5653151fe..f7623842498909badea29bf32cbb70118e34cf8b 100644
--- a/pykat/utilities/optics/gaussian_beams.py
+++ b/pykat/utilities/optics/gaussian_beams.py
@@ -1,11 +1,16 @@
import pykat.exceptions as pkex
-import numpy
+import numpy as np
import math
import copy
+import warnings
+import cmath
+from scipy.special import hermite
from pykat.SIfloat import SIfloat
class gauss_param(object):
"""
+ Use beam_param instead, prefer that as a name
+
Gaussian beam complex parameter
gauss_param is effectively a complex number with extra
@@ -24,8 +29,10 @@ class gauss_param(object):
"""
def __init__(self, wavelength=1064e-9, nr=1, *args, **kwargs):
+ if self.__class__ != beam_param:
+ warnings.warn("Name changed. Use beam_param instead of gauss_param.")
+
self.__q = None
- print "wavelength", wavelength, kwargs
self.__lambda = SIfloat(wavelength)
self.__nr = SIfloat(nr)
@@ -153,49 +160,141 @@ class gauss_param(object):
# reverse beam direction
def reverse(self):
self.__q = -1.0 * self.__q.real + 1j * self.__q.imag
+
+
+class beam_param(gauss_param):
+ pass
+
-class HG_gauss_beam(object):
+class HG_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)
+ if qy.__class__ == beam_param:
+ self._qy = copy.deepcopy(qx)
else:
- self._q0 = copy.deepcopy(qy)
+ self._qy = copy.deepcopy(qy)
- self.n = n
- self.m = m
+ self._n = n
+ self._m = m
+ self._hn = hermite(n)
+ self._hm = hermite(m)
self._calc_constants()
@property
def n(self): return self._n
@n.setter
def n(self,value):
- self._n = float(value)
+ self._n = int(value)
self._calc_constants()
+ self._hn = hermite(self._n)
@property
def m(self): return self._m
@m.setter
def m(self,value):
- self._m = float(value)
+ self._m = int(value)
self._calc_constants()
+ self._hm = hermite(self._m)
+
+ @property
+ def q(self):
+ if self._qx.q == self._qy.q:
+ return self._qx.q
+ else:
+ return (self._qx.q, self._qy.q)
+ @q.setter
+ def q(self, value):
+ if value.__class__ == beam_param:
+ self._qx = copy.deepcopy(value)
+ self._qy = copy.deepcopy(value)
+ else:
+ self._qx = beam_param(q=complex(value))
+ self._qy = beam_param(q=complex(value))
+
+ @property
+ def qx(self):
+ return self._qx.q
+
+ @qx.setter
+ def qx(self, value):
+ if value.__class__ == beam_param:
+ self._qx = copy.deepcopy(value)
+ else:
+ self._qx = beam_param(q=complex(value))
+
+ @property
+ def qy(self):
+ return self._qy.q
+ @qy.setter
+ def qy(self, value):
+ if value.__class__ == beam_param:
+ self._qy = copy.deepcopy(value)
+ else:
+ self._qy = beam_param(q=complex(value))
+
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.__xpre_const *= math.sqrt(1.0/(2**self._n * math.factorial(self._n)))
+ self.__xpre_const *= cmath.sqrt(1j*self._qx.imag / self._qx.q)
+ self.__xpre_const *= ((1j*self._qx.imag * self._qx.q.conjugate())/(-1j*self._qx.imag * self._qx.q)) ** ( 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):
- # need to finish...
- return self.__xpre_const
+ self.__ypre_const *= math.sqrt(1.0/(2**self._m * math.factorial(self._m)))
+ self.__ypre_const *= cmath.sqrt(1j*self._qy.imag / self._qy.q)
+ self.__ypre_const *= ((1j*self._qy.imag * self._qy.conjugate())/(-1j*self._qy.imag * self._qy.q)) **(self._m/2.0)
+
+ self.__sqrt2_wxz = math.sqrt(2) / self._qx.w
+ self.__sqrt2_wyz = math.sqrt(2) / self._qy.w
+
+ self.__kx = 2*math.pi / self._qx.wavelength
+ self.__ky = 2*math.pi / self._qy.wavelength
+
+ self.__invqx = 1/ self._qx.q
+ self.__invqy = 1/ self._qy.q
+
+ def _Un(self, x):
+ return self.__xpre_const * self._hn(self.__sqrt2_wxz * x) * np.exp(-0.5j * self.__kx * x*x * self.__invqx)
+
+ def _Um(self, y):
+ return self.__ypre_const * self._hm(self.__sqrt2_wyz * y) * np.exp(-0.5j * self.__ky * y*y * self.__invqy)
+
+ def Un(self, x):
+ vec = np.vectorize(self._Un, otypes=[np.complex64])
+ return vec(x=x)
+
+ def Um(self, y):
+ vec = np.vectorize(self._Um, otypes=[np.complex64])
+ return vec(y=y)
+
+ def _unm(self, x, y):
+ return self._Un(x) * self._Um(y)
+
+ def Unm(self, x, y):
+ vec = np.vectorize(self._unm, otypes=[np.complex64])
+ return vec(x=x,y=y)
+
+ def plot(self, ndx=100, ndy=100, xscale=4, yscale=4):
+ import pylab
+
+ xrange = xscale * np.linspace(-self._qx.w, self._qx.w, ndx)
+ yrange = yscale * np.linspace(-self._qy.w, self._qy.w, ndy)
+
+ dx = xrange[1]-xrange[0]
+ dy = yrange[1]-yrange[0]
+
+ xx, yy = np.meshgrid(xrange,yrange)
+
+ data = self.Unm(xx, yy)
+
+ fig = pylab.figure()
+ axes = pylab.imshow(np.abs(data), aspect=dx/dy, extent=[min(xrange),max(xrange),min(yrange),max(yrange)])
+ pylab.xlabel('x [m]')
+ pylab.ylabel('y [m]')
+ cbar = fig.colorbar(axes)
+ pylab.show()