Commit 4989df75 authored by Daniel Brown's avatar Daniel Brown

fixing some optivis node convention issues. Adding other components and trying...

fixing some optivis node convention issues. Adding other components and trying out aLIGO plotting, nearly working. Other ROM updates too
parent 8df654b8
......@@ -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
......
%--------------------------------------------------------------------------
%
% 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
import pykat
kat = pykat.finesse.kat(kat_file="LLO_matched.kat")
kat.optivis().show()
......@@ -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()
......
......@@ -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
......
......@@ -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()
......
......@@ -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:
......
......@@ -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):
......
......@@ -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]