From 4989df75a39d70a218c8dc7257fa1aa718f784b0 Mon Sep 17 00:00:00 2001
From: Daniel Brown <ddb@star.sr.bham.ac.uk>
Date: Mon, 16 Feb 2015 11:57:17 +0000
Subject: [PATCH] fixing some optivis node convention issues. Adding other
components and trying out aLIGO plotting, nearly working. Other ROM updates
too
---
bin/test_profile.py | 2 +
examples/LLO_matched.kat | 283 +++++++++++++++++++++++++++++++++++
examples/optivis_aligo_ex.py | 7 +
examples/optivis_ex.py | 21 ++-
pykat/components.py | 57 ++++++-
pykat/finesse.py | 5 +
pykat/optics/knm.py | 102 ++++++++-----
pykat/optics/maps.py | 30 ++--
pykat/optics/romhom.py | 233 ++++++++++++++++++----------
9 files changed, 591 insertions(+), 149 deletions(-)
create mode 100644 examples/LLO_matched.kat
create mode 100644 examples/optivis_aligo_ex.py
diff --git a/bin/test_profile.py b/bin/test_profile.py
index 5c420ff..588b68d 100644
--- a/bin/test_profile.py
+++ b/bin/test_profile.py
@@ -13,6 +13,8 @@ attr m2 Rc 2
pd circ n2
xaxis m1 phi lin 0 180 1000
yaxis log abs
+maxtem 0
+cav c1 m1 n1 m2 n2
""")
kat.timeCode = True
diff --git a/examples/LLO_matched.kat b/examples/LLO_matched.kat
new file mode 100644
index 0000000..7d54c8b
--- /dev/null
+++ b/examples/LLO_matched.kat
@@ -0,0 +1,283 @@
+%--------------------------------------------------------------------------
+%
+% Finesse model for the full dual recycled inteferometer, with Fabry-Perot
+% arm cavities and OMC.
+%
+% Currently plots the DARM coupled cavity transfer function for determining
+% cavity pole.
+%
+% Sources of 'real' numbers, unless otherwise stated:
+% - Mirror parameters (curvatures, losses etc.) from
+% https://galaxy.ligo.caltech.edu/optics/ If no number for the loss is
+% given we use 0.
+% - Distances between optics from LIGO-E1200274, the L1 master coordinate
+% list.
+%
+% 3/12/14 - Updated with Galaxy page values for installed SRM-w
+% - Includes locks for DARM (10pm offset), PRCL, MICH, CARM and SRCL
+% 4/12/14 - Added in more details about which Galaxy mirror data is taken from
+%
+% To-do:
+% - Check frequency and modulation index for side-bands
+% - AR surface for PRM
+%
+%
+# Charlotte Bond, Paul Fulda, Daniel Brown
+%--------------------------------------------------------------------------
+
+
+# %%% FTblock Laser
+# ###########################################################################
+# # Laser
+# l L0 2 0 n0
+# s lmod1 1 n0 n1
+# # Modulation indices updated from LLO logbook denis.martynov@LIGO.ORG - 02:01,
+# # Thursday 03 October 2013 (8940)
+# mod mod1 $f1 0.25 1 pm n1 n2
+# s lmod2 1 n2 n3
+# mod mod2 $f2 0.25 1 pm n3 nin
+# # Input beam mode-matched to case with 50km lens in itm substrates (matched
+# # to arm cavities) for LLO design case
+# gauss gparam L0 n0 1.3498243m 4.3433202 1.3619432m 4.5706208
+# ###########################################################################
+# %%% FTend Laser
+
+# %%% FTblock PR
+# ###########################################################################
+# # Distance to power recycling mirror
+# s lin 1 nin nREFL
+# # Power recycling mirror PRM-02
+# m1 PRM $T_PRM $L_PRM $phi_PRM nREFL nPRMb
+# attr PRM Rc 11.009
+# # Distance between PRM and PR2
+# s lp1 16.6107 nPRMb nPR2a
+# # PR2 PR2-02
+# bs1 PR2 $T_PR2 $L_PR2 0 -0.79 nPR2a nPR2b dump nPOP
+# attr PR2 Rc -4.545
+# # Distance from PR2 to PR3
+# s lp2 16.1647 nPR2b nPR3a
+# # PR3 PR3-03
+# bs1 PR3 $T_PR3 $L_PR3 0 0.615 nPR3a nPR3b dump dump
+# attr PR3 Rc 36.027
+# # Distance from PR3 to BS
+# s lp3 19.5381 nPR3b nHRBS_PR
+# ###########################################################################
+# %%% FTend PR
+
+l L0 2 0 n0
+s s1 100 n0 nHRBS_PR
+%%% FTblock BS
+###########################################################################
+# BS beamsplitter BS-02
+##------------------------------------------------------------
+## BS
+## ^
+## to IMY |
+## | ,'-.
+## | + `.
+## nYBS | ,' :'
+## nPR3b | +i1 +
+## ----------------> ,:._ i2 ,'
+## from the PRC nPRBS + \ `-. + nXBS
+## ,' i3\ ,' --------------->
+## + \ + to IMX
+## ,' i4.'
+## `._ ..
+## `._ ,' |nSRBS
+## - |
+## |to the SRC
+## |
+## v
+##------------------------------------------------------------
+# BS BS-02
+bs1 HRBS 0.5 8.6u $phi_BS 45 nHRBS_PR nHRBS_Y nHRBS_X nHRBS_SR
+s sHRBStoARBSX 5 $nsilica nHRBS_X nARBSX_sub
+bs2 ARBSX 30u 1.7u $phi_ARBSX 29.1951 nARBSX_sub dump nARBSX_X dump
+s sHRBStoARBSSR 5 $nsilica nHRBS_SR nARBSSR_sub
+bs2 ARBSSR 30u 1.7u $phi_ARBSSR -29.1951 nARBSSR_sub dump nARBSSR_SR dump
+###########################################################################
+%%% FTend BS
+
+
+# %%% FTblock SR
+# ###########################################################################
+# # Distance from BS to SR3
+# s ls3 19.3661 nARBSSR_SR nSR3a
+#
+# # SR3 SR3-01
+# bs1 SR3 $T_SR2 $L_SR3 0 0.785 nSR3a nSR3b dump dump
+# attr SR3 Rc 35.97
+#
+# # Distance from SR3 to SR2
+# s ls2 15.4435 nSR3b nSR2a
+#
+# # SR2 SR2-04
+# bs1 SR2 $T_SR2 $L_SR2 0 0.87 nSR2a nSR2b dump dump
+# attr SR2 Rc -6.406
+#
+# # Distance from SR2 to SRMHR
+# s ls1 15.7566 nSR2b nSRMHRa
+#
+# # Signal recycling mirror SRM-08
+# m1 SRMHR $T_SRM $L_SRM $phi_SRM nSRMHRa nSRMHRb
+# s SRMsub 0.0749 $nsilica nSRMHRb nSRMARa
+# attr SRMHR Rc -5.677
+# m2 SRMAR 50u 0 $phi_SRM nSRMARa nSRMARb
+# ###########################################################################
+# %%% FTend SR
+
+
+%%% FTblock Yarm
+###########################################################################
+# Using values from E1200274
+s ly1 4.847 nHRBS_Y nCPYar1
+# Y arm compensation plate CP-08
+m2 CPYar1 48.9u 0.4u 0 nCPYar1 nCPYar1s
+s sCPY 0.10032 $nsilica nCPYar1s nCPYar2s
+m2 CPYar2 30.5u 0.3u 0 nCPYar2s nCPYar2
+s sCPYtoITMYar 0.02 nCPYar2 nITMYTLin
+# Y arm input mirror ITM-08
+# Thermal lens
+lens ITMYTL $TL_f nITMYTLin nITMYTLtrans
+s ITMYTL_null 10 nITMYTLtrans nITMYconstL_in
+# Constant ITMY substrate lens
+lens ITMYconstL inf nITMYconstL_in nITMYconstL_trans
+# s ITMYTL_null2 0 nITMYconstL_trans nITMYar_in
+# m2 ITMYar 250u 0.6u $phi_ITMY nITMYar_in nITMYs1
+# s lITMY 0.19961 $nsilica nITMYs1 nITMYs2
+# m1 ITMY $T_ITMY $L_ITMY $phi_ITMY nITMYs2 nITMY2
+# attr ITMY Rc -1940.7
+# # Y-arm
+# s Ly 10 nITMY2 nETMY1
+# # ETMY
+# m1 ETMY $T_ETMY $L_ETMY $phi_ETMY nETMY1 nTY
+# attr ETMY Rc 2242.4
+###########################################################################
+%%% FTend Yarm
+
+# %%% FTblock Xarm
+# ###########################################################################
+# # Now using length taken from E1200616
+# s lx1 4.829 nARBSX_X nCPXar1
+# # X arm compensation plate CP-06
+# m2 CPXar1 33u 0.6u 0 nCPXar1 nCPXar1s
+# s sCPX 0.10031 $nsilica nCPXar1s nCPXar2s
+# m2 CPXar2 8u 0.6u 0 nCPXar2s nCPXar2
+# s sCPXtoITMXar 0.02 nCPXar2 nITMXTLin
+# # X arm input mirror ITM-04
+# # Thermal lens
+# lens ITMXTL $TL_f nITMXTLin nITMXTLtrans
+# s ITMXtl_null 0 nITMXTLtrans nITMXconstL_in
+# # Non-thermal ITM lens
+# lens ITMXconstL inf nITMXconstL_in nITMXconstL_trans
+# s ITMXTL_null2 0 nITMXconstL_trans nITMXar_in
+# m2 ITMXar 164u 0.5u $phi_ITMX nITMXar_in nITMXs1
+# s lITMX1 0.20027 $nsilica nITMXs1 nITMXs2
+# m1 ITMX $T_ITMX $L_ITMX $phi_ITMX nITMXs2 nITMX2
+# # default Rc from nebula page
+# attr ITMX Rc -1937.9
+# # X-arm
+# s Lx 10 nITMX2 nETMX1
+# # ETMX
+# m1 ETMX $T_ETMX $L_ETMX $phi_ETMX nETMX1 nTX
+# attr ETMX Rc 2239.7
+# ###########################################################################
+# %%% FTend Xarm
+
+
+%%% FTblock Reflectivities
+###########################################################################
+# Galaxy number 100
+const T_PR2 243u
+const L_PR2 8.6u
+
+# Galaxy number 105
+const T_PR3 5.3u
+const L_PR3 17u
+
+# Galaxy number 107
+const T_PRM 0.021
+const L_PRM 5.9u
+
+# Galaxy number 113
+const T_SR2 18.3u
+const L_SR2 6.1u
+
+# Galaxy number 114
+const T_SR3 5u
+const L_SR3 19.1u
+
+# Galaxy number 126
+const T_SRM 0.3688
+const L_SRM 8.3u
+
+# Galaxy number 75
+const T_ITMX 0.0148
+const L_ITMX 10.4u
+
+# Galaxy number 79
+const T_ITMY 0.0148
+const L_ITMY 14.3u
+
+# Galaxy number 30
+const T_ETMX 3.7u
+const L_ETMX 100.9u
+
+# Galaxy number 32
+const T_ETMY 3.6u
+const L_ETMY 9.3u
+###########################################################################
+%%% FTend Reflectivities
+
+
+%%% FTblock Constants
+###########################################################################
+const nsilica 1.44963098985906
+% Sidebands tuned to be resonant for PRC in this file (lprc = 57.6564)
+% Design sidebands
+%const f1 9099471
+%const mf1 -9099471
+%const f2 45497355
+%const mf2 -45497355
+% Measured?
+const f1 9099055
+const mf1 9099055
+const f2 45495275
+const mf2 -45495275
+const TL_f 34.5k
+###########################################################################
+%%% FTend Constants
+
+
+%%% FTblock Tunings
+###########################################################################
+# offset computed with zero_locks.py
+const phi_ETMX 0.106992111770048
+const phi_ETMY -0.109282418715425
+const phi_ITMX 0.106389383390564
+const phi_ITMY -0.106389383390564
+const phi_PRM -0.024603392339103
+const phi_SRM 90.108393822523524
+
+const phi_BS 0
+const phi_ARBSX 0
+const phi_ARBSSR 0
+const phi_IC 0
+###########################################################################
+%%% FTend Tunings
+
+
+%%% FTblock HOMs
+###########################################################################
+# Set PRY and PRX to use cavity eigenmodes for mode calculations
+cav SRCY SRMHR nSRMHRa ITMY nITMYs2
+cav SRCX SRMHR nSRMHRa ITMX nITMXs2
+cav PRCY PRM nPRMb ITMY nITMYs2
+cav PRCX PRM nPRMb ITMX nITMXs2
+cav Xarm ITMX nITMX2 ETMX nETMX1
+cav Yarm ITMY nITMY2 ETMY nETMY1
+cav OMC IC nICtrans IC nICout
+maxtem 4
+phase 2
+###########################################################################
+%%% FTend HOMs
diff --git a/examples/optivis_aligo_ex.py b/examples/optivis_aligo_ex.py
new file mode 100644
index 0000000..0e49704
--- /dev/null
+++ b/examples/optivis_aligo_ex.py
@@ -0,0 +1,7 @@
+import pykat
+
+kat = pykat.finesse.kat(kat_file="LLO_matched.kat")
+
+
+kat.optivis().show()
+
diff --git a/examples/optivis_ex.py b/examples/optivis_ex.py
index 130686d..e6638e0 100644
--- a/examples/optivis_ex.py
+++ b/examples/optivis_ex.py
@@ -4,20 +4,19 @@ kat = pykat.finesse.kat()
kat.parseCommands("""
l l1 1 0 n0
-s s1 100 n0 n1
-m m1 1 0 0 n1 n2
-s s1 50 n2 n3
+s s1 100 n0 n3
-bs bs1 1 0 45 45 n3 n4 n5 n6
+bs bs1 1 0 0 45 n3 n4 n5 n6
-s s2 200 n5 n5a
-m m2 1 0 0 n5a dump
+s s3 50 n4 n4a
+m m1 1 0 0 n4a n4b
+s s3a 200 n4b n7a
+m m2 1 0 0 n7a n7b
-s s3 200 n4 n4a
-m m3 1 0 0 n4a dump
-
-s s4 50 n6 n6a
-m m4 1 0 0 n6a n6b
+s s4 50 n5 n5a
+m m3 1 0 0 n5a n5b
+s s4a 200 n5b n8a
+m m4 1 0 0 n8a n8b
""")
kat.optivis().show()
diff --git a/pykat/components.py b/pykat/components.py
index 413a012..7593f1e 100644
--- a/pykat/components.py
+++ b/pykat/components.py
@@ -510,7 +510,7 @@ class beamSplitter(AbstractMirrorComponent):
def getOptivisComponent(self):
if self._optivis_component is None:
- self._optivis_component = optivis_components.BeamSplitter(name=self.name, aoi=self.alpha)
+ self._optivis_component = optivis_components.BeamSplitter(name=self.name, aoi=-self.alpha)
return self._optivis_component
@@ -526,14 +526,14 @@ class beamSplitter(AbstractMirrorComponent):
elif kat_node is self.nodes[1]:
return self._optivis_component.getInputNode("frB")
elif kat_node is self.nodes[2]:
- return self._optivis_component.getInputNode("bkA")
- elif kat_node is self.nodes[3]:
return self._optivis_component.getInputNode("bkB")
+ elif kat_node is self.nodes[3]:
+ return self._optivis_component.getInputNode("bkA")
elif mode == "output":
if kat_node is self.nodes[0]:
- return self._optivis_component.getOutputNode("frA")
- elif kat_node is self.nodes[1]:
return self._optivis_component.getOutputNode("frB")
+ elif kat_node is self.nodes[1]:
+ return self._optivis_component.getOutputNode("frA")
elif kat_node is self.nodes[2]:
return self._optivis_component.getOutputNode("bkA")
elif kat_node is self.nodes[3]:
@@ -915,7 +915,30 @@ class lens(Component):
rtn.extend(p.getFinesseText())
return rtn
+
+ def getOptivisComponent(self):
+ if self._optivis_component is None:
+ self._optivis_component = optivis_components.ConvexLens(name=self.name)
+ return self._optivis_component
+
+ def getOptivisNode(self, mode, kat_node):
+ mode = mode.lower()
+
+ if mode != "input" and mode.lower() != "output":
+ raise pkex.BasePyKatException("Mode must be either input or output")
+
+ if mode == "input":
+ if kat_node is self.nodes[0]:
+ return self._optivis_component.getInputNode("fr")
+ elif kat_node is self.nodes[1]:
+ return self._optivis_component.getInputNode("bk")
+ elif mode == "output":
+ if kat_node is self.nodes[0]:
+ return self._optivis_component.getnOutputNode("fr")
+ elif kat_node is self.nodes[1]:
+ return self._optivis_component.getOutputNode("bk")
+
def getQGraphicsItem(self):
if not USE_GUI:
raise NoGUIException
@@ -997,6 +1020,30 @@ class modulator(Component):
return rtn
+ def getOptivisComponent(self):
+ if self._optivis_component is None:
+ #self._optivis_component = optivis_components.Modulator(name=self.name)
+ self._optivis_component = optivis_components.ConvexLens(name=self.name)
+
+ return self._optivis_component
+
+ def getOptivisNode(self, mode, kat_node):
+ mode = mode.lower()
+
+ if mode != "input" and mode.lower() != "output":
+ raise pkex.BasePyKatException("Mode must be either input or output")
+
+ if mode == "input":
+ if kat_node is self.nodes[0]:
+ return self._optivis_component.getInputNode("fr")
+ elif kat_node is self.nodes[1]:
+ return self._optivis_component.getInputNode("bk")
+ elif mode == "output":
+ if kat_node is self.nodes[0]:
+ return self._optivis_component.getnOutputNode("fr")
+ elif kat_node is self.nodes[1]:
+ return self._optivis_component.getOutputNode("bk")
+
def getQGraphicsItem(self):
if not USE_GUI:
raise NoGUIException
diff --git a/pykat/finesse.py b/pykat/finesse.py
index 25a12df..20354e8 100644
--- a/pykat/finesse.py
+++ b/pykat/finesse.py
@@ -1480,6 +1480,11 @@ class kat(object):
a = c.connectingComponents()
+ # Need to handle where spaces don't connect two components but there is a loose
+ # node, which may or may not have detectors attached
+ if a[0] is None or a[1] is None:
+ continue
+
c1 = a[0].getOptivisComponent()
c2 = a[1].getOptivisComponent()
diff --git a/pykat/optics/knm.py b/pykat/optics/knm.py
index b75ad51..69435a8 100644
--- a/pykat/optics/knm.py
+++ b/pykat/optics/knm.py
@@ -5,6 +5,7 @@ from pykat.optics.romhom import u_star_u
from pykat.external.progressbar import ProgressBar, ETA, Percentage, Bar
from scipy.interpolate import interp2d
from scipy.integrate import dblquad
+from pykat.optics.romhom import ROMWeights
import time
import pykat.optics.maps
@@ -271,65 +272,84 @@ def ROM_HG_knm(weights, mode_in, mode_out, q1, q2, q1y=None, q2y=None, cache=Non
npr = mode_in[1]
mpr = mode_out[1]
- if cache == None:
- u_x_nodes = u_star_u(q1.z, q2.z, q1.w0, q2.w0, n, m, weights.EI["xm"].nodes)
- u_y_nodes = u_star_u(q1y.z, q2y.z, q1y.w0, q2y.w0, npr, mpr, weights.EI["ym"].nodes)
+ if isinstance(weights, ROMWeights):
+ if cache == None:
+ u_x_nodes = u_star_u(q1.z, q2.z, q1.w0, q2.w0, n, m, weights.EI["x"].nodes)
+ u_y_nodes = u_star_u(q1y.z, q2y.z, q1y.w0, q2y.w0, npr, mpr, weights.EI["y"].nodes)
- w_ij_Q1Q3 = weights.w_ij_Q1 + weights.w_ij_Q3
- w_ij_Q2Q4 = weights.w_ij_Q2 + weights.w_ij_Q4
- w_ij_Q1Q2 = weights.w_ij_Q1 + weights.w_ij_Q2
- w_ij_Q1Q4 = weights.w_ij_Q1 + weights.w_ij_Q4
- w_ij_Q2Q3 = weights.w_ij_Q2 + weights.w_ij_Q3
- w_ij_Q3Q4 = weights.w_ij_Q3 + weights.w_ij_Q4
- w_ij_Q1Q2Q3Q4 = weights.w_ij_Q1 + weights.w_ij_Q2 + weights.w_ij_Q3 + weights.w_ij_Q4
+ w_ij_Q1Q3 = weights.w_ij_Q1 + weights.w_ij_Q3
+ w_ij_Q2Q4 = weights.w_ij_Q2 + weights.w_ij_Q4
+ w_ij_Q1Q2 = weights.w_ij_Q1 + weights.w_ij_Q2
+ w_ij_Q1Q4 = weights.w_ij_Q1 + weights.w_ij_Q4
+ w_ij_Q2Q3 = weights.w_ij_Q2 + weights.w_ij_Q3
+ w_ij_Q3Q4 = weights.w_ij_Q3 + weights.w_ij_Q4
+ w_ij_Q1Q2Q3Q4 = weights.w_ij_Q1 + weights.w_ij_Q2 + weights.w_ij_Q3 + weights.w_ij_Q4
- else:
- strx = "x[%i,%i]" % (mode_in[0], mode_out[0])
- stry = "y[%i,%i]" % (mode_in[1], mode_out[1])
+ else:
+ strx = "x[%i,%i]" % (mode_in[0], mode_out[0])
+ stry = "y[%i,%i]" % (mode_in[1], mode_out[1])
- u_x_nodes = cache[strx]
- u_y_nodes = cache[stry]
+ u_x_nodes = cache[strx]
+ u_y_nodes = cache[stry]
- w_ij_Q1Q3 = cache["w_ij_Q1Q3"]
- w_ij_Q2Q4 = cache["w_ij_Q2Q4"]
- w_ij_Q1Q2 = cache["w_ij_Q1Q2"]
- w_ij_Q1Q4 = cache["w_ij_Q1Q4"]
- w_ij_Q2Q3 = cache["w_ij_Q2Q3"]
- w_ij_Q3Q4 = cache["w_ij_Q3Q4"]
- w_ij_Q1Q2Q3Q4 = cache["w_ij_Q1Q2Q3Q4"]
+ w_ij_Q1Q3 = cache["w_ij_Q1Q3"]
+ w_ij_Q2Q4 = cache["w_ij_Q2Q4"]
+ w_ij_Q1Q2 = cache["w_ij_Q1Q2"]
+ w_ij_Q1Q4 = cache["w_ij_Q1Q4"]
+ w_ij_Q2Q3 = cache["w_ij_Q2Q3"]
+ w_ij_Q3Q4 = cache["w_ij_Q3Q4"]
+ w_ij_Q1Q2Q3Q4 = cache["w_ij_Q1Q2Q3Q4"]
- u_xy_nodes = np.outer(u_x_nodes, u_y_nodes)
+ u_xy_nodes = np.outer(u_x_nodes, u_y_nodes)
- n_mod_2 = n % 2
- m_mod_2 = m % 2
- npr_mod_2 = npr % 2
- mpr_mod_2 = mpr % 2
+ n_mod_2 = n % 2
+ m_mod_2 = m % 2
+ npr_mod_2 = npr % 2
+ mpr_mod_2 = mpr % 2
- if n_mod_2 == m_mod_2 and npr_mod_2 == mpr_mod_2:
- k_ROQ = np.einsum('ij,ij', u_xy_nodes, w_ij_Q1Q2Q3Q4)
+ if n_mod_2 == m_mod_2 and npr_mod_2 == mpr_mod_2:
+ k_ROQ = np.einsum('ij,ij', u_xy_nodes, w_ij_Q1Q2Q3Q4)
- elif n_mod_2 != m_mod_2:
- if npr_mod_2 == mpr_mod_2:
- k_ROQ = np.einsum('ij,ij', u_xy_nodes, w_ij_Q1Q4) - np.einsum('ij,ij', u_xy_nodes, w_ij_Q2Q3)
- else:
- k_ROQ = np.einsum('ij,ij', u_xy_nodes, w_ij_Q2Q4) - np.einsum('ij,ij', u_xy_nodes, w_ij_Q1Q3)
+ elif n_mod_2 != m_mod_2:
+ if npr_mod_2 == mpr_mod_2:
+ k_ROQ = np.einsum('ij,ij', u_xy_nodes, w_ij_Q1Q4) - np.einsum('ij,ij', u_xy_nodes, w_ij_Q2Q3)
+ else:
+ k_ROQ = np.einsum('ij,ij', u_xy_nodes, w_ij_Q2Q4) - np.einsum('ij,ij', u_xy_nodes, w_ij_Q1Q3)
- elif npr_mod_2 != mpr_mod_2:
- if n_mod_2 == m_mod_2:
- k_ROQ = np.einsum('ij,ij', u_xy_nodes, w_ij_Q3Q4) - np.einsum('ij,ij', u_xy_nodes, w_ij_Q1Q2)
+ elif npr_mod_2 != mpr_mod_2:
+ if n_mod_2 == m_mod_2:
+ k_ROQ = np.einsum('ij,ij', u_xy_nodes, w_ij_Q3Q4) - np.einsum('ij,ij', u_xy_nodes, w_ij_Q1Q2)
+ else:
+ k_ROQ = np.einsum('ij,ij', u_xy_nodes, w_ij_Q2Q4) - np.einsum('ij,ij', u_xy_nodes, w_ij_Q1Q3)
+
+ else:
+ if cache == None:
+ u_x_nodes = u_star_u(q1.z, q2.z, q1.w0, q2.w0, n, m, weights.EI["x"].nodes)
+ u_y_nodes = u_star_u(q1y.z, q2y.z, q1y.w0, q2y.w0, npr, mpr, weights.EI["y"].nodes)
+
+ w_ij = weights.w_ij
else:
- k_ROQ = np.einsum('ij,ij', u_xy_nodes, w_ij_Q2Q4) - np.einsum('ij,ij', u_xy_nodes, w_ij_Q1Q3)
+ strx = "x[%i,%i]" % (mode_in[0], mode_out[0])
+ stry = "y[%i,%i]" % (mode_in[1], mode_out[1])
+
+ u_x_nodes = cache[strx]
+ u_y_nodes = cache[stry]
-
+ u_xy_nodes = np.outer(u_x_nodes, u_y_nodes)
+
+ k_ROQ = np.einsum('ij,ij', u_xy_nodes, w_ij)
+
return k_ROQ
__fac_cache = []
def fac(n):
global __fac_cache
- return __fac_cache[n]
- #return math.factorial(int(n))
+ if len(__fac_cache) == 0:
+ return math.factorial(int(n))
+ else:
+ return __fac_cache[n]
def m_1_pow(n):
if n % 2 == 0:
diff --git a/pykat/optics/maps.py b/pykat/optics/maps.py
index 0eb5422..a399727 100644
--- a/pykat/optics/maps.py
+++ b/pykat/optics/maps.py
@@ -143,7 +143,7 @@ class surfacemap(object):
raise BasePyKatException("Map type needs handling")
- def generateROMWeights(self, isModeMatched=True, verbose=False, interpolate=False, interpolate_N=None, tolerance = 1e-12, sigma = 1, sort=False, greedyfile=None):
+ def generateROMWeights(self, isModeMatched=True, verbose=False, interpolate=False, interpolate_N=None, tolerance = 1e-12, sigma = 1, sort=False, greedyfile=None, useSymmetry=True):
if interpolate:
from scipy.interpolate import interp2d
@@ -175,27 +175,31 @@ class surfacemap(object):
self.center = (np.array(data.shape)+1)/2.0
self.data = data
-
- xm = self.x[self.x<0]
- ym = self.y[self.y<0]
-
- symm = False
+
+ if useSymmetry:
+ xm = self.x[self.x<0]
+ ym = self.y[self.y<0]
+ else:
+ xm = self.x
+ ym = self.y
if min(xm) == min(ym) and max(xm) == max(ym) and len(xm) == len(ym):
symm = True
-
+ else:
+ symm = False
+
EI = {}
- EI["xm"] = makeEmpiricalInterpolant(makeReducedBasis(xm, isModeMatched=isModeMatched, tolerance = tolerance, sigma = sigma, greedyfile=greedyfile), sort=sort)
+ EI["x"] = makeEmpiricalInterpolant(makeReducedBasis(xm, isModeMatched=isModeMatched, tolerance = tolerance, sigma = sigma, greedyfile=greedyfile), sort=sort)
if symm:
- EI["ym"] = EI["xm"]
+ EI["y"] = EI["x"]
else:
- EI["ym"] = makeEmpiricalInterpolant(makeReducedBasis(ym, isModeMatched=isModeMatched, tolerance = tolerance, sigma = sigma, greedyfile=greedyfile), sort=sort)
+ EI["y"] = makeEmpiricalInterpolant(makeReducedBasis(ym, isModeMatched=isModeMatched, tolerance = tolerance, sigma = sigma, greedyfile=greedyfile), sort=sort)
- EI["limits"] = EI["xm"].limits
+ EI["limits"] = EI["x"].limits
- self._rom_weights = makeWeights(self, EI, verbose=verbose)
+ self._rom_weights = makeWeights(self, EI, verbose=verbose, useSymmetry=useSymmetry)
return self.ROMWeights, EI
@@ -306,7 +310,7 @@ class tiltmap(surfacemap):
xx, yy = np.meshgrid(self.x, self.y)
- self.data = (xx * self.tilt[1] + yy * self.tilt[0])/self.scaling
+ self.data = (yy * self.tilt[1] + xx * self.tilt[0])/self.scaling
class zernikemap(surfacemap):
diff --git a/pykat/optics/romhom.py b/pykat/optics/romhom.py
index 1766b83..958e3e7 100644
--- a/pykat/optics/romhom.py
+++ b/pykat/optics/romhom.py
@@ -41,14 +41,14 @@ class ROMWeights:
f.write("w0max=%16.16e\n" % self.limits.w0max)
f.write("maxorder=%i\n" % self.limits.max_order)
- f.write("xnodes=%i\n" % len(self.EI["xm"].nodes))
+ f.write("xnodes=%i\n" % len(self.EI["x"].nodes))
- for v in self.EI["xm"].nodes.flatten():
+ for v in self.EI["x"].nodes.flatten():
f.write("%s\n" % repr(float(v)))
- f.write("ynodes=%i\n" % len(self.EI["ym"].nodes))
+ f.write("ynodes=%i\n" % len(self.EI["y"].nodes))
- for v in self.EI["ym"].nodes.flatten():
+ for v in self.EI["y"].nodes.flatten():
f.write("%s\n" % repr(float(v)))
f.write(repr(self.w_ij_Q1.shape) + "\n")
@@ -73,6 +73,45 @@ class ROMWeights:
f.close()
+class ROMWeightsFull:
+
+ def __init__(self, w_ij, EI, limits):
+ self.w_ij = w_ij
+
+ self.EI = EI
+ self.limits = limits
+
+ def writeToFile(self, filename):
+ """
+ Writes this map's ROM weights to a file
+ that can be used with Finesse. The filename
+ is appended with '.rom' internally.
+ """
+ f = open(filename + ".rom", 'w+')
+
+ f.write("zmin=%16.16e\n" % self.limits.zmin)
+ f.write("zmax=%16.16e\n" % self.limits.zmax)
+ f.write("w0min=%16.16e\n" % self.limits.w0min)
+ f.write("w0max=%16.16e\n" % self.limits.w0max)
+ f.write("maxorder=%i\n" % self.limits.max_order)
+
+ f.write("xnodes=%i\n" % len(self.EI["x"].nodes))
+
+ for v in self.EI["x"].nodes.flatten():
+ f.write("%s\n" % repr(float(v)))
+
+ f.write("ynodes=%i\n" % len(self.EI["y"].nodes))
+
+ for v in self.EI["y"].nodes.flatten():
+ f.write("%s\n" % repr(float(v)))
+
+ f.write(repr(self.w_ij.shape) + "\n")
+
+ for v in self.w_ij.flatten():
+ f.write("%s\n" % repr(complex(v))[1:-1])
+
+ f.close()
+
def combs(a, r):
"""
Return successive r-length combinations of elements in the array a.
@@ -273,100 +312,136 @@ def makeEmpiricalInterpolant(RB, sort=False):
return EmpiricalInterpolant(B=B, nodes=nodes, node_indices=node_indices, limits=RB.limits, x=RB.x)
-def makeWeights(smap, EI, verbose=True):
+def makeWeights(smap, EI, verbose=True, useSymmetry=True):
- # get full A_xy
- A_xy = smap.z_xy()
+ if useSymmetry:
+ # get full A_xy
+ A_xy = smap.z_xy().transpose()
- xm = smap.x[smap.x < 0]
- xp = smap.x[smap.x > 0]
+ xm = smap.x[smap.x < 0]
+ xp = smap.x[smap.x > 0]
- ym = smap.y[smap.y < 0]
- yp = smap.y[smap.y > 0]
+ ym = smap.y[smap.y < 0]
+ yp = smap.y[smap.y > 0]
- Q1xy = np.ix_(smap.x < 0, smap.y > 0)
- Q2xy = np.ix_(smap.x > 0, smap.y > 0)
- Q3xy = np.ix_(smap.x > 0, smap.y < 0)
- Q4xy = np.ix_(smap.x < 0, smap.y < 0)
+ Q1xy = np.ix_(smap.x < 0, smap.y > 0)
+ Q2xy = np.ix_(smap.x > 0, smap.y > 0)
+ Q3xy = np.ix_(smap.x > 0, smap.y < 0)
+ Q4xy = np.ix_(smap.x < 0, smap.y < 0)
- # get A_xy in the four quadrants of the x-y plane
- A_xy_Q1 = A_xy[Q1xy]
- A_xy_Q2 = A_xy[Q2xy]
- A_xy_Q3 = A_xy[Q3xy]
- A_xy_Q4 = A_xy[Q4xy]
+ # get A_xy in the four quadrants of the x-y plane
+ A_xy_Q1 = A_xy[Q1xy]
+ A_xy_Q2 = A_xy[Q2xy]
+ A_xy_Q3 = A_xy[Q3xy]
+ A_xy_Q4 = A_xy[Q4xy]
- # Not sure if there is a neater way to select the x=0,y=0 cross elements
- A_xy_0 = np.hstack([A_xy[:,smap.x==0].flatten(), A_xy[smap.y==0,:].flatten()])
+ # Not sure if there is a neater way to select the x=0,y=0 cross elements
+ A_xy_0 = np.hstack([A_xy[:,smap.x==0].flatten(), A_xy[smap.y==0,:].flatten()])
- full_x = smap.x
- full_y = smap.y
+ full_x = smap.x
+ full_y = smap.y
- dx = full_x[1] - full_x[0]
- dy = full_y[1] - full_y[0]
+ dx = full_x[1] - full_x[0]
+ dy = full_y[1] - full_y[0]
- if verbose:
- count = 4*len(EI["xm"].B) * len(EI["ym"].B)
- p = ProgressBar(maxval=count, widgets=["Computing weights: ", Percentage(), Bar(), ETA()])
+ if verbose:
+ count = 4*len(EI["x"].B) * len(EI["y"].B)
+ p = ProgressBar(maxval=count, widgets=["Computing weights: ", Percentage(), Bar(), ETA()])
- n = 0
+ n = 0
- # make integration weights
- Bx = EI["xm"].B
- By = EI["ym"].B[:,::-1]
- w_ij_Q1 = np.zeros((len(Bx),len(By)), dtype = complex)
+ # make integration weights
+ Bx = EI["x"].B
+ By = EI["y"].B[:,::-1]
+ w_ij_Q1 = np.zeros((len(Bx),len(By)), dtype = complex)
- for i in range(len(Bx)):
- for j in range(len(By)):
- B_ij_Q1 = np.outer(Bx[i], By[j])
- w_ij_Q1[i][j] = dx*dy*np.einsum('ij,ij', B_ij_Q1, A_xy_Q1)
+ for i in range(len(Bx)):
+ for j in range(len(By)):
+ B_ij_Q1 = np.outer(Bx[i], By[j])
+ w_ij_Q1[i][j] = dx*dy*np.einsum('ij,ij', B_ij_Q1, A_xy_Q1)
- if verbose:
- p.update(n)
- n+=1
-
- Bx = EI["xm"].B[:,::-1]
- By = EI["ym"].B[:,::-1]
- w_ij_Q2 = np.zeros((len(Bx),len(By)), dtype = complex)
-
- for i in range(len(Bx)):
- for j in range(len(By)):
- B_ij_Q2 = np.outer(Bx[i], By[j])
- w_ij_Q2[i][j] = dx*dy*np.einsum('ij,ij', B_ij_Q2, A_xy_Q2)
+ if verbose:
+ p.update(n)
+ n+=1
+
+ Bx = EI["x"].B[:,::-1]
+ By = EI["y"].B[:,::-1]
+ w_ij_Q2 = np.zeros((len(Bx),len(By)), dtype = complex)
+
+ for i in range(len(Bx)):
+ for j in range(len(By)):
+ B_ij_Q2 = np.outer(Bx[i], By[j])
+ w_ij_Q2[i][j] = dx*dy*np.einsum('ij,ij', B_ij_Q2, A_xy_Q2)
- if verbose:
- p.update(n)
- n+=1
+ if verbose:
+ p.update(n)
+ n+=1
+
+ Bx = EI["x"].B[:,::-1]
+ By = EI["y"].B
+ w_ij_Q3 = np.zeros((len(Bx),len(By)), dtype = complex)
+
+ for i in range(len(Bx)):
+ for j in range(len(By)):
+ B_ij_Q3 = np.outer(Bx[i], By[j])
+ w_ij_Q3[i][j] = dx*dy*np.einsum('ij,ij', B_ij_Q3, A_xy_Q3)
+
+ if verbose:
+ p.update(n)
+ n+=1
+
+ Bx = EI["x"].B
+ By = EI["y"].B
+ w_ij_Q4 = np.zeros((len(Bx),len(By)), dtype = complex)
+
+ for i in range(len(Bx)):
+ for j in range(len(By)):
+ B_ij_Q4 = np.outer(Bx[i], By[j])
+ w_ij_Q4[i][j] = dx*dy*np.einsum('ij,ij', B_ij_Q4, A_xy_Q4)
+
+ if verbose:
+ p.update(n)
+ n+=1
+ 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, EI=EI, limits=EI["x"].limits)
+
+ else:
+ # get full A_xy
+ A_xy = smap.z_xy()
- Bx = EI["xm"].B[:,::-1]
- By = EI["ym"].B
- w_ij_Q3 = np.zeros((len(Bx),len(By)), dtype = complex)
-
- for i in range(len(Bx)):
- for j in range(len(By)):
- B_ij_Q3 = np.outer(Bx[i], By[j])
- w_ij_Q3[i][j] = dx*dy*np.einsum('ij,ij', B_ij_Q3, A_xy_Q3)
+ full_x = smap.x
+ full_y = smap.y
- if verbose:
- p.update(n)
- n+=1
+ dx = full_x[1] - full_x[0]
+ dy = full_y[1] - full_y[0]
- Bx = EI["xm"].B
- By = EI["ym"].B
- w_ij_Q4 = np.zeros((len(Bx),len(By)), dtype = complex)
-
- for i in range(len(Bx)):
- for j in range(len(By)):
- B_ij_Q4 = np.outer(Bx[i], By[j])
- w_ij_Q4[i][j] = dx*dy*np.einsum('ij,ij', B_ij_Q4, A_xy_Q4)
+ if verbose:
+ count = len(EI["x"].B) * len(EI["y"].B)
+ p = ProgressBar(maxval=count, widgets=["Computing weights: ", Percentage(), Bar(), ETA()])
- if verbose:
- p.update(n)
- n+=1
- if verbose:
- p.finish()
+ n = 0
+
+ # make integration weights
+ Bx = EI["x"].B
+ By = EI["y"].B
- 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, EI=EI, limits=EI["xm"].limits)
-
+ w_ij = np.zeros((len(Bx), len(By)), dtype = complex)
+
+ for i in range(len(Bx)):
+ for j in range(len(By)):
+ B_ij = np.outer(Bx[i], By[j])
+ w_ij[i][j] = dx*dy*np.einsum('ij,ij', B_ij, A_xy)
+ if verbose:
+ p.update(n)
+ n+=1
+
+ if verbose:
+ p.finish()
+
+ return ROMWeightsFull(w_ij=w_ij, EI=EI, limits=EI["x"].limits)
+
--
GitLab