ifo.py 47.4 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
from __future__ import absolute_import
from __future__ import division
from __future__ import print_function
from __future__ import unicode_literals

import numpy as np
import math
import copy
import warnings
import cmath
11
12
13
import inspect
import six 

14
from pykat import finesse
15
from pykat.finesse import BlockedKatFile
16

17
18
import pykat.components
import pykat.exceptions as pkex
Andreas Freise's avatar
Andreas Freise committed
19
20
import pykat.external.peakdetect as peak
import matplotlib.pyplot as plt
21
import pkg_resources
22

23
from scipy.optimize import fmin
24

25
26
27
global nsilica, clight
nsilica = 1.44963098985906
clight = 299792458.0
28

29
class aLIGO(object):
30
31
32
    """
    Object storing the Finesse input file and some auxilliary information
    about Avanced LIGO interferometers.
33
34
35
36
37
38
39
40
41
42
43
44
45

    References:
    [1] C. Bond `How to stay in shape: Overcoming beam and mirror distortions
        in advanced gravitational wave interferometers', PhD thesis, University
        of Birmingham, 2014, http://etheses.bham.ac.uk/5223/2/Bond14PhD.pdf
    [2] D. Martynov, `Lock Acquisition and Sensitivity Analysis of Advanced
        LIGO Interferometers', PhD thesis, Caltech, 2015
        http://thesis.library.caltech.edu/8899/1/DenisMartynovThesis.pdf
    [3] A. Staley, `Locking the Advanced LIGO Gravitational Wave Detector:
        with a focus on the Arm Length Stabilization Technique', PhD thesis,
        Columbia University, 2015
        https://academiccommons.columbia.edu/catalog/ac:189457
        
46
47
    """

48
    def __init__(self, _name="default", katfile=None, verbose = False, debug=False):
49
        names = ['default', 'LLO', 'LHO']
50
51
52
53
        if debug:
            self.kat = finesse.kat(tempdir=".",tempname="test")
        else:
            self.kat = finesse.kat()
54
        self.kat.verbose=verbose
Andreas Freise's avatar
Andreas Freise committed
55
56
        self._data_path=pkg_resources.resource_filename('pykat.gw_detectors','finesse_files/')

57
        self.rawBlocks = BlockedKatFile()
58
59
        if katfile:
            self.kat.loadKatFile(katfile)
60
            self.rawBlocks.read(katfile)
61
62
63
64
        else:
            if _name not in names: # TODO different files not yet implemented
                printf("aLIGO name `{}' not recognised, must be 'default', 'LLO' or 'LHO'",_name)
            #print(data_path)
Andreas Freise's avatar
Andreas Freise committed
65
            self.kat.loadKatFile(self._data_path+"aLIGO.kat")
66
            self.rawBlocks.read(self._data_path+"aLIGO.kat")
67

68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
        # ----------------------------------------------------------------------
        # set variables to zero first
        self.DCoffset = 0.0
        self.DCoffsetW = 0.0

        # ----------------------------------------------------------------------
        # get and derive parameters from the kat file

        # get main frequencies
        if "f1" in self.kat.constants.keys():
            self.f1 = float(self.kat.constants["f1"].value)
        else:
            self.f1 = 9099471.0
        if "f2" in self.kat.constants.keys():
            self.f2 = float(self.kat.constants["f2"].value)
        else:
            self.f2 = 5.0 * self.f1
        if "f3" in self.kat.constants.keys():
            self.f3 = float(self.kat.constants["f3"].value)
        # TODO add else here!
88
89
90
        # check modultion frequencies
        if (5 * self.f1 != self.f2):
            print(" ** Warning: modulation frequencies do not match: 5*f1!=f2")
91
92
93
94
95
        
        # defining a dicotionary for the main mirror positions (tunings)
        self.tunings = {}
        self.tunings = self.get_tunings(self.kat)
        self.compute_derived_lengths(self.kat)
96

97
            
98
99
100
        # ----------------------------------------------------------------------
        # define ports and signals 
        
101
        # useful ports
102
103
104
105
        self.POP_f1  = port("POP_f1",   "nPOP",  self.f1, phase=101)
        self.POP_f2  = port("POP_f2",   "nPOP",  self.f2, phase=13)
        self.REFL_f1 = port("REFL_f1",  "nREFL", self.f1, phase=101)
        self.REFL_f2 = port("REFL_f2",  "nREFL", self.f2, phase=14)
106
107
108
109
110
111
112
113
        self.AS_DC   = port("AS_DC", "nSRM2")
        self.POW_BS  = port("PowBS", "nPRBS*")
        self.POW_X   = port("PowX",  "nITMX2")
        self.POW_Y   = port("PowY",  "nITMY2")

        # pretune DOF
        self.preARMX =  DOF("ARMX", self.POW_X,   "", "ETMX", 1, 1.0)
        self.preARMY =  DOF("ARMY", self.POW_Y,   "", "ETMY", 1, 1.0)
114
        self.preMICH =  DOF("AS"  , self.AS_DC,   "", ["ITMX", "ETMX", "ITMY", "ETMY"], [1,1,-1,-1], 6.0)
115
        self.prePRCL =  DOF("PRCL", self.POW_BS,  "", "PRM",  1, 10.0)
116
        self.preSRCL =  DOF("SRCL", self.AS_DC,   "", "SRM",  1, 10.0)
117
        
118
        # control scheme as in [1] Table C.1  
119
        self.PRCL =  DOF("PRCL", self.POP_f1,  "I", "PRM", 1, 100.0)
120
        self.MICH =  DOF("MICH", self.POP_f2,  "Q", ["ITMX", "ETMX", "ITMY", "ETMY"], [1,1,-1,-1], 100.0)
121
122
123
        self.CARM =  DOF("CARM", self.REFL_f1, "I", ["ETMX", "ETMY"], [1, 1], 1.5)
        self.DARM =  DOF("DARM", self.AS_DC,   "",  ["ETMX", "ETMY"], [1,-1], 1.0)
        self.SRCL =  DOF("SRCL", self.REFL_f2, "I", "SRM", 1, 1e2)
124
125
126
127
128
129
130
131
132
133
        
        self.__DOFs = {}
        
        for _ in inspect.getmembers(self, lambda x: isinstance(x, DOF)):
            self.__DOFs[_[0]] = _[1]
    
    @property
    def DOFs(self):
        return copy.copy(self.__DOFs)
    
134
    def adjust_PRC_length(self, kat, verbose=False):
135
136
137
138
139
        """
        Adjust PRC length so that it fulfils the requirement
        lPRC = (N+1/2) * c/(2*f1), see [1] equation C.1
        In the current design N=3
        """
140
        vprint(verbose, "-- adjusting PRC length")
141
142
        ltmp = 0.5 * clight / self.f1
        delta_l = 3.5 * ltmp - self.lPRC
143
        vprint(verbose, "   adusting kat.lp1.L by {:.4g}m".format(delta_l))
144
145
146
147
148
149
150
151
152
153
154
155
        kat.lp1.L += delta_l
        self.compute_derived_lengths(kat)

    def check_f1_PRC_resonance(self, _kat):
        """
        Plot the sideband amplitudes for mdulation frequecy
        f1 (~ 9MHz) in the PRC, to check the resonance
        condition.
        """
        kat = _kat.deepcopy()
        fig = plt.figure()
        ax = fig.add_subplot(1,1,1)
156
157
158
159
160
161
        startf = self.f1-4000.0
        stopf  = self.f1+4000.0
        if _kat.maxtem == "off":
            nmStr = None
        else:
            nmStr = "0 0"
162
        code = """
163
164
165
166
167
168
ad f1p {0} {1} nPRM2
ad f1m {0} -{1} nPRM2
xaxis mod1 f lin {2} {3} 200
put f1p f $x1 
put f1m f $mx1 
""".format(nmStr, self.f1, startf, stopf)
169
170
171
172
173
174
175
        kat.parseCommands(code)
        out = kat.run()
        ax.plot(out.x-self.f1,np.abs(out["f1p"]), label=" f1")
        ax.plot(out.x-self.f1,np.abs(out["f1m"]), label="-f1")
        ax.set_xlim([np.min(out.x-self.f1), np.max(out.x-self.f1)])
        ax.set_xlabel("delta_f1 [Hz]")
        ax.set_ylabel('sqrt(W) ')
176
        ax.grid(True)
177
178
179
180
181
182
183
184
185
186
187
188
189
190
        ax.legend()
        plt.tight_layout()
        plt.show(block=0)

        
    def compute_derived_lengths(self, kat, verbose=False):
        """
        Compute derived length from individual space components.
        Design values are currently:
        lPRC = 57.656, lSRC = 56.008, lSchnupp = 0.08
        and the individual lengths:
        PRC: Lp1 16.6107, Lp2 16.1647, Lp3 19.5381 
        SRC: Ls1 15.7586, Ls2 15.4435, Ls3 19.3661
        """
191
192
        # distances between HR surfaces:
        self.lpr = kat.lp1.L + kat.lp2.L + kat.lp3.L
193
194
        self.lx = kat.lx1.L + kat.BSsub1.L * kat.BSsub1.n + kat.ITMXsub.L * kat.ITMXsub.n
        self.ly = kat.ly1.L + kat.ITMYsub.L * kat.ITMYsub.n
195
196
197
198
199
200
        self.lsr = kat.ls1.L + kat.ls2.L + kat.ls3.L + kat.BSsub2.L * kat.BSsub2.n
        # resulting combined distances (single, not roundtrip)
        self.lMI =  0.5 * (self.lx + self.ly)
        self.lPRC = self.lpr + self.lMI
        self.lSRC = self.lsr + self.lMI
        self.lSchnupp = self.lx - self.ly
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
        self.compute_derived_resonances(kat)

    def compute_derived_resonances(self, kat):
        self.fsrX = 0.5 * clight / float(kat.LX.L)
        self.fsrY = 0.5 * clight / float(kat.LY.L)
        self.fsrPRC = 0.5 * clight / self.lPRC
        self.fsrSRC = 0.5 * clight / self.lSRC
        self.f1_PRC = 3.5 * self.fsrPRC
        
    def lengths_status(self, kat):
        print(" .--------------------------------------------------.")
        print("| - arm length:                                     |")
        print("| Lx   = {:11.7}m, Ly   = {:11.7}m          |".format(float(kat.LX.L), float(kat.LY.L)))
        print("| - small MI and recycling lengths:                 | ")
        print("| lx   = {:11.7}m, ly   = {:11.7}m          |".format(self.lx, self.ly))
        print("| lpr  = {:11.7}m, lsr  = {:11.7}m          |".format(self.lpr, self.lsr))
        print("| lMI  = {:11.7}m, lSchnupp = {:11.5}m      |".format(self.lMI, self.lSchnupp))
        print("| lPRC = {:11.7}m, lSRC = {:11.7}m          |".format(self.lPRC, self.lSRC))
        print("+---------------------------------------------------+")
        print("| - associated cavity frequencies [Hz]:             |")
        print("| fsrx   = {:11.8}, fsry   = {:11.8}        |".format(self.fsrX, self.fsrY))
        print("| fsrPRC = {:11.8}, fsrSRC = {:11.8}        |".format(self.fsrPRC, self.fsrSRC))
        print("| f1_PRC = {:11.8}                              |".format(self.f1_PRC))
        print("| f1     = {:11.8}, f2     = {:11.9}        |".format(self.f1, self.f2))
        print(" `--------------------------------------------------'")
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
        
    def get_tunings(self, kat):
        self.tunings["maxtem"] = kat.maxtem
        self.tunings["PRM"]    = kat.PRM.phi
        self.tunings["ITMX"]   = kat.ITMX.phi
        self.tunings["ETMX"]   = kat.ETMX.phi
        self.tunings["ITMY"]   = kat.ITMY.phi
        self.tunings["ETMY"]   = kat.ETMY.phi
        self.tunings["BS"]     = kat.BS.phi
        self.tunings["SRM"]    = kat.SRM.phi
        return self.tunings
    
    def set_tunings(self, kat, tunings):
        kat.maxtem   = tunings["maxtem"]
        kat.PRM.phi  = tunings["PRM"]  
        kat.ITMX.phi = tunings["ITMX"]  
        kat.ETMX.phi = tunings["ETMX"]  
        kat.ITMY.phi = tunings["ITMY"]  
        kat.ETMY.phi = tunings["ETMY"]  
        kat.BS.phi   = tunings["BS"]  
        kat.SRM.phi  = tunings["SRM"]  

248
249
    def apply_lock_feedback(self, kat, out):
        tuning = self.get_tunings(kat)
250
        if "ETMX_lock" in out.ylabels:
251
252
253
            tuning["ETMX"] += float(out["ETMX_lock"])
        else:
            print(" ** Warning: could not find ETMX lock")
254
        if "ETMY_lock" in out.ylabels:
255
256
257
            tuning["ETMY"] += float(out["ETMY_lock"])
        else:
            print(" ** Warning: could not find ETMY lock")
258
        if "PRCL_lock" in out.ylabels:
259
260
261
            tuning["PRM"]  += float(out["PRCL_lock"])
        else:
            print(" ** Warning: could not find PRCL lock")
262
        if ("MICH_lock" in out.ylabels) and ("ITMY_lock" in out.ylabels):
263
264
265
266
            tuning["ITMX"] += float(out["MICH_lock"])
            tuning["ITMY"] += float(out["ITMY_lock"])
        else:
            print(" ** Warning: could not find MICH (ITMY) lock")
267
        if "SRCL_lock" in out.ylabels:
268
269
270
            tuning["SRM"]  += float(out["SRCL_lock"])
        else:
            print(" ** Warning: could not find SRCL lock")
271
        self.set_tunings(kat, tuning)
272
273

        
274
    def pretune(self, _kat, pretune_precision=1.0e-4, verbose=False):
275
276
        print("-- pretuning interferometer to precision {0:2g} deg = {1:2g} m".format(pretune_precision, pretune_precision*_kat.lambda0/360.0))
        kat=_kat.deepcopy()
277
278
        #_kat.ITMX.phi=0.0
        #_kat.ITMY.phi=0.0
279
        vprint(verbose, "   scanning X arm (maximising power)")
280
281
282
283
284
        kat1 = kat.deepcopy()
        make_transparent(kat1,["PRM","SRM"])
        make_transparent(kat1,["ITMY", "ETMY"])
        kat1.BS.setRTL(0.0,1.0,0.0) # set BS refl. for X arm
        phi, precision = self.scan_to_precision(kat1, self.preARMX, pretune_precision)
285
286
        phi=round(phi/pretune_precision)*pretune_precision
        phi=round_to_n(phi,5)
287
        vprint(verbose, "   found max/min at: {} (precision = {:2g})".format(phi, precision))
288
        self.preARMX.apply_tuning(_kat,phi)
289
    
290
        vprint(verbose, "   scanning Y arm (maximising power)")
291
292
293
294
295
296
297
        kat = _kat.deepcopy()
        make_transparent(kat,["PRM","SRM"])
        make_transparent(kat,["ITMX", "ETMX"])
        kat.BS.setRTL(1.0,0.0,0.0) # set BS refl. for Y arm
        phi, precision = self.scan_to_precision(kat, self.preARMY, pretune_precision)
        phi=round(phi/pretune_precision)*pretune_precision
        phi=round_to_n(phi,5)
298
        vprint(verbose, "   found max/min at: {} (precision = {:2g})".format(phi, precision))
299
        self.preARMY.apply_tuning(_kat,phi)
300
    
301
        vprint(verbose, "   scanning MICH (minimising power)")
302
303
304
305
306
        kat = _kat.deepcopy()
        make_transparent(kat,["PRM","SRM"])
        phi, precision = self.scan_to_precision(kat, self.preMICH, pretune_precision, minmax="min", precision=30.0)
        phi=round(phi/pretune_precision)*pretune_precision
        phi=round_to_n(phi,5)
307
        vprint(verbose, "   found max/min at: {} (precision = {:2g})".format(phi, precision))
308
        self.preMICH.apply_tuning(_kat,phi, add=True)
309

310
        vprint(verbose, "   scanning PRCL (maximising power)")
311
312
313
314
315
        kat = _kat.deepcopy()
        make_transparent(kat,["SRM"])
        phi, precision = self.scan_to_precision(kat, self.prePRCL, pretune_precision)
        phi=round(phi/pretune_precision)*pretune_precision
        phi=round_to_n(phi,5)
316
        vprint(verbose, "   found max/min at: {} (precision = {:2g})".format(phi, precision))
317
        self.prePRCL.apply_tuning(_kat,phi)
318

319
        vprint(verbose, "   scanning SRCL (maximising carrier power, then adding 90 deg)")
320
321
322
323
        kat = _kat.deepcopy()
        phi, precision = self.scan_to_precision(kat, self.preSRCL, pretune_precision, phi=0)
        phi=round(phi/pretune_precision)*pretune_precision
        phi=round_to_n(phi,4)-90.0
324
        vprint(verbose, "   found max/min at: {} (precision = {:2g})".format(phi, precision))
325
        self.preSRCL.apply_tuning(_kat,phi)
326
        print("   ... done")
327
        
328
329
330
331
332
333
334
    def scan_to_precision(self, kat, DOF, pretune_precision, minmax="max", phi=0.0, precision=60.0):
        while precision>pretune_precision*DOF.scale:
            out = scan_DOF(kat, DOF, xlimits = [phi-1.5*precision, phi+1.5*precision])
            phi, precision = find_peak(out, DOF.port.portName, minmax=minmax)
            #print("** phi= {}".format(phi))
        return phi, precision

335
336
337
338
339
340
341
342
343
344
345
346
347
348
    def pretune_status(self, _kat):
        kat = _kat.deepcopy()
        kat.verbose = False
        kat.noxaxis = True
        pretune_DOFs = [self.preARMX, self.preARMY, self.prePRCL, self.preMICH, self.preSRCL]
        _detStr=""
        for p in pretune_DOFs:
            _sigStr = p.port.signal(kat)
            _detStr = "\n".join([_detStr, _sigStr])
        kat.parseCommands(_detStr)
        out = kat.run()
        Pin = float(kat.L0.P)
    
        tunings = self.get_tunings(kat)
349
350
351
352
353
        _maxtemStr = "{:3}".format(tunings["maxtem"])
        if tunings["maxtem"] == -1:
            _maxtemStr="off"
        print(" .--------------------------------------------------.")
        print(" | pretuned for maxtem = {}, phase = {:2}            |".format(_maxtemStr, int(kat.phase)))
354
355
        keys_t = list(tunings.keys())
        keys_t.remove("maxtem")
356
357
358
        print(" .--------------------------------------------------.")
        print(" | port   power[W] pow. ratio | optics   tunings    |")
        print(" +----------------------------|---------------------+")
359
360
        idx_p = 0
        idx_t = 0
361
        while (idx_p < len(pretune_DOFs) or idx_t < len(keys_t)):
362
363
            if idx_p < len(pretune_DOFs):
                p = pretune_DOFs[idx_p]
364
                print(" | {:5}: {:9.4g} {:9.4g} |".format(p.name, float(out[p.port.name]), float(out[p.port.name])/Pin),end="")
365
366
367
368
369
                idx_p +=1
            else:
                print(" |                            |", end="")
            if idx_t < len(keys_t):
                t=keys_t[idx_t]
370
                print(" {:5}: {:9.3g}    |".format(t, float(self.tunings[t])))
371
372
                idx_t +=1
            else:
373
374
375
376
                print("                     |")
        print(" `--------------------------------------------------'")

    # probably extra and can be removed
377
378
379
380
381
382
383
384
    def power_ratios(self, _kat):
        kat = _kat.deepcopy()
        kat.verbose = False
        kat.noxaxis = True
    
        ports = [self.POW_X, self.POW_Y, self.AS_DC, self.POW_BS]
        _detStr = ""
        for p in ports:
385
            _sigStr = p.signal(kat)
386
387
388
            _detStr = "\n".join([_detStr, _sigStr])
        kat.parseCommands(_detStr)
        out = kat.run()
389
        
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
        Pin = float(kat.L0.P)
    
        print("-- power ratios (Pin = {0:.3g} W)".format(Pin))
        for p in ports:
            print(" {0:6} = {1:8.3g} W ({0:6}/Pin = {2:8.2g})" .format(p.name, float(out[p.name]), float(out[p.name])/Pin))
        return

    def plot_pretuning_powers(self, _kat, xlimits=[-10,10]):
        kat = _kat.deepcopy()
        kat.verbose = False
        kat.noxaxis = True
        dofs = [self.preARMX, self.preARMY, self.preMICH, self.prePRCL]
        idx=1
        fig = plt.figure()
        for d in dofs:
            ax = fig.add_subplot(2,2,idx)
            idx+=1
            out = scan_DOF(kat, d, xlimits = np.multiply(d.scale, xlimits), relative = True)
408
            ax.semilogy(out.x,out[d.signal_name(kat)])
409
410
            ax.set_xlim([np.min(out.x), np.max(out.x)])
            ax.set_xlabel("phi [deg] {}".format(d.optics[0]))
411
            ax.set_ylabel('{} [W] '.format(d.signal_name(kat)))
412
            ax.grid(True)
413
414
415
        plt.tight_layout()
        plt.show(block=0)

416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
    def _strToDOFs(self, DOFs):
        dofs = []
        
        for _ in DOFs:
            if isinstance(_, six.string_types):
                if _ in self.__DOFs:
                    dofs.append(self.__DOFs[_])
                else:
                    raise pkex.BasePyKatException("Could not find DOF called `%s`. Possible DOF options: %s" % (_, str(list(self.__DOFs.keys()))))
            else:
                raise pkex.BasePyKatException("'%s' not possible DOF options: %s" % (_, str(list(self.__DOFs.keys()))))
        
        return dofs

    def plot_error_signals(self, _kat, xlimits=[-1,1], DOFs=None, plotDOFs=None,
                                replaceDOFSignals=False, block=True, fig=None, legend=None):
        """
        Displays error signals for a given kat file. Can also be used to plot multiple
        DOF's error signals against each other for visualising any cross coupling.
        
        _kat: LIGO based kat object.
        xlimits: Range of DOF to plot in degrees
        DOFs: list, DOF names to compute. Default: DARM, CARM, PRCL, SRCL, MICH
        plotDOFs: list, DOF names to plot against each DOF. If None the same DOF as in DOFs is plotted.
        block: Boolean, for plot blocking terminal or not if being shown
        replaceDOFSignals: Bool, replaces already present signals for any DOF if already defined in kat. Regardless of this value, it will add default signals if none found.
        fig: figure, uses predefined figure, when defined it won't be shown automatically
        legend: string, if no plotDOFs is defined this legend is shown
        
        Example:
            import pykat
            from pykat.gw_detectors import ifo

            ligo = ifo.aLIGO()
            
            # Plot default
            ligo.plot_error_signals(ligo.kat, block=True)
            # plot MICH and CARM against themselves
            ligo.plot_error_signals(ligo.kat, DOFs=["MICH", "CARM"], block=True)
            # plot DARM and CARM against MICH
            ligo.plot_error_signals(ligo.kat, DOFs=["MICH"], plotDOFs=["DARM", "CARM"], block=True)
        """
        
459
460
461
        kat = _kat.deepcopy()
        kat.verbose = False
        kat.noxaxis = True
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
        
        if DOFs is None:
            dofs = [self.DARM, self.CARM, self.PRCL, self.SRCL, self.MICH]
        else:
            dofs = self._strToDOFs(DOFs)
        
        # add in signals for those DOF to plot
        for _ in dofs:
            if not (not replaceDOFSignals and hasattr(kat, _.signal_name(kat))):
                kat.parseCommands(_.signal(kat))
                
        toShow = None
        
        if plotDOFs is not None:
            toShow = self._strToDOFs(plotDOFs)
        
            # Check if other DOF signals we need to include for plotting
            for _ in toShow:
                if not (not replaceDOFSignals and hasattr(kat, _.signal_name(kat))):
                    kat.parseCommands(_.signal(kat))
                    
        if fig is not None:
            _fig = fig
        else:
            _fig = plt.figure()
        
        nrows = 2
        ncols = 3
        
        if DOFs is not None:
            n = len(DOFs)
            
            if n < 3:
                nrows = 1
                ncols = n
        
        for d, idx in zip(dofs, range(1, len(dofs)+1)):
            ax = _fig.add_subplot(nrows, ncols, idx)
            
            scan_cmd = scan_optics_string(d.optics, d.factors, "scan", linlog="lin",
                                            xlimits=np.multiply(d.scale, xlimits), steps=200,
                                            axis=1, relative=True)
            kat.parseCommands(scan_cmd)
            out = kat.run()
            
            if toShow is None:
                ax.plot(out.x, out[d.signal_name(kat)], label=legend)
            else:
                for _ in toShow:
Daniel Brown's avatar
Daniel Brown committed
511
512
513
514
                    if legend is None:
                        legend = _.name
                        
                    ax.plot(out.x, out[_.signal_name(kat)], label=legend)
515
                
516
517
            ax.set_xlim([np.min(out.x), np.max(out.x)])
            ax.set_xlabel("{} [deg]".format(d.name))
Daniel Brown's avatar
Daniel Brown committed
518
519
520
521
            
            if plotDOFs is None:
                ax.set_ylabel('{} [W] '.format(d.port.name))
            else:
Daniel Brown's avatar
Daniel Brown committed
522
                ax.set_ylabel('Error signal [W]')
523
524
525
526
527
528
            
            ax.grid(True)
        
        if toShow is not None or legend is not None:
            plt.legend(loc=0)
           
529
        plt.tight_layout()
530
531
532
        
        if fig is None:
            plt.show(block=block)
533

534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
    def set_DC_offset(self, _kat, DCoffset=None, verbose=False):
        if DCoffset:
            self.DCoffset=DCoffset
            print("-- applying user-defined DC offset:")
            pretuning = self.get_tunings(_kat)
            pretuning["ETMY"] += self.DCoffset
            pretuning["ETMX"] -= self.DCoffset
            self.set_tunings(_kat, pretuning)        
            kat = _kat.deepcopy()
            sigStr = self.AS_DC.signal(kat)
            signame = self.AS_DC.signal_name(kat)
            kat.parseCommands(sigStr)
            kat.noxaxis=True
            out = kat.run()
            self.DCoffsetW=float(out[signame])
        else:
            # Finding light power in AS port (mostly due to RF sidebands now
            kat = _kat.deepcopy()
            sigStr = self.AS_DC.signal(kat)
            signame = self.AS_DC.signal_name(kat)
            kat.parseCommands(sigStr)
            kat.noxaxis=True
            out = kat.run()
            print("-- adjusting DCoffset based on light in dark port:")
            waste_light = round(float(out[signame]),1)
            print("   waste light in AS port of {:2} W".format(waste_light))
            kat_lock = _kat.deepcopy()
            self.find_DC_offset(kat_lock, 2*waste_light)
            pretuning = self.get_tunings(kat_lock)
            pretuning["ETMY"] += self.DC_offset
            pretuning["ETMX"] -= self.DC_offset
            self.set_tunings(_kat, pretuning)
        self.DCoffset_meter = self.DCoffset / 360.0 * _kat.lambda0 
        vprint(verbose, "   DCoffset = {:6.4} deg ({:6.4}m)".format(self.DCoffset, self.DCoffset_meter))
        vprint(verbose, "   at dark port power: {:6.4}W".format(self.DCoffsetW))
    

    def find_DC_offset(self, _kat, AS_power, precision=1e-4, verbose=False):
572
573
574
575
        """
        Returns the DC offset of DARM that corrponds to the
        specified power in the AS power.
        """
576
        vprint(verbose, "   finding DC offset for AS power of {:3g} W".format(AS_power))
577
578
579
        kat = _kat.deepcopy()
        kat.verbose = False
        kat.noxaxis = True
580
        _sigStr = self.AS_DC.signal(kat)
581
582
583
584
585
586
587
588
589
590
591
        kat.parseCommands(_sigStr)
        Xphi = float(kat.ETMX.phi)
        Yphi = float(kat.ETMY.phi)

        def powerDiff(phi, kat, Xphi, Yphi, AS_power):
            kat.ETMY.phi = Yphi + phi
            kat.ETMX.phi = Xphi - phi
            out=kat.run()
            #print(out[self.AS_DC.name]-AS_power)
            return np.abs(out[self.AS_DC.name]-AS_power)

592
593
594
595
        vprint(verbose, "   starting peak search...")
        out=fmin(powerDiff,0,xtol=precision,ftol=1e-3,args=(kat, Xphi, Yphi, AS_power), disp=verbose)
        vprint(verbose, "   ... done")
        vprint(verbose, "   DC offset for AS_DC={} W is: {}".format(AS_power, out[0]))
596
597
598
599
        self.DCoffset = round(out[0],6)
        self.DCoffsetW = AS_power
        return self.DCoffset

600
    def generate_tuning_block(self, kat):
601
        code1 = """###########################################################################
602
603
604
605
606
const phi_ITMX {:.8}
const phi_ITMY {:.8}
const phi_ETMX {:.8}
const phi_ETMY {:.8}""".format(float(kat.ITMX.phi), float(kat.ITMY.phi), float(kat.ETMX.phi),float(kat.ETMY.phi))

607
        code2 = """
608
609
const phi_BS   {:.8}
const phi_PRM  {:.8}
610
611
const phi_SRM  {:.8}
###########################################################################""".format(float(kat.BS.phi), float(kat.PRM.phi), float(kat.SRM.phi))
612

613
        return "".join([code1, code2])
614

615
    def generate_errsig_block(self, kat, noplot=False):
616
        code1 = "###########################################################################\n" 
617
618
619
620
621
        sigDARM = self.DARM.signal(kat)
        sigCARM = self.CARM.signal(kat)
        sigPRCL = self.PRCL.signal(kat)
        sigMICH = self.MICH.signal(kat)
        sigSRCL = self.SRCL.signal(kat)
622
        code2 = "\n".join([sigDARM, sigCARM, sigPRCL, sigMICH, sigSRCL])
623

624
        code3= ""
625
626
627
628
629
630
        if noplot:
            nameDARM = self.DARM.signal_name(kat)
            nameCARM = self.CARM.signal_name(kat)
            namePRCL = self.PRCL.signal_name(kat)
            nameMICH = self.MICH.signal_name(kat)
            nameSRCL = self.SRCL.signal_name(kat)
631
632
633
634
635
636
            code3 = """
noplot {}
noplot {}
noplot {}
noplot {}
noplot {}""".format(nameDARM, nameCARM, namePRCL, nameMICH, nameSRCL)
637
        code4="\n###########################################################################"
638
639
        return "".join([code1, code2, code3, code4])

640
            
641
        
642
    def generate_locks(self, _kat, gainsAdjustment = [0.5, 0.005, 1.0, 0.5, 0.025], gains=None, accuracies=None, verbose=True):
643
        """
644
645
646
        gainsAdjustment: factors to apply to loop gains computed from optical gains
        gains: override loop gain [W per deg]
        accuracies: overwrite error signal threshold [W]
647
648
        """
        kat = _kat.deepcopy()
649
650
651
652
653
654
655
        # optical gains in W/rad
        ogDARM = optical_gain(kat, self.DARM, self.DARM)
        ogCARM = optical_gain(kat, self.CARM, self.CARM)
        ogPRCL = optical_gain(kat, self.PRCL, self.PRCL)
        ogMICH = optical_gain(kat, self.MICH, self.MICH)
        ogSRCL = optical_gain(kat, self.SRCL, self.SRCL)

656
        if gains == None:            
657
658
            # manually tuning relative gains
            factor = -1.0 * 180 / math.pi # convert from rad/W to -1 * deg/W
659
660
661
662
663
            gainDARM = round_to_n(gainsAdjustment[0] * factor / ogDARM, 2) # manually tuned
            gainCARM = round_to_n(gainsAdjustment[1] * factor / ogCARM, 2) # factor 0.005 for better gain hirarchy with DARM
            gainPRCL = round_to_n(gainsAdjustment[2] * factor / ogPRCL, 2) # manually tuned
            gainMICH = round_to_n(gainsAdjustment[3] * factor / ogMICH, 2) # manually tuned
            gainSRCL = round_to_n(gainsAdjustment[4] * factor / ogSRCL, 2) # gain hirarchy with MICH
664
            gains = [ gainDARM, gainCARM, gainPRCL, gainMICH, gainSRCL]
665
        self.lockGains = copy.deepcopy(gains)
666
        
667
668
669
670
671
672
673
        # rms: loop accuracies in meters (manually tuned for the loops to work
        # with the default file)
        # to compute accuracies from rms, we convert
        # rms to radians as rms_rad = rms * 2 pi/lambda
        # and then multiply by the optical gain.
        rms = [1e-13, 1e-13, 1e-12, 1e-11, 50e-11] # default accuracies in meters
        factor = 2.0 * math.pi / kat.lambda0 # convert from m to radians
674
        if accuracies == None:
675
676
677
678
679
            accDARM = round_to_n(np.abs(factor * rms[0] * ogDARM),2) 
            accCARM = round_to_n(np.abs(factor * rms[1] * ogCARM),2) 
            accPRCL = round_to_n(np.abs(factor * rms[2] * ogPRCL),2) 
            accMICH = round_to_n(np.abs(factor * rms[3] * ogMICH),2)
            accSRCL = round_to_n(np.abs(factor * rms[4] * ogSRCL),2) 
680
            accuracies = [accDARM, accCARM, accPRCL, accMICH, accSRCL]
681
        self.lockAccuracies = copy.deepcopy(accuracies)
682

683
684
685
686
687
688

        nameDARM = self.DARM.signal_name(kat)
        nameCARM = self.CARM.signal_name(kat)
        namePRCL = self.PRCL.signal_name(kat)
        nameMICH = self.MICH.signal_name(kat)
        nameSRCL = self.SRCL.signal_name(kat)
689
        self.lockNames = [nameDARM, nameCARM, namePRCL, nameMICH, nameSRCL]
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
        factor1 = 2.0 * math.pi / 360.0 
        factor2 = 2.0 * math.pi / kat.lambda0 
        factor3 = 360.0  / kat.lambda0
        factor4 = -1.0 * 180 / math.pi 

        if verbose:
            print(" .--------------------------------------------------.")
            print(" | Parameters for locks:                            |")
            print(" +--------------------------------------------------+")
            print(" | -- optical gains [W/rad], [W/deg] and [W/m]:     |")
            print(" | DARM: {:12.5}, {:12.5}, {:12.5}   |".format(ogDARM, ogDARM*factor1, ogDARM*factor2))
            print(" | CARM: {:12.5}, {:12.5}, {:12.5}   |".format(ogCARM, ogCARM*factor1, ogCARM*factor2))
            print(" | PRCL: {:12.5}, {:12.5}, {:12.5}   |".format(ogPRCL, ogPRCL*factor1, ogPRCL*factor2))
            print(" | MICH: {:12.5}, {:12.5}, {:12.5}   |".format(ogMICH, ogMICH*factor1, ogMICH*factor2))
            print(" | SRCL: {:12.5}, {:12.5}, {:12.5}   |".format(ogSRCL, ogSRCL*factor1, ogSRCL*factor2))
            print(" +--------------------------------------------------+")
            print(" | -- defult loop accuracies [deg], [m] and [W]:    |")
            print(" | DARM: {:12.6}, {:12.6}, {:12.6}   |".format(factor3*rms[0], rms[0], np.abs(rms[0]*ogDARM*factor2)))
            print(" | CARM: {:12.6}, {:12.6}, {:12.6}   |".format(factor3*rms[1], rms[1], np.abs(rms[1]*ogCARM*factor2)))
            print(" | PRCL: {:12.6}, {:12.6}, {:12.6}   |".format(factor3*rms[2], rms[2], np.abs(rms[2]*ogPRCL*factor2)))
            print(" | MICH: {:12.6}, {:12.6}, {:12.6}   |".format(factor3*rms[3], rms[3], np.abs(rms[3]*ogMICH*factor2)))
            print(" | SRCL: {:12.6}, {:12.6}, {:12.6}   |".format(factor3*rms[4], rms[4], np.abs(rms[4]*ogSRCL*factor2)))
            print(" +--------------------------------------------------+")
            print(" | -- extra gain factors (factor * 1/optical_gain): |")
714
715
716
717
718
            print(" | DARM: {:5.4} * {:12.6} = {:12.6}        |".format(gainsAdjustment[0],factor4/ogDARM, gainsAdjustment[0]*factor4/ogDARM))
            print(" | CARM: {:5.4} * {:12.6} = {:12.6}        |".format(gainsAdjustment[1],factor4/ogCARM, gainsAdjustment[1]*factor4/ogCARM))
            print(" | PRCL: {:5.4} * {:12.6} = {:12.6}        |".format(gainsAdjustment[2],factor4/ogPRCL, gainsAdjustment[2]*factor4/ogPRCL))
            print(" | MICH: {:5.4} * {:12.6} = {:12.6}        |".format(gainsAdjustment[3],factor4/ogMICH, gainsAdjustment[3]*factor4/ogMICH))
            print(" | SRCL: {:5.4} * {:12.6} = {:12.6}        |".format(gainsAdjustment[4],factor4/ogSRCL, gainsAdjustment[4]*factor4/ogSRCL))
719
720
721
722
723
724
725
            print(" `--------------------------------------------------'")


        
    def generate_lock_block(self, kat, verbose=False):
        if self.lockNames == None or self.lockAccuracies == None or self.lockGains == None:
            raise pkex.BasePyKatException("run gerate_locks before generate_lock_block")            
726
        code1 = """###########################################################################
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
set AS_f2_I_re {} re
set CARM_err {} re
set PRCL_err {} re
set MICH_err {} re
set SRCL_err {} re
func DARM_err = $AS_f2_I_re - {}
""".format(self.lockNames[0], self.lockNames[1], self.lockNames[2], self.lockNames[3], self.lockNames[4], self.DCoffsetW)
       
        code2 = """
lock DARM_lock $DARM_err {:8.2} {:8.2}
lock CARM_lock $CARM_err {:8.2g} {:8.2g} 
lock PRCL_lock $PRCL_err {:8.2g} {:8.2g}
lock MICH_lock $MICH_err {:8.2g} {:8.2g} 
lock SRCL_lock $SRCL_err {:8.2g} {:8.2g}""".format(self.lockGains[0], self.lockAccuracies[0], self.lockGains[1], self.lockAccuracies[1], self.lockGains[2], self.lockAccuracies[2], self.lockGains[3], self.lockAccuracies[3], self.lockGains[4], self.lockAccuracies[4])
        
        code3 = """
noplot ITMY_lock
func ITMY_lock = (-1.0) * $MICH_lock 
func ETMX_lock = $CARM_lock + $MICH_lock + $DARM_lock
func ETMY_lock = $CARM_lock - $MICH_lock - $DARM_lock

put* PRM     phi     $PRCL_lock
put* ITMX    phi     $MICH_lock
put* ITMY    phi     $ITMY_lock
put* ETMX    phi     $ETMX_lock
put* ETMY    phi     $ETMY_lock
put* SRM     phi     $SRCL_lock

noplot PRCL_lock
noplot SRCL_lock
noplot MICH_lock
noplot DARM_lock
noplot CARM_lock
noplot ETMX_lock
noplot ETMY_lock
###########################################################################
763
"""
764
765
766
        if verbose:
            print(" .--------------------------------------------------.")
            print(" | Lock commands used:                              |")
767
768
769
770
            print(" +--------------------------------------------------+")
            for l in code2.splitlines():
                print (" | {:49}|".format(l))
            print(" `--------------------------------------------------'")
771
772
        return "".join([code1, code2, code3])

773
774
775
776
# ---------------------------------------------------------------------------------
# ---------------------------------------------------------------------------------
# ---------------------------------------------------------------------------------
# ---------------------------------------------------------------------------------
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
class DOF(object):
    """
    Defining a degree of freedom for the interferometer, includes the
    objects and how to move them, and the default output port to read
    out the DOF signal.
    """
    def __init__(self, _DOFName, _port, _quad, _optics, _factors, _scale, sigtype="z"):
        self.name = _DOFName
        self.port = _port
        self.quad = _quad
        self.sigtype = sigtype
        self.optics=make_list_copy(_optics)
        self.factors=make_list_copy(_factors)
        # scaling factor, to compensate for lower sensitivity compared
        # to DARM (in tuning plots for example)
        # Thus DARM has a scale of 1, all other DOFs a scale >1
        self.scale = _scale

795
796
797
798
799
800
801
    def apply_tuning(self, kat, phi, add=False):
        for idx, o in enumerate(self.optics):
            if add:
                kat.components[o].phi += phi * self.factors[idx]
            else:
                kat.components[o].phi = phi * self.factors[idx]
            
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
    def signal(self, kat):
        return self.port.signal(kat, self.quad, sigtype=self.sigtype)
    def signal_name(self, kat):
        return self.port.signal_name(kat, self.quad, sigtype=self.sigtype)

    def transfer(self, kat, fsig, phase2=None):
        return self.port.transfer(kat, self.quad, fsig=fsig, phase2=phase2, sigtype=self.sigtype)
    def transfer_name(self, kat):
        return self.port.transfer_name(kat, self.quad)

    def fsig(self, _fsigName, fsig=1.0):
        _fsigStr= ""
        for idx, o in enumerate(self.optics):
            phase = 0.0
            if self.factors[idx] == -1:
                phase = 180.0
            _fsigStr = "\n".join([_fsigStr, "fsig {} {} {} {} ".format(_fsigName, o, fsig, phase)])
        return _fsigStr

821
822
823
# ---------------------------------------------------------------------------------
# ---------------------------------------------------------------------------------
# ---------------------------------------------------------------------------------
824
825
826
827
828
class port(object):
    """
    Defining an output port for the interferometer, can be either a
    pd or a pd1 detector (for error signal generation).
    """
829
    def __init__(self, _portName, _nodeNames, f=None, phase=None):
830
        self.portName = _portName
831
        self.nodeNames = make_list_copy(_nodeNames)
832
        self.f=f            # demodulation frequency, float
833
        self.phase = phase  # demodulation frequency for I quadrature, float
834
        self.name = self.portName            
835
836
837
838
839
840
841
842
843
844
845
846
847
    
    def check_nodeName(self, kat):
        self.nodeName = None
        for node in self.nodeNames:
            _node = node
            if _node[-1] == "*":
                _node = _node[:-1]
            if _node in kat.nodes:
                self.nodeName=node
                break
        if self.nodeName==None:
            raise pkex.BasePyKatException("port {}: cannot find any of these nodes: '{}'".format(self.name,self.nodeNames))

848
849
850
851
852
    def amplitude_name(self, kat, f, n=None, m=None, sigtype="z"):
        name = self.name + "_ad"
        return name
    
    def amplitude(self, kat, f, n=None, m=None, sigtype="z"):
853
        self.check_nodeName(kat)
854
855
856
857
858
859
860
        name = self.amplitude_name(kat, f, n=n, m=m, sigtype=sigtype)
        if n==None and m==None:
            return "ad {} {} {}".format(name, f, self.nodeName)
        else:
            return "ad {} {} {} {} {}".format(name, f, n, m, self.nodeName)
    
    def signal_name(self, kat, quad="I", sigtype="z"):
861
        name = self.name
862
863
864
865
866
867
868
        if self.f!=None:
            name = self.name+"_"+quad
        return name
    
    def signal(self, kat, quad="I", sigtype="z"):
        self.check_nodeName(kat)
        name = self.signal_name(kat, quad=quad, sigtype=sigtype)
869
870
871
        if sigtype != "z":
                raise pkex.BasePyKatException("alignment signals are not implemented yet")            
        if self.f==None:
872
            return "pd {} {}".format(name, self.nodeName)
873
874
875
        else:
            if quad !="I" and quad != "Q":
                raise pkex.BasePyKatException("quadrature must be 'I' or 'Q'")            
876
            phase = self.IQ_phase(quad, self.phase)
877
            return "pd1 {} {} {} {}".format(name, self.f, phase, self.nodeName)
878
        
879
    def IQ_phase(self, quad, phase):
880
881
882
883
884
885
        if quad== "Q":
            phase = phase + 90.0
            if phase >=360.0 :
                phase -= 360.0
        return phase
        
886
    def transfer_name(self, kat, quad="I"):
887
        name = self.name
888
889
890
891
892
        if self.f!=None:
            name = self.name+"_"+quad
        return name
        
    def transfer(self, kat, quad="I", fsig=1.0, phase2=None, sigtype="z"):
893
        self.check_nodeName(kat)
894
        name = self.transfer_name(kat, quad=quad)
895
896
897
        if sigtype!="z":
                raise pkex.BasePyKatException("alignment signals are not implemented yet")            
        if self.f==None:
898
            if phase2 == None:
899
                return "pd1 {} {} {}".format(name, fsig, self.nodeName)
900
            else:
901
                return "pd1 {} {} {} {}".format(name, fsig, phase2, self.nodeName)
902
903
904
        else:
            if quad !="I" and quad != "Q":
                raise pkex.BasePyKatException("quadrature must be 'I' or 'Q'")            
905
            phase = self.IQ_phase(quad, self.phase)
906
            if phase2 == None:
907
                return "pd2 {} {} {} {} {}".format(name, self.f , phase, fsig, self.nodeName)
908
            else:
909
                return "pd2 {} {} {} {} {} {}".format(name, self.f, phase, fsig, phase2, self.nodeName)
910
911
912
913
914
                                
        
def scan_optics_string(_optics, _factors, _varName, linlog="lin", xlimits=[-100, 100], steps=200, axis=1,relative=False):
    optics=make_list_copy(_optics)
    factors=make_list_copy(_factors)
915
    
916
917
    if len(optics) != len(factors):
        raise pkex.BasePyKatException("you must provide a factor for each optics")
918

919
920
    if linlog not in ["lin", "log"]: 
        raise pkex.BasePyKatException("linlog must be 'lin' or 'log'")
921
        
922
    _tuneStr  = "var {} 0\n".format(_varName)
923
    
924
925
926
927
928
929
    if axis==1:
        _tuneStr += "xaxis {} phi {} {} {} {}".format(_varName, linlog, xlimits[0], xlimits[1], steps)
    elif (axis==2 or axis==3): 
        _tuneStr += "x{}axis {} phi {} {} {} {}".format(axis, _varName, linlog, xlimits[0], xlimits[1], steps)
    else:
        raise pkex.BasePyKatException("axis must be 1, 2 or 3")
930
        
931
    _putStr = ""
932
    
933
934
935
936
937
938
939
940
941
942
943
944
945
946
    for idx, o in enumerate(optics):
        if factors[idx] == 1:
            _xStr="$x"
        elif factors[idx] == -1:
            _xStr="$mx"
        else:
            raise pkex.BasePyKatException("optics factors must be 1 or -1")
        if (relative):
            _putCmd = "put*"
        else:
            _putCmd = "put"
                
        _putStr = "\n".join([_putStr, "{} {} phi {}{}".format(_putCmd, o, _xStr, axis)])            
        _tuneStr += _putStr
947
        
948
949
    return _tuneStr
        
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
def make_transparent(kat, _components):
    """
    Function to make certain mirror or beamsplitter objects transparent
    Usage:
    make_transparent(kat, "PRM")
    or
    make_transparent(kat, ["PRM", "SRM"])
    The function modifies the kat object. Use something like 
    kat = copy.deepcopy(basekat)
    before if you want to preseve the original kat object
    """
    components=make_list_copy(_components)
    for o in kat.components.values():
        if o.name in components:
            if isinstance(o, pykat.components.AbstractMirrorComponent): 
                o.setRTL(0.0,1.0,0.0)
                components = [c for c in components if c != o.name]
            else:
                raise pkex.BasePyKatException("{} is not a mirror or beamsplitter".format(o))
969
            #print('  made {} transparent'.format(o))
970
971
972
    if len(components) != 0:
        raise pkex.BasePyKatException("Cannot find component {}".format(components))
    return kat
Andreas Freise's avatar
Andreas Freise committed
973

974
def reconnect_nodes(kat, component1, idx1, node_name, verbose=False):
Andreas Freise's avatar
Andreas Freise committed
975
976
977
978
979
980
981
982
    c_string = component1.getFinesseText()
    c = c_string[0].split()
    new_string = " ".join(c[:-2])
    nodes = {}
    nodes[0] = c[-2]
    nodes[1] = c[-1]
    nodes[idx1]=node_name
    new_string = new_string + " " + nodes[0] + " " + nodes[1]
983
    vprint(verbose, "   new string ='{}'".format(new_string))
Andreas Freise's avatar
Andreas Freise committed
984
985
    kat.parseCommands(new_string)

986
def remove_commands(kat, _commands, verbose=False):
Andreas Freise's avatar
Andreas Freise committed
987
988
989
990
991
992
    commands=make_list_copy(_commands)
    # removing commands
    for o in kat.commands.values():
        if o.name in commands:
            o.remove()
            commands = [c for c in commands if c != o.name]
993
            vprint(verbose, '   {} removed'.format(o))
Andreas Freise's avatar
Andreas Freise committed
994
995
996
997
    if len(commands) != 0:
        raise pkex.BasePyKatException("Cannot find command(s) {}".format(commands))
    return kat
    
998
def remove_components(kat, _components, component_in=None, component_out=None, verbose=False):
Andreas Freise's avatar
Andreas Freise committed
999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
    components=make_list_copy(_components)
    if  kat.components[components[-1]].nodes[1]:
        node_in  = kat.components[components[-1]].nodes[1].name
    else:
        node_in = None
    node_out = kat.components[components[0]].nodes[0].name
    # removing components
    for o in kat.components.values():
        if o.name in components:
            o.remove()
            components = [c for c in components if c != o.name]
1010
            vprint(verbose, '  {} removed'.format(o))
Andreas Freise's avatar
Andreas Freise committed
1011
1012
1013
1014
1015
1016
1017
1018
    if len(components) != 0:
        raise pkex.BasePyKatException("Cannot find component(s) {}".format(components))
    # reconnecting nodes if requested
    if component_in:
        reconnect_nodes(kat, kat.components[component_in],1, node_in)
    if component_out:
        reconnect_nodes(kat, kat.components[component_out],0, node_out)
    return kat
1019
1020
        
def BS_optical_path(thickness, n=nsilica, angle=45.0):
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
    """
    Compute optical path length in BS substrate, default
    parameters assume angle of incidence of 45 deg and fused
    silica substrate.
    thickness: substrate thickness [m]
    n: index of refraction of substrate
    angle: angle of incidence (in vacuum) [deg]
    """
    
    angle_subst = math.asin(math.sin(math.radians(angle))/n)
    L = thickness / math.cos(angle_subst) 
    
    return math.degrees(angle_subst), L
1034

1035
1036
def scan_DOF(_kat, DOF, xlimits=[-100, 100], steps=200, relative=False): 
    kat = _kat.deepcopy()
1037
1038
1039
1040
    
    scan_string = scan_optics_string(DOF.optics, DOF.factors, "scan", linlog="lin",
                                     xlimits=xlimits, steps=steps, axis=1, relative=relative)
    
1041
    kat.parseCommands(scan_string)
1042
1043
1044
    kat.parseCommands(DOF.signal(kat))
    
    return kat.run()
1045

1046
def scan_optics(_kat, _optics, _factors, xlimits=[-100, 100], steps=200,relative=False): 
Andreas Freise's avatar
Andreas Freise committed
1047
    """
1048
    Scans one or more optics (by changing its tuning).
Andreas Freise's avatar
Andreas Freise committed
1049
1050
1051
1052
1053
1054
    Parameters:
    optics: list of names of components to be tuned
    factors: list of scaling factors for the tuning for each optics, first element must be 1.0
    xlimits: limits of scan, defaults to [-100, 100]
    steps: number of steps to use in scan, default is 200
    Usage:
1055
1056
    scan_optis(kat, "PRM", 1)
    scan_optis(kat, ["ETMX", "ETMY", [1, -1])
Andreas Freise's avatar
Andreas Freise committed
1057
1058
    """
    kat = _kat.deepcopy()
1059
    
Andreas Freise's avatar
Andreas Freise committed
1060
1061
    optics=make_list_copy(_optics)
    factors=make_list_copy(_factors)
1062
1063
1064
1065
    
    scan_string = scan_optics_string(optics, factors, "scan", linlog="lin", xlimits=xlimits,
                                      steps=steps, axis=1,relative=relative)
    
1066
    kat.parseCommands(scan_string)
1067
1068
    
    return kat.run()
1069

1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
def optical_gain(_kat, DOF_sig, DOF_det, f=10.0):
    kat = _kat.deepcopy()
    _fsigStr = DOF_sig.fsig("sig1", fsig=f)
    _detStr = DOF_det.transfer(kat, fsig=f)
    _detName = DOF_det.transfer_name(kat)
    kat.parseCommands(_fsigStr)
    kat.parseCommands(_detStr)
    kat.noxaxis = True
    kat.parseCommands("yaxis lin abs:deg")
    out = kat.run()
1080
    return float(np.real(out[_detName]))
1081

1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
def find_peak(out, detector, minmax='max', debug=False): 
    """
    Expects an output of a kat scan and find the max/min output on
    a sepcific detector. Useful for pre-tuning cavities or interferometers.
    Returns the tuning of the maximum (or minimum) and the precision
    Parameters:
    detector: name of detetor to use
    minmax, string to indicate maximum or minum detection, default 'max'
    Usage:
    find_peak(out, "pdout")
    find_peak(out, "pdout", minmax='min')
    """
Andreas Freise's avatar
Andreas Freise committed
1094
    stepsize = out.x[1]-out.x[0]
1095
1096
    if debug:
        print("  stepsize (precision) of scan: {0:g}".format(stepsize))
Andreas Freise's avatar
Andreas Freise committed
1097
1098
1099

    _max, _min = peak.peakdetect( out[detector],out.x, 1)
    
1100
    if debug:
1101
1102
        plt.figure()
        plt.plot(out.x,out[detector])
Andreas Freise's avatar
Andreas Freise committed
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
        print("max: ")
        print(_max)
        print("min: ")
        print(_min)
        
    if minmax == 'max':
        X = [p[0] for p in _max]
        Y = [p[1] for p in _max]
        X_out = X[np.argmax(Y)]
        Y_out = np.max(Y)
    elif minmax == 'min':
        X = [p[0] for p in _min]
        Y = [p[1] for p in _min]
        X_out = X[np.argmin(Y)]
        Y_out = np.min(Y)
    else:
        raise pkex.BasePyKatException("maxmin must be 'max' or 'min'")
        
1121
    if debug:
Andreas Freise's avatar
Andreas Freise committed
1122
        plt.plot(X_out,Y_out,'o')
1123
        plt.xlabel('tuning [deg]')
Andreas Freise's avatar
Andreas Freise committed
1124
1125
1126
1127
        plt.ylabel('{0} output'.format(detector))
        plt.show()
    return X_out, stepsize

1128
1129
1130
1131
1132
1133
1134
1135
1136
1137
def make_list_copy(_l):
    """
    Utility function, takes a list of strings or single string
    and returns a copy of a list, e.g.
    "string" copy to ["string"]
    ["string1", "string2"] copy to ["string1", "string2"]
    """
    if not isinstance(_l, list):
        _l = [_l]
    return _l[:] # copy the list, just to be save
1138
1139
1140
1141
1142
1143

def round_to_n(x, n):
    if not x: return 0
    power = -int(math.floor(math.log10(abs(x)))) + (n - 1)
    factor = (10 ** power)
    return round(x * factor) / factor
1144
1145
1146
1147
1148


def vprint(verbose, printstr):
    if verbose:
        print(printstr)