diff --git a/bin/test_2d_plot.py b/bin/test_2d_plot.py index 82d95e7feb27e6a882a0c15e5b568eeff554c2d8..30786b86c270dfcb3f4c0233dfd5866a3fcef1b6 100644 --- a/bin/test_2d_plot.py +++ b/bin/test_2d_plot.py @@ -1,5 +1,4 @@ from pykat import finesse -from pykat.utilities.optics.gaussian_beams import gauss_param from pykat.detectors import * from pykat.components import * from pykat.commands import * @@ -26,7 +25,7 @@ x2axis b1 y lin -6 6 100 kat = finesse.kat() kat.parseCommands(code) -kat.s1.n1.q = gauss_param(w0=1e-3, z=0) +kat.s1.n1.q = pykat.beam_param(w0=1e-3, z=0) out = kat.run(printout=0,printerr=0) diff --git a/pykat/finesse.py b/pykat/finesse.py index 2d1338594806fe86cff1c858a1893f06e63057ba..e06f322b31c9eeb117aed32dfae2f6829aeb1fee 100644 --- a/pykat/finesse.py +++ b/pykat/finesse.py @@ -249,7 +249,7 @@ class katRun(object): class katRun2D(object): def __init__(self): - self.runtime + self.runtime = None self.startDateTime = datetime.datetime.now() self.x = None self.y = None @@ -1692,10 +1692,13 @@ class kat(object): if self.noxaxis == True: out.append("noxaxis\n") - if self.yaxis != None: + if self.yaxis is not None: out.append("yaxis {0}\n".format(self.yaxis)) - if self.printmatrix != None and self.printmatrix == True: + if self.retrace is not None: + out.append("retrace {0}\n".format(str(self.retrace))) + + if self.printmatrix is not None and self.printmatrix == True: out.append("printmatrix\n") if self.lambda0 != 1064e-9: diff --git a/pykat/maths/hermite.py b/pykat/maths/hermite.py index 84f729ccd554046db57bb96cf353a680d2e5d9bc..f300ad42f7746b24c88fdee764d38adba1a97a8e 100644 --- a/pykat/maths/hermite.py +++ b/pykat/maths/hermite.py @@ -10,36 +10,36 @@ def hermite(n, x): elif n == 3: return (8.0 * x * x * x - 12.0 * x) elif n == 4: - x_sq = x * x + x_sq = x * x return (16.0 * x_sq*x_sq - 48.0 * x_sq + 12.0) elif n == 5: - x_sq = x * x - x_cb = x * x_sq + x_sq = x * x + x_cb = x * x_sq return (32.0 * x_cb * x_sq - 160.0 * x_cb + 120.0 * x) elif n == 6: - x_sq = x * x - x_four = x_sq * x_sq + x_sq = x * x + x_four = x_sq * x_sq return (64.0 * x_four*x_sq - 480.0 * x_four + 720.0 * x_sq - 120.0) elif n == 7: x_sq = x*x - x_cb = x_sq*x - x_four = x_cb*x + x_cb = x_sq*x + x_four = x_cb*x return (128.0 * x_cb*x_four - 1344.0 * x_cb*x_sq + 3360.0 * x_cb - 1680.0 * x) elif n == 8: - x_sq = x*x - x_four = x_sq*x_sq - x_six = x_four*x_sq + x_sq = x*x + x_four = x_sq*x_sq + x_six = x_four*x_sq return (256.0 * x_four*x_four - 3584.0 * x_six + 13440.0 * x_four - 13440.0 * x_sq + 1680.0) elif n == 9: - x_sq = x*x - x_cb = x_sq*x - x_four = x_cb*x + x_sq = x*x + x_cb = x_sq*x + x_four = x_cb*x return (512.0 * x_cb*x_cb*x_cb - 9216.0 * x_four*x_cb + 48384.0 * x_cb*x_sq - 80640.0 * x_cb + 30240.0 * x) elif n == 10: - x_sq = x*x - x_cb = x_sq*x - x_four = x_cb*x - x_five = x_four * x + x_sq = x*x + x_cb = x_sq*x + x_four = x_cb*x + x_five = x_four * x return (1024.0 * x_five*x_five - 23040.0 * x_four*x_four + 161280.0 * x_cb*x_cb - 403200.0 * x_four + 302400.0 * x_sq - 30240.0) elif n == 11: return -665280 * x + 2217600 * x**3 - 1774080 * x**5 + 506880 * x**7 - 56320 * x**9 + 2048 * x**11 diff --git a/pykat/optics/maps.py b/pykat/optics/maps.py index 86b2e066947e308ce398ac66e54e98960cf37b36..cfea9374bb210d23c449aacc38b78a52509ca47d 100644 --- a/pykat/optics/maps.py +++ b/pykat/optics/maps.py @@ -303,7 +303,182 @@ class surfacemap(object): pylab.show() return fig + +class mergedmap: + """ + A merged map combines multiple surfaces map to form one. Such a map can be used + for computations of coupling coefficients but it cannot be written to a file to + be used with Finesse. For this you must output each map separately. + + """ + + def __init__(self, name, size, center, step_size, scaling): + + self.name = name + self.center = center + self.step_size = step_size + self.scaling = scaling + self.__interp = None + self._rom_weights = None + self.__maps = [] + + def addMap(self, m): + self.__maps.append(m) + + @property + def center(self): + return self.__center + + @center.setter + def center(self, value): + self.__center = value + self.__interp = None + + @property + def step_size(self): + return self.__step_size + + @step_size.setter + def step_size(self, value): + self.__step_size = value + self.__interp = None + + @property + def scaling(self): + return self.__scaling + + @scaling.setter + def scaling(self, value): + self.__scaling = value + self.__interp = None + + @property + def x(self): + return self.step_size[0] * (np.array(range(1, self.size[0]+1)) - self.center[0]) + + @property + def y(self): + return self.step_size[1] * (np.array(range(1, self.size[1]+1))- self.center[1]) + + @property + def size(self): + return self.__maps[0].data.shape + + @property + def offset(self): + return np.array(self.step_size)*(np.array(self.center) - 1/2. - np.array(self.size)/2.0) + + @property + def ROMWeights(self): + return self._rom_weights + + def z_xy(self, wavelength=1064e-9, direction="reflection", nr1=1.0, nr2=1.0): + + z_xy = np.ones(self.size, dtype=np.complex128) + + for m in self.__maps: + z_xy *= m.z_xy(wavelength=wavelength, direction=direction, nr1=nr1, nr2=nr2) + + return z_xy + + def generateROMWeights(self, EIxFilename, EIyFilename=None, verbose=False, interpolate=False): + if interpolate == True: + # Use EI nodes to interpolate if we + with open(EIxFilename, 'rb') as f: + EIx = pickle.load(f) + + if EIyFilename is None: + EIy = EIx + else: + with open(EIyFilename, 'rb') as f: + EIy = pickle.load(f) + + x = EIx.x + x.sort() + nx = np.unique(np.hstack((x, -x[::-1]))) + + y = EIy.x + y.sort() + ny = np.unique(np.hstack((y, -y[::-1]))) + + self.interpolate(nx, ny) + + self._rom_weights = makeWeightsNew(self, EIxFilename, EIyFilename, verbose=verbose) + return self.ROMWeights + + def interpolate(self, nx, ny, **kwargs): + """ + Interpolates all the maps that are used to fc + + Uses scipy.interpolate.interp2d and any keywords arguments are + passed on to it, thus settings like interpolation type and + fill values can be set. + + The range of nx and ny must contain the value zero so that the + center point of the map can be set. + """ + + for m in self.__maps: + m.interpolate(nx, ny) + + def plot(self, mode="absorption", show=True, clabel=None, xlim=None, ylim=None, wavelength=1064e-9): + + import pylab + + if xlim is not None: + _x = np.logical_and(self.x<=max(xlim)/100.0, self.x>=min(xlim)/100.0) + xmin = np.min(np.where(_x == True)) + xmax = np.max(np.where(_x == True)) + else: + xmin = 0 + xmax = len(self.x)-1 + xlim = [self.x.min()*100, self.x.max()*100] + + if ylim is not None: + _y = np.logical_and(self.y<=max(ylim)/100.0, self.y>=min(ylim)/100.0) + ymin = np.min(np.where(_y == True)) + ymax = np.max(np.where(_y == True)) + else: + ymin = 0 + ymax = len(self.y)-1 + ylim = [self.y.min()*100, self.y.max()*100] + + if mode == "absorption": + # plots how much of field is absorbed + data = 1-np.abs(self.z_xy()) + elif mode == "meter": + # plot the phase in terms of meters of displacement + k = 2*np.pi/wavelength + data = np.angle(self.z_xy()) / (2*k) + + zmin = data[xmin:xmax,ymin:ymax].min() + zmax = data[xmin:xmax,ymin:ymax].max() + + # 100 factor for scaling to cm + xrange = 100*self.x + yrange = 100*self.y + + fig = pylab.figure() + axes = pylab.pcolormesh(xrange, yrange, data, vmin=zmin, vmax=zmax) + pylab.xlabel('x [cm]') + pylab.ylabel('y [cm]') + + if xlim is not None: pylab.xlim(xlim) + if ylim is not None: pylab.ylim(ylim) + + pylab.title('Merged map {0}, mode {1}'.format(self.name, mode)) + + cbar = fig.colorbar(axes) + cbar.set_clim(zmin, zmax) + if clabel is not None: + cbar.set_label(clabel) + + if show: + pylab.show() + + return fig + class aperturemap(surfacemap): def __init__(self, name, size, step_size, R): diff --git a/pykat/optics/romhom.py b/pykat/optics/romhom.py index df26f13b37be02eafaa964627b10545c0b966c4f..1f6a9e635f908e7d41b1384b6970a22921178aac 100644 --- a/pykat/optics/romhom.py +++ b/pykat/optics/romhom.py @@ -323,7 +323,7 @@ def makeEmpiricalInterpolant(RB, sort=False): nodes = np.array(x_nodes, dtype=np.float64) if sort: - print sort + print (sort) ix = np.argsort(nodes) nodes = nodes[ix] node_indices = node_indices[ix] @@ -513,7 +513,7 @@ def _write_TS(queue, filename, tssize): i += 1 if i % 1000 == 0: - print i/float(tssize) + print(i/float(tssize)) hfile.flush() finally: hfile.close() @@ -542,8 +542,8 @@ def CreateTrainingSetHDF5(filename, maxOrder, z, w0, R, halfMapSamples, newtonCo nModes = 0 - for n in xrange(0, maxOrder+1): - for m in xrange(0, maxOrder+1): + for n in range(0, maxOrder+1): + for m in range(0, maxOrder+1): if n+m <= maxOrder and n <= m: nModes += 1 @@ -563,10 +563,10 @@ def CreateTrainingSetHDF5(filename, maxOrder, z, w0, R, halfMapSamples, newtonCo hfile.close() # make sure it's closed before - print "Starting processes..." + print("Starting processes...") # Have a bunch of processes doing the computation and one doing the writing iprocesses = [Process(target=_compute_TS, name="irom%i" % i, args=(iq, oq, x, w)) for i in range(NProcesses)] - oprocess = Process(target=_write_TS, name="orom%i" % i, args=(oq, filename, tssize)) + oprocess = Process(target=_write_TS, name="orom", args=(oq, filename, tssize)) oprocess.start() @@ -574,10 +574,10 @@ def CreateTrainingSetHDF5(filename, maxOrder, z, w0, R, halfMapSamples, newtonCo for P in iprocesses: P.start() - print "Filling queue..." + print("Filling queue...") curr = 0 - for n in xrange(0, maxOrder+1): - for m in xrange(0, maxOrder+1): + for n in range(0, maxOrder+1): + for m in range(0, maxOrder+1): if n+m <= maxOrder and n <= m: for _z in z: for _w0 in w0: @@ -760,10 +760,10 @@ def MakeROMFromHDF5(hdf5Filename, greedyFilename=None, EIFilename=None, tol=1e-1 next_RB_index = max_idx[np.argmax(max_res)] if worst_error <= tol: - print "Final basis size = %d, Final error = %e, Tolerance=%e" % (k, worst_error, tol) + print( "Final basis size = %d, Final error = %e, Tolerance=%e" % (k, worst_error, tol) ) break - print "worst error = %e at %i on iteration %d" % (worst_error, next_RB_index, k) + print ("worst error = %e at %i on iteration %d" % (worst_error, next_RB_index, k)) epsilon = TS['%i/data' % next_RB_index].value res = epsilon - emp_interp(B, epsilon, EI_indices) @@ -782,7 +782,7 @@ def MakeROMFromHDF5(hdf5Filename, greedyFilename=None, EIFilename=None, tol=1e-1 _TS = TS[str(next_RB_index)] f.write("%15.15e %15.15e %i %i\n" % (_TS["w0"].value, _TS["z"].value, _TS["n1"].value, _TS["n2"].value)) - print time.time() - start, "Seconds" + print (time.time() - start, "Seconds") if NProcesses > 1: for P in procs: P.terminate() @@ -791,7 +791,7 @@ def MakeROMFromHDF5(hdf5Filename, greedyFilename=None, EIFilename=None, tol=1e-1 greedyFilenameBase = os.path.splitext(greedyFilename)[0] - print "Writing to %s" % greedyFilename + print ("Writing to %s" % greedyFilename) EI = EmpiricalInterpolant(B=np.matrix(B).real, nodes=np.array(x_nodes).squeeze(), @@ -803,7 +803,7 @@ def MakeROMFromHDF5(hdf5Filename, greedyFilename=None, EIFilename=None, tol=1e-1 with open("%s.p" % EIFilename, 'wb') as f: pickle.dump(EI, f) - print "Writing to %s.p" % EIFilename + print ("Writing to %s.p" % EIFilename) return EI