Skip to content
Snippets Groups Projects
Select Git revision
  • 90dad473a344d34bdd01f1dc07c11af9580fd6df
  • master default protected
  • 72-improve-docs-for_optimal_setup
  • os-path-join
  • develop-GA
  • add-higher-spindown-components
  • v1.3
  • v1.2
  • v1.1.2
  • v1.1.0
  • v1.0.1
11 results

mcmc_based_searches.py

Blame
  • maps.py 81.08 KiB
    """
    ------------------------------------------------------
    Utility functions for handling mirror surface
    maps. Some functions based on earlier version
    in Matlab (http://www.gwoptics.org/simtools/)
    Work in progress, currently these functions are
    untested!
    
    http://www.gwoptics.org/pykat/
    ------------------------------------------------------
    """
    from __future__ import absolute_import
    from __future__ import division
    from __future__ import print_function
    
    from pykat.optics.romhom import makeWeightsNew
    from scipy.interpolate import interp2d, interp1d
    from scipy.optimize import minimize
    from pykat.maths.zernike import *        
    from pykat.exceptions import BasePyKatException
    from copy import deepcopy
    
    import numpy as np
    import math
    import pickle
    
    
    class MirrorROQWeights:
        
        def __init__(self, rFront, rBack, tFront, tBack):
            self.rFront = rFront
            self.rBack = rBack
            self.tFront = tFront
            self.tBack = tBack
        
        def writeToFile(self, romfilename):
            with open(romfilename + ".rom", "w+") as f:
                if self.rFront is not None: self.rFront.writeToFile(f=f)
                if self.rBack  is not None: self.rBack.writeToFile(f=f)
                if self.tFront is not None: self.tFront.writeToFile(f=f)
                if self.tBack  is not None: self.tBack.writeToFile(f=f)
                
    class surfacemap(object):
        
        def __init__(self, name, maptype, size, center, step_size, scaling, data=None,
                     notNan=None, Rc=None, zOffset=None, xyOffset=(.0,.0)):
            '''
            size, center, step_size, xyOffset are all tuples of the form (x, y),
            i.e., (col, row).
            '''
            
            self.name = name
            self.type = maptype
            # Currently "beam center", i.e., mirror_center + xyOffset. 
            self.center = center
            self.step_size = step_size
            self.scaling = scaling
            self.notNan = notNan
            self.Rc = Rc
            # Offset of fitted sphere. Proably unnecessary to have here.
            self.zOffset = zOffset
            self.__interp = None
            self._zernikeRemoved = {}
            self._betaRemoved = None
            self._xyOffset = xyOffset
    		
            if data is None:
                self.data = np.zeros(size[::-1])
            else:
                self.data = data
    
            if notNan is None:
                self.notNan = np.ones(size[::-1], dtype=bool)
            else:
                self.notNan = notNan
    
            self._rom_weights = None
            
        def write_map(self, filename):
            with open(filename,'w') as mapfile:
                
                mapfile.write("% Surface map\n")
                mapfile.write("% Name: {0}\n".format(self.name))
                mapfile.write("% Type: {0}\n".format(self.type))
                mapfile.write("% Size: {0} {1}\n".format(self.data.shape[0], self.data.shape[1]))
                mapfile.write("% Optical center (x,y): {0} {1}\n".format(self.center[0], self.center[1]))
                mapfile.write("% Step size (x,y): {0} {1}\n".format(self.step_size[0], self.step_size[1]))
                mapfile.write("% Scaling: {0}\n".format(float(self.scaling)))
                mapfile.write("\n\n")
                
                for i in range(0, self.data.shape[0]):
                    for j in range(0, self.data.shape[1]):
                        mapfile.write("%.15g " % self.data[i,j])
                    mapfile.write("\n")
    
        @property
        def xyOffset(self):
            '''
            The offset from the mirror center measured in meters.
            '''
            return self._xyOffset
    
        @xyOffset.setter
        def xyOffset(self,offset):
            '''
            Give 'offset' is in meters.
            '''
            
            x0new = self.center[0] + (-self.xyOffset[0] + offset[0])/self.step_size[0]
            y0new = self.center[1] + (-self.xyOffset[1] + offset[1])/self.step_size[1]
    
            # Checking so new suggested "center" is within the the mirror surface
            if (x0new>0 and x0new<self.size[0]-1 and y0new>0 and y0new<self.size[1]-1 and
                self.notNan[round(x0new),round(y0new)]):
                
                self.center = (x0new, y0new)
                self._xyOffset = (offset[0], offset[1])
            
            else:
                print(('Error in xyOffset: ({:.2e}, {:.2e}) m --> pos = ({:.2f}, {:.2f}) '+
                      'is outside mirror surface.').format(offset[0],offset[1],x0new,y0new))
                
    
        @property
        def betaRemoved(self):
            return self._betaRemoved
        
        @betaRemoved.setter
        def betaRemoved(self, v):
            self._betaRemoved = v
            
        @property
        def zernikeRemoved(self):
            return self._zernikeRemoved
        
        @zernikeRemoved.setter
        def zernikeRemoved(self, v):
            '''
            v = tuple(m,n,amplitude)
            '''
            self._zernikeRemoved["%i%i" % (v[0],v[1])] = v
    
        @property
        def data(self):
            return self.__data
        
        @data.setter
        def data(self, value):
            self.__data = value
            self.__interp = None
        
        @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
    
        # CHANGED! I swapped: self.data.shape[0]  <----> self.data.shape[1]. Because
        # the rows/columns in the data matrix are supposed to be y/x-values(?).
        # This may mess things up in other methods though... I fixed the plot-method.
        # /DT
        @property
        def x(self):
            return self.step_size[0] * (np.array(range(0, self.data.shape[1])) - self.center[0])
        
        @property
        def y(self):
            return self.step_size[1] * (np.array(range(0, self.data.shape[0]))- self.center[1])
            
        # CHANGED! Since everything else (step_size, center) are given as (x,y), and not
        # as (row, column), I changed this to the same format. /DT
        @property
        def size(self):
            return self.data.shape[::-1]
            
        @property
        def offset(self):
            return np.array(self.step_size)*(np.array(self.center) - (np.array(self.size)-1)/2.0)
        
        @property
        def ROMWeights(self):
            return self._rom_weights
        
        def z_xy(self, x=None, y=None, wavelength=1064e-9, direction="reflection_front", nr1=1.0, nr2=1.0):
            """
            For this given map the field perturbation is computed. This data
            is used in computing the coupling coefficient. It returns a grid
            of complex values representing the change in amplitude or phase
            of the field.
            
                x, y      : Points to interpolate at, 'None' for no interpolation.
                
                wavelength: Wavelength of light in vacuum [m]
                
                direction : Sets which distortion to return, as beams travelling
                            in different directions will see different distortions.
                            Options are:
                                    "reflection_front"
                                    "transmission_front" (front to back)
                                    "transmission_back" (back to front)
                                    "reflection_back"
                                    
                nr1       : refractive index on front side
                
                nr2       : refractive index on back side
                
            """
            
            assert(nr1 >= 1)
            assert(nr2 >= 1)
            
            if x is None and y is None:
                data = self.scaling * self.data
            else:
                if self.__interp is None:
                    self.__interp = interp2d(self.x, self.y, self.data * self.scaling)
                    
                data = self.__interp(x, y)
            
            if direction == "reflection_front" or direction == "reflection_back":
                if "phase" in self.type:
                    k = math.pi * 2 / wavelength
                    
                    if direction == "reflection_front":
                        return np.exp(-2j * nr1 * k * data)
                    else:
                        return np.exp(2j * nr2 * k * data[:,::-1])
                    
                elif "absorption" in self.type:
                    if direction == "reflection_front":
                        return np.sqrt(1.0 - data)
                    else:
                        return np.sqrt(1.0 - data[:, ::-1])
                elif "surface_motion" in self.type:
                    if direction == "reflection_front":
                        return data
                    else:
                        return data[:, ::-1]
                else:
                    raise BasePyKatException("Map type needs handling")
                    
            elif direction == "transmission_front" or direction == "transmission_back":
                if "phase" in self.type:
                    k = math.pi * 2 / wavelength
                    
                    if direction == "transmission_front":
                        return np.exp((nr1-nr2) * k * data)
                    else:
                        return np.exp((nr2-nr1) * k * data[:, ::-1])
                    
                elif "absorption" in self.type:
                    if direction == "transmission_front":
                        return np.sqrt(1.0 - data)
                    else:
                        return np.sqrt(1.0 - data[:, ::-1])
                        
                elif "surface_motion" in self.type:
                    return np.ones(data.shape)
                else:
                    raise BasePyKatException("Map type needs handling")
                    
            else:
                raise ValueError("Direction not valid")
            
    
        
        def generateROMWeights(self, EIxFilename, EIyFilename=None, nr1=1.0, nr2=1.0, verbose=False, interpolate=False, newtonCotesOrder=8):
            
            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)
            
            w_refl_front, w_refl_back, w_tran_front, w_tran_back = (None, None, None, None)
            
            if "reflection" in self.type or "both" in self.type:
                w_refl_front = makeWeightsNew(self, EIxFilename, EIyFilename,
                                          verbose=verbose, newtonCotesOrderMapWeight=newtonCotesOrder,
                                          direction="reflection_front")
                
                w_refl_front.nr1 = nr1
                w_refl_front.nr2 = nr2
                
                w_refl_back = makeWeightsNew(self, EIxFilename, EIyFilename,
                                          verbose=verbose, newtonCotesOrderMapWeight=newtonCotesOrder,
                                          direction="reflection_back")
                
                w_refl_back.nr1 = nr1
                w_refl_back.nr2 = nr2
    
            if "transmission" in self.type or "both" in self.type:                                      
                w_tran_front = makeWeightsNew(self, EIxFilename, EIyFilename,
                                          verbose=verbose, newtonCotesOrderMapWeight=newtonCotesOrder,
                                          direction="transmission_front")
    
                w_refl_front.nr1 = nr1
                w_refl_front.nr2 = nr2
                                                
                w_tran_back  = makeWeightsNew(self, EIxFilename, EIyFilename,
                                          verbose=verbose, newtonCotesOrderMapWeight=newtonCotesOrder,
                                          direction="transmission_back")
                
                w_refl_back.nr1 = nr1
                w_refl_back.nr2 = nr2
                
            self._rom_weights = MirrorROQWeights(w_refl_front, w_refl_back, w_tran_front, w_tran_back)
            
            return self._rom_weights
                
        def interpolate(self, nx, ny, **kwargs):
            """
            Interpolates the map for some new x and y values.
            
            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.
            """
    
            D = interp2d(self.x, self.y, self.data, **kwargs)
            
            data = D(nx-self.offset[0], ny-self.offset[1])
            
            Dx = interp1d(nx, np.arange(0,len(nx)))
            Dy = interp1d(ny, np.arange(0,len(ny)))
            
            self.center = (Dx(0), Dy(0))
            self.step_size = (nx[1]-nx[0], ny[1]-ny[0])
            self.data = data
    
        # xlim and ylim given in centimeters
        def plot(self, show=True, clabel=None, xlim=None, ylim=None):
            import pylab
            
            if xlim is not None:
                # Sorts out the x-values within xlim
                _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:
                # Uses the whole available x-range
                xmin = 0
                xmax = len(self.x)-1
                xlim = [self.x.min()*100, self.x.max()*100]
        
            if ylim is not None:
                # Sorts out the y-values within ylim
                _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:
                # Uses the whole available y-range
                ymin = 0
                ymax = len(self.y)-1
                ylim = [self.y.min()*100, self.y.max()*100]
                
            # CHANGED! y corresponds to rows and x to columns, so this should be
            # correct now. Switch the commented/uncommented things to revert. /DT
            # ------------------------------------------------------
            # min and max of z-values
            zmin = self.data[ymin:ymax,xmin:xmax].min()
            zmax = self.data[ymin:ymax,xmin:xmax].max()
            # zmin = self.data[xmin:xmax,ymin:ymax].min()
            # zmax = self.data[xmin:xmax,ymin:ymax].max()
            # ------------------------------------------------------
            
            # 100 factor for scaling to cm
            xRange = 100*self.x
            yRange = 100*self.y
            # xRange, yRange = np.meshgrid(xRange,yRange)
            
            fig = pylab.figure()
    
            # Important to remember here is that xRange corresponds to column and
            # yRange to row indices of the matrix self.data.
            pcm = pylab.pcolormesh(xRange, yRange, self.data)
            pcm.set_rasterized(True)
            
            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('Surface map {0}, type {1}'.format(self.name, self.type))
    
            cbar = pylab.colorbar()
            #cbar.set_clim(zmin, zmax)
            
            if clabel is not None:
                cbar.set_label(clabel)
        
            if show:
                pylab.show()
            
            return fig
    
        def rescaleData(self,scaling):
            '''
            Rescaling data. Assumes that the current map is scaled correctly, i.e.,
            surfacemap.scaling*surfacemap.data would give the data in meters.
    
            Input: scaling
            scaling  - float number that representing the new scaling. E.g.,
                       scaling = 1.0e-9 means that surfacemap.data is converted to
                       being given in nano meters.
            '''
            if isinstance(scaling, float):
                self.data[self.notNan] = self.data[self.notNan]*self.scaling/scaling
                self.scaling = scaling
            else:
                print('Error: Scaling factor must be of type float, {:s} found'.
                      format(type(scaling)))
                
    
        def rms(self, w=None):
            '''
            Returns the root mean square value of the mirror map, expressed in the unit
            given by self.scaling. If w [m] is specified, only the area inside the radius
            is considered. The area is r<w is centered around the self.center, i.e., the
            beam center (not mirror center).
            '''
            if w is None:
                return math.sqrt((self.data[self.notNan]**2).sum())/self.notNan.sum()
            else:
                R = self.find_radius(unit='meters')
                if w>=R:
                    return math.sqrt((self.data[self.notNan]**2).sum())/self.notNan.sum()
                else:
                    rho = self.createPolarGrid()[0]
                    inside = np.zeros(self.data.shape,dtype=bool)
                    tmp = rho<w
                    inside[tmp] = self.notNan[tmp]
                    return math.sqrt((self.data[inside]**2).sum())/inside.sum()
                
        def avg(self, w=None):
            '''
            Returns the average height of the mirror map, expressed in the unit
            given by self.scaling. If w [m] is specified, only the area inside the radius
            is considered. The area is r<w is centered around the self.center, i.e., the
            beam center (not mirror center).
            '''
            if w is None:
                tot = self.data[self.notNan].sum()
                sgn = np.sign(tot)
                return sgn*math.sqrt(sgn*tot)/self.notNan.sum()
            else:
                R = self.find_radius(unit='meters')
                if w>=R:
                    tot = self.data[self.notNan].sum()
                    sgn = np.sign(tot)
                    return sgn*math.sqrt(sgn*tot)/self.notNan.sum()
                else:
                    rho = self.createPolarGrid()[0]
                    inside = np.zeros(self.data.shape,dtype=bool)
                    tmp = rho<w
                    inside[tmp] = self.notNan[tmp]
                    tot = self.data[inside].sum()
                    sgn = np.sign(tot)
                    return sgn*math.sqrt(sgn*tot)/inside.sum()
                
                
        def find_radius(self, method = 'max', unit='points'):
            '''
            Estimates the radius of the mirror in the xy-plane. If xyOffset is different from
            (0,0), then the radius will be the minimum distance from the beam center to the
            nearest mirror edge.
            
            method   - 'max' gives maximal distance from beam centre to the edge.
                     - 'xy' returns the distance from centre to edge along x- and y-directions.
                       Doesn't make sense if the beam center is offset from the mirror center.
                     - 'area' calculates the area and uses this to estimate the mean radius.
                     - 'min' gives minimal distance from beam center to the edge
            unit     - 'points' gives radius in data points.
                       'meters' gives radius in meters.
    
            Based on Simtools function 'FT_find_map_radius.m' by Charlotte Bond.
            '''
    
            # Row and column indices of non-NaN
            rIndx, cIndx = self.notNan.nonzero()
            if unit == 'meters' or unit == 'metres' or unit=='Meters' or unit=='Metres':
                # Offsets so that (0,0) is in the center.
                x = self.step_size[0]*(cIndx - self.center[0])
                y = self.step_size[1]*(rIndx - self.center[1])
            else:
                # Offsets so that (0,0) is in the center.
                x = cIndx - self.center[0]
                y = rIndx - self.center[1]
                
            # Maximum distance from center to the edge of the mirror.
            if method=='max' or method=='Max' or method=='MAX':
                r = math.sqrt((x**2 + y**2).max())
            # x and y radii
            elif method=='xy' or method=='XY' or method == 'yx' or method == 'YX':
                r = []
                r.append(abs(x).max()+0.5)
                r.append(abs(y).max()+0.5)
            # Mean radius by dividing the area by pi. Should change this in case
            # x and y step sizes are different.
            elif method=='area' or method=='Area' or method=='AREA':
                if unit == 'meters' or unit == 'metres' or unit=='Meters' or unit=='Metres':
                    r = step_size[0]*math.sqrt(len(cIndx)/math.pi)
                else:
                    r = math.sqrt(len(cIndx)/math.pi)
            elif method=='min':
                # Arrays of row and column indices where data is NaN.
                r,c = np.where(self.notNan == False)
                # Sets the map center to index [0,0]
                x = c-self.center[0]
                y = r-self.center[1]
                if unit == 'meters' or unit == 'metres' or unit=='Meters' or unit=='Metres':
                    x = x*self.step_size[0]
                    y = y*self.step_size[1]
                # Minimum radius squared
                r = math.sqrt((x**2+y**2).min())
            return r
    
    
        def zernikeConvol(self,n,m=None):
            '''
            Performs a convolution with the map surface and Zernike polynomials.
    
            Input: n, m
            n  -  If m (see below) is set to None, n is the maximum order of Zernike
                  polynomials that will be convolved with the map. If m is an integer
                  or a list, only the exact order n will be convolved.
            m  -  = None. All Zernike polynomials up to and equal to order n will
                  be convolved.
                  = value [int]. The Zernike polynomial characterised by (n,m) will
                  be convolved witht the map.
                  = list of integers. The Zernike polynomials of order n with m-values
                  in m will be convolved with the map.
                 
            Returns: A
            A  -  If m is None:
                    List with lists of amplitude coefficients. E.g. A[n] is a list of
                    amplitudes for the n:th order with m ranging from -n to n, i.e.,
                    A[n][0] is the amplitude of the Zernike polynomial Z(m,n) = Z(-n,n).
                  If m is list or int:
                    Returns list with with Zernike amplitudes ordered in the same way as
                    in the input list m. If m == integer, there is only one value in the
                    list A.
            
            Based on the simtool function 'FT_zernike_map_convolution.m' by Charlotte
            Bond.
            '''
            
            mVals = []
            if m is None:
                nVals = range(0,n+1)
                for k in nVals:
                    mVals.append(range(-k,k+1,2))
            elif isinstance(m,list):
                nVals = range(n,n+1)
                mVals.append(m)
            elif isinstance(m,int):
                nVals = range(n,n+1)
                mVals.append(range(m,m+1))
                
            # Amplitudes
            A = []
            # Radius
            R = self.find_radius(unit='meters')
            
            # Grid for for Zernike-polynomials (same size as map).
            rho, phi = self.createPolarGrid()
            rho = rho/R
            # Loops through all n-values up to n_max
            for k in range(len(nVals)):
                n = nVals[k]
                # Loops through all possible m-values for each n.
                tmp = []
                for m in mVals[k]:
                    # (n,m)-Zernike polynomial for this grid.
                    Z = zernike(m, n, rho, phi)
                    # Z = znm(n, m, rho, phi)
                    # Convolution (unnecessary conjugate? map.data is real numbers.)
                    c = (Z[self.notNan]*np.conjugate(self.data[self.notNan])).sum()
                    cNorm = (Z[self.notNan]*np.conjugate(Z[self.notNan])).sum()
                    c = c/cNorm
                    tmp.append(c)
                A.append(tmp)
            if len(A) == 1:
                A = A[0]
                if len(A)==1:
                    A = A[0]
                    
            return A
    
        def rmZernikeCurvs(self, zModes='all',interpol=False):
            '''
            Removes curvature from mirror map by fitting Zernike polynomials to the 
            mirror surface. Also stores the estimated radius of curvature
            and the Zernike polynomials and amplitudes.
            
            zModes   - 'defocus' uses only Zernike polynomial (n,m) = (2,0), which has a
                       parabolic shape.
                       'astigmatism' uses Zernike polynomials (n,m) = (2,-2) and (n,m) = (2,2).
                       'all' uses both defocus and the astigmatic modes.
            interpol - Boolean where 'true' means that the surface map is interpolated to get
                       more data points.
    
            Returns [Rc,A]
            Rc       - Estimated radius of curvature removed [m].
            A        - Amplitudes of Zernike polynomials removed [surfacemap.scaling].
    
            Based on the Simtools functions 'FT_remove_zernike_curvatures_from_map.m',
            'FT_zernike_map_convolution.m', etc. by Charlotte Bond.
            '''
            
            tmp = deepcopy(self)
            ny,nx = self.data.shape
            # Interpolating for more precise convolution, in case size of map is small. Not sure
            # this is necessary. It isn't used in the latest Simtools tutorials.
            if interpol and (nx<1500 or ny<1500):
                # Number of extra steps inserted between two adjacent data points.
                N = math.ceil(1500.0/min(nx,ny))
                # New arrays of x-values and y-values
                x = np.arange(tmp.x[0],tmp.x[-1]+tmp.step_size[0]/(N+1),tmp.step_size[0]/N)
                y = np.arange(tmp.y[0],tmp.y[-1]+tmp.step_size[1]/(N+1),tmp.step_size[1]/N)
    
                '''
                # --------------------------------------------------------------
                # Interpolating data
                tmp.interpolate(x,y)
                tmp.plot()
                # Interpolating for notNan.
                g = interp2d(self.x, self.y, self.notNan, kind='linear')
                tmp.notNan = np.round(g(x,y)) # round() or floor() here?! Can't decide...
                # Converting binary to boolean
                tmp.notNan = tmp.notNan==1
                # --------------------------------------------------------------
                '''
                # Later when compatible with interpolation function, switch code
                # below for code above (but need to change 
                # --------------------------------------------------------------
                # Interpolation functions, for the data (f) and for notNan (g).
                f = interp2d(tmp.x, tmp.y, tmp.data, kind='linear')
                g = interp2d(tmp.x, tmp.y, tmp.notNan, kind='linear')
                # Interpolating
                tmp.data = f(x,y)
    
                tmp.notNan = np.round(g(x,y)) # round() or floor() here?! Can't decide...
                # Converting binary to boolean
                tmp.notNan = tmp.notNan==1
    
                # Setting new step size
                tmp.step_size = (x[1]-x[0],y[1]-y[0])
                # Setting new center
                fx = interp1d(x, np.arange(0,len(x)))
                fy = interp1d(y, np.arange(0,len(y)))
                tmp.center = (fx(0.0), fy(0.0))
                # --------------------------------------------------------------
    
            # Convolution between the surface map and Zernike polynomials
            # to get amplitudes of the polynomials. Only interested in n=2.
            A = tmp.zernikeConvol(2)[2]
            
            # Radius of mirror.
            R = self.find_radius(unit='meters')
            # Grid in polar coordinates
            rho,phi = self.createPolarGrid()
            rho = rho/R
            
            # Creates the choosen Zernike polynomials and removes them from the
            # mirror map.
            if zModes=='all' or zModes=='All':
                for k in range(0,3):
                    m = (k-1)*2
                    Z = A[k]*zernike(m, 2, rho, phi)
                    self.data[self.notNan] = self.data[self.notNan]-Z[self.notNan]
                    self.zernikeRemoved = (m, 2, A[k])
                # Estimating radius of curvature
                Rc = znm2Rc([a*self.scaling for a in A], R)
            elif zModes=='astigmatism' or zModes=='Astigmatism':
                for k in range(0,3,2):
                    m = (k-1)*2
                    Z = A[k]*zernike(m, 2, rho, phi)
                    smap.data[self.notNan] = self.data[self.notNan]-Z[self.notNan]
                    self.zernikeRemoved = (m, 2, A[k])
                Rc = znm2Rc([a*self.scaling for a in A[::2]], R)
            elif zModes=='defocus' or zModes=='Defocus':
                Z = A[1]*zernike(0, 2, rho, phi)
                self.data[self.notNan] = self.data[self.notNan]-Z[self.notNan]
                self.zernikeRemoved = (0, 2, A[1])
                Rc = znm2Rc(A[1]*self.scaling, R)
            
            self.Rc = Rc
            
            return self.Rc, self.zernikeRemoved
    
    
        def rmTilt(self, method='fitSurf', w=None, xbeta=None, ybeta=None, zOff=None):
            
            R = self.find_radius(unit='meters')
            
            if method=='zernike' or method=='Zernike':
                A1 = self.rmZernike(1, m = 'all')
                # Calculating equivalent tilts in radians
                xbeta = np.arctan(A1[1]*self.scaling/R)
                ybeta = np.arctan(A1[0]*self.scaling/R)
                self.betaRemoved = (xbeta,ybeta)
                return A1, xbeta, ybeta
    
            else:
                X,Y = np.meshgrid(self.x,self.y)
                r2 = X**2 + Y**2
    
                if w is not None:
                    weight = 2/(math.pi*w**2) * np.exp(-2*r2[self.notNan]/w**2)
    
                def f(p):
                    # This is used in simtools, but gives less good results from what I have seen, also in simtools.
                    # I would guess the problem is that the initial guess is 0, and the minimising algorithm would
                    # need to go quite far away from zero to find a minumum even slightly away way from zero.
                    # I think the idea is to make the algorithm to search for small angles. Of course, if used, this
                    # must be converted back after minimisation.
                    # p[0] = p[0]*1.0e-9
                    # p[1] = p[1]*1.0e-9
                    Z = self.createSurface(0,X,Y,p[2],0,0,p[0],p[1])
                    if w is None:
                        res = math.sqrt(((self.data[self.notNan]-Z[self.notNan])**2).sum())/self.notNan.sum()
                    else:
                        # weight = 2/(math.pi*w**2) * np.exp(-2*r2[self.notNan]/w**2)
                        res = math.sqrt((weight*(self.data[self.notNan]-Z[self.notNan])**2).sum())/weight.sum()
    
                    return res
    
                if xbeta is None:
                    xbeta = 0
                if ybeta is None:
                    ybeta = 0
                if zOff is None:
                    zOff = 0
    
                params = [xbeta,ybeta,zOff]
    
                opts = {'xtol': 1.0e-9, 'ftol': 1.0e-9, 'maxiter': 1000, 'maxfev': 1000, 'disp': False}
                out = minimize(f, params, method='Nelder-Mead', options=opts)
                if not out['success']:
                    msg = '  Warning: ' + out['message'].split('.')[0] + ' (nfev={:d}).'.format(out['nfev'])
                    print(msg)
    
                xbeta = out['x'][0]
                ybeta = out['x'][1]
                zOff = out['x'][2]
    
                Z = self.createSurface(0,X,Y,zOff,0,0,xbeta,ybeta)
                self.data[self.notNan] = self.data[self.notNan] - Z[self.notNan]
                self.betaRemoved = (xbeta, ybeta)
    
                # Equivalent Zernike amplitude
                A1 = R*np.tan(np.array([ybeta,xbeta]))/self.scaling
                self.zernikeRemoved = (-1,1,A1[0])
                self.zernikeRemoved = (1,1,A1[1])
                
                return A1,xbeta,ybeta,zOff
            
    
        def rmSphericalSurf(self, Rc0, w=None, zOff=None, isCenter=[False,False]):
            '''
            Fits spherical surface to the mirror map and removes it.
    
            Inputs: Rc0, w, zOff, isCenter
            Rc       - Initial guess of the Radius of curvature. [m]
            w        - Gaussian weighting parameter. The distance from the center where the
                       weigths have decreased by a factor of exp(-2) compared to the center.
                       Should preferrably be equal to the beam radius at the mirror. [m]
            zOff     - Initial guess of the z-offset. Generally unnecessary. [surfacemap.scaling]
            isCenter - 2D-list with booleans. isCenter[0] Determines if the center of the
                       sphere is to be fitted (True) or not (False, recommended). If the center is
                       fitted, then isCenter[1] determines if the weights (in case w!=None) are
                       centered around the fitted center (True) or centered around the center of
                       the data-grid (False, highly recommended).
                       
            if isCenter[0] == False
            Returns: Rc, zOff
            Rc       - Radius of curvature of the sphere removed from the mirror map. [m]
            zOff     - z-offset of the sphere removed. [surfacemap.scaling]
            isCenter[0] == True
            Returns Rc, zOff, center
            center   - [x0, y0] giving the coordinates of the center of the fitted sphere.
            
            Based on the file 'FT_remove_curvature_from_mirror_map.m' by Charlotte Bond.
            '''
            
            # z-value at centre of the data-grid. Serves as initial guess for deviation
            # from z=0, if no other first guess is given.
            if zOff is None:
                zOff = self.data[round(self.center[1]), round(self.center[0])]
    
            # If fitting center of the sphere, four variables are fitted. Initial guess
            # of deviation from notNan-data-grid-center: (x0,y0) = (0,0).
            if isCenter[0]:
                p = [Rc0, zOff, 0.0, 0.0]
            # Otherwise two.
            else:
                p = [Rc0, zOff]
    
            # Grids with X,Y and r2 values. X and Y crosses zero in the center
            # of the xy-plane.
            X,Y = np.meshgrid(self.x, self.y)
            r2 = X**2 + Y**2
            
            # Cost-function to minimize.
            def costFunc(p):
                # If the center of the mirror is fitted, four variables are fitted.
                if isCenter[0]:
                    Rc = p[0]
                    zOff = p[1]
                    x0 = p[2]
                    y0 = p[3]
                # Otherwise 2 variables.
                else:
                    Rc = p[0]
                    zOff = p[1]
                    x0 = 0
                    y0 = 0
    
                Z = self.createSurface(Rc,X,Y,zOff,x0,y0)
    
                if w is None:
                    # Mean squared difference between map and the created sphere.
                    res = math.sqrt( ((self.data[self.notNan] -
                                       Z[self.notNan])**2).sum() )/self.notNan.sum()
                else:
                    if isCenter[0] and isCenter[1]:
                        # Weights centered around fitting spehere center. May give weird
                        # results if the mirror deviates much from a sphere.
                        weight = (2/(math.pi*w**2))*np.exp(-2*( (X-x0)**2 + (Y-y0)**2 )/w**2)
                    else:
                        # Weights centered around the center of the mirror xy-plane.
                        weight = (2/(math.pi*w**2))*np.exp(-2*r2/w**2)
                    # Weighted mean squared difference between map and the created sphere.
                    res = math.sqrt( ( weight[self.notNan]*
                                      ( (self.data[self.notNan] - Z[self.notNan])**2 )
                                      ).sum() ) \
                                      /weight[self.notNan].sum()
                return res
    
            # Using the simplex Nelder-Mead method. This is the same or very
            # similar to the method used in 'FT_remove_curvature_from_mirror_map.m',
            # but there are probably better methods to use.
            opts = {'xtol': 1.0e-5, 'ftol': 1.0e-9, 'maxiter': 1000, 'maxfev': 1000, 'disp': False}
            out = minimize(costFunc, p, method='Nelder-Mead', options=opts)
            if not out['success']:
                msg = '  Warning: ' + out['message'].split('.')[0] + ' (nfev={:d}).'.format(out['nfev'])
                print(msg)
                
            # Assigning values to the instance variables
            self.Rc = out['x'][0]
            if self.zOffset is None:
                self.zOffset = 0
            self.zOffset = self.zOffset + out['x'][1]
    
            # Equivalent Zernike (n=2,m=0) amplitude.
            R = self.find_radius(unit='meters')
            A20 = Rc2znm(self.Rc,R)/self.scaling
            self.zernikeRemoved = (0,2,A20)
    
            # If center was fitted, assign new values to instance variable center, and
            # subtract the fitted sphere from the mirror map.
            if isCenter[0]:
                x0 = out['x'][2]
                y0 = out['x'][3]
                # Converts the deviation into a new absolut center in data points.
                self.center = (self.center[0] + x0/self.step_size[0],
                               self.center[1] + y0/self.step_size[1])
                # Creating fitted sphere
                Z = self.createSurface(self.Rc, X, Y, self.zOffset, x0, y0)
                # Subtracting sphere from map
                self.data[self.notNan] = self.data[self.notNan]-Z[self.notNan]
                return self.Rc, self.zOffset, self.center, A20
            # Subtracting fitted sphere from mirror map.
            else:
                # Creating fitted sphere
                Z = self.createSurface(self.Rc,X,Y,self.zOffset)
                # Subtracting sphere from map
                self.data[self.notNan] = self.data[self.notNan]-Z[self.notNan]
                return self.Rc, self.zOffset, A20
    
        def remove_curvature(self, method='zernike', w=None, zOff=None,
                             isCenter=[False,False], zModes = 'all'):
            '''
            Removes curvature from mirror map by fitting a sphere to the mirror map, or
            by convolving with second order Zernike polynomials and then extracting
            these polynomials with the obtained amplitudes. 
    
            Inputs: method, w, zOff, isCenter, zModes
            method   - 'zernike' (default) convolves second order Zernike polynomials with the
                       map, and then removes the polynomial with the obtained amplitude from the
                       surface map. Which of the three second order modes that are fitted is
                       determined by zMods (see below).
                     - 'sphere' fits a sphere to the mirror map. The fitting is Gaussian
                       weighted if w (see below) is specifed. 
            w        - Gaussian weighting parameter. The distance from the center where the
                       weigths have decreased by a factor of exp(-2) compared to the center.
                       Should preferrably be equal to the beam radius at the mirror. [m]
            zOff     - Initial guess of the z-offset. Only used together with method='sphere'.
                       Generally unnecessary. [surfacemap.scaling]
            isCenter - 2D-list with booleans. isCenter[0] Determines if the center of the
                       sphere is to be fitted (True) or not (False, recommended). If the center is
                       fitted, then isCenter[1] determines if the weights (in case w!=None) are
                       centered around the fitted center (True) or centered around the center of
                       the data-grid (False, highly recommended).
            zModes   - 'defocus' uses only Zernike polynomial (n,m) = (2,0), which has a
                       parabolic shape.
                       'astigmatism' uses Zernike polynomials (n,m) = (2,-2) and (n,m) = (2,2).
                       'all' uses both the defocus and the astigmatic modes.
    
            if method == 'zernike'
            Returns: Rc, znm
            Rc       - Approximate spherical radius of curvature removed.
            znm      - Dictionary specifying which Zernike polynomials that have been removed
                       as well as their amplitudes.
            if method == 'sphere' and isCenter[0] == False
            Returns: Rc, zOff
            Rc       - Radius of curvature of the sphere removed from the mirror map. [m]
            zOff     - z-offset of the sphere removed. [surfacemap.scaling]
            if method == 'sphere' and isCenter[0] == True
            Returns Rc, zOff, center
            center   - [x0, y0] giving the coordinates of the center of the fitted sphere. 
            '''
            
            
            if method == 'zernike' or method == 'Zernike':
                Rc, znm = self.rmZernikeCurvs(zModes)
                return Rc, znm
            elif method == 'sphere' or method == 'Sphere':
                smap = deepcopy(self)
                # Using Zernike polynomial (m,n) = (0,2) to estimate radius of curvature.
                Rc0 = smap.rmZernikeCurvs('defocus')[0]
                if isCenter[0]:
                    Rc, zOff, center, A20 = self.rmSphericalSurf(Rc0, w, zOff, isCenter)
                    return Rc, zOff, self.center, A20
                else:
                    Rc, zOff, A20 = self.rmSphericalSurf(Rc0, w, zOff, isCenter)
                    return Rc, zOff, A20
    
                
        def recenter(self):
            '''
            Setting center of mirror to be in the middle of the notNan field. 
            
            Based on 'FT_recenter_mirror_map.m' by Charlotte Bond.
            '''
            # Row and column indices with non-NaN elements
            rIndx, cIndx = self.notNan.nonzero()
            # Finding centres
            x0 = float(cIndx.sum())/len(cIndx)
            y0 = float(rIndx.sum())/len(rIndx)
            self.center = tuple([x0,y0])
            self.xyOffset = (.0,.0)
            return self.center
            # -------------------------------------------------
            
        def removePeriphery(self):
            '''
            Finds the NaN elements closest to the map center, and sets all
            elements further out to NaN. Then cropping and recentering.
            
            Based on FT_remove_elements_outside_map.m by Charlotte Bond.
            '''
            # Saves the xyOffset to restore it afterwards
            xyOffset = self.xyOffset
            self.recenter()
            # Arrays of row and column indices where data is NaN.
            r,c = np.where(self.notNan == False)
            # Sets the map center to index [0,0]
            x = c-self.center[0]
            y = r-self.center[1]
            # Minimum radius squared measured in data points.
            r2 = (x**2+y**2).min()
            # Grid with the same dimensions as the map.
            X,Y = np.meshgrid(self.x, self.y)
            # Grid containing radial distances squared measured from the center.
            R2 = (X/self.step_size[0])**2 + (Y/self.step_size[1])**2
            # Matrix with True=='outside mirror surface'
            outs = R2>=r2
            # Setting notNan-matrix to false outside the mirror.
            self.notNan[outs] = False
            # Setting data matrix to 0 outside the mirror.
            self.data[outs] = 0.0
    
            # Removing unnecessary data points.
            # --------------------------------------
            # Radius inside data is kept. Don't know why R is set this way,
            # but trusting Dr. Bond.
            R = round(math.ceil(2.0*math.sqrt(r2)+6.0)/2.0)
            if (2*R<self.data.shape[0] and 2*R<self.data.shape[1] and
                R<self.center[0] and R<self.center[1] and
                R<(self.size[0]-self.center[0]) and R<(self.size[1]-self.center[1])):
                
                x0 = round(self.center[0])
                y0 = round(self.center[1])
    
                self.data = self.data[y0-R:y0+R+1, x0-R:x0+R+1]
                self.notNan = self.notNan[y0-R:y0+R+1, x0-R:x0+R+1]
    
            # Centering the new cropped map and restoring offset
            self.recenter()
            self.xyOffset = xyOffset
            
        
        def createSurface(self,Rc,X,Y,zOffset=0,x0=0,y0=0,xTilt=0,yTilt=0,isPlot=False):
            '''
            Creating surface.
            
            Inputs: Rc, X, Y, zOffset=0, x0=0, y0=0, xTilt=0, yTilt=0, isPlot=False
            Rc      - Radius of curvature, and center of sphere on z-axis in case zOffset=0. [m]
            X, Y    - Matrices of x,y-values generated by 'X,Y = numpy.meshgrid(xVec,yVec)'. [m]
            zOffset - Surface center offset from 0. Positive value gives the surface center
                      positive z-coordinate. [self.scaling]
            x0,y0   - The center of the sphere in the xy-plane.
            x/yTilt - Tilt of x/y-axis [rad]. E.g. xTilt rotates around the y-axis.
            isPlot  - Plot or not [boolean]. Not recommended when used within optimization
                      algorithm.
    
            Returns: Z
            Z       - Matrix of surface height values. [self.scaling]
            
            Based on Simtools function 'FT_create_sphere_for_map.m' by Charlotte Bond.
            '''
            
            # Adjusting for tilts and offset
            Z = zOffset + (X*np.tan(xTilt) + Y*np.tan(yTilt))/self.scaling
            # Adjusting for spherical shape.
            if Rc !=0 and Rc is not None:
                Z = Z + (Rc - np.sign(Rc)*np.sqrt(Rc**2 - (X-x0)**2 - (Y-y0)**2))/self.scaling
            
            if isPlot:
                import pylab
                pylab.clf()
                fig = pylab.figure()
                axes = pylab.pcolormesh(X, Y, Z) #, vmin=zmin, vmax=zmax)
                cbar = fig.colorbar(axes)
                cbar.set_clim(Z.min(), Z.max())
                pylab.title('Z_min = ' + str(Z.min()) + '   Z_max = ' + str(Z.max()))
                pylab.show()
            return Z
        
        def createPolarGrid(self):
            '''
            Creating polar grid from the rectangular x and y coordinates in
            surfacemap. 
    
            Returns: rho, phi
            rho  - matrix with radial distances from centre of mirror map.
            phi  - matrix with polar angles.
            '''
            X,Y = np.meshgrid(self.x,self.y)
            phi = np.arctan2(Y,X)
            rho = np.sqrt(X**2 + Y**2)
            return rho, phi
    
        def preparePhaseMap(self, w=None, xyOffset=None, verbose=False):
            '''
            Prepares the phase mirror map for being used in Finesse. The map is cropped, and
            the curvature, offset and tilt are removed, in that order.
    
            Input: w, xyOffset, verbose
            w        - radius of Gaussian weights [m]. If specified, surfaces are fitted
                       to the mirror map when removing curvature and tilt. Otherise, the
                       whole mirror map is treated the same, thus covolutions with Zernike
                       polynomials are used instead.
            xyOffset - Beam center offset compared to the mirror center.
            verbose  - True makes a chatty method, while False only prints the end results.
    
            Output: amap
            amap     - Aperture map (absorption map), i.e. zeros within the map surface and
                       ones outside it.
                       
            Based on Charlotte Bond's Simtools function 'FT_prepare_phase_map_for_finesse.m'.
            '''
            if verbose:
                print('--------------------------------------------------')
                print('Preparing phase map for Finesse...')
                print('--------------------------------------------------')
            w_max = self.find_radius(method='min', unit='meters')
            if w is None:
                if verbose:
                    print(' No weighting used')
            elif w >= w_max:
                if verbose:
                    print(' Weighting radius ({:.3f} cm) larger than mirror radius ({:.3f} cm).'.format(w*100, w_max*100))
                    print(' Thus, no weighting used.')
                w = None
            elif verbose:
                print(' Gaussian weights used with radius: {:.2f} cm'.format(w*100.0))
            if verbose:
                print(' (rms, avg) = ({:.3e}, {:.3e}) nm'.format(self.rms(w),self.avg(w)))
                print('--------------------------------------------------')
            
                print(' Centering...')
            self.recenter()
            if verbose:
                print('  New center (x0, y0) = ({:.2f}, {:.2f})'.
                      format(self.center[0],self.center[1]))
            if verbose:
                print(' Cropping map...')
            before = self.size
            self.removePeriphery()
            if verbose:
                print('  Size (rows, cols): ({:d}, {:d}) ---> ({:d}, {:d})'.
                      format(before[1], before[0], self.size[1], self.size[0]))
                print('  New center (x0, y0) = ({:.2f}, {:.2f})'.
                      format(self.center[0],self.center[1]))
                print('  (rms, avg) = ({:.3e}, {:.3e}) nm'.
                      format(self.rms(w),self.avg(w)))
    
            # Radius of mirror in xy-plane.
            R = self.find_radius(unit='meters')
            
            if verbose:
                print(' Removing curvatures...')
            # --------------------------------------------------------
            if w is None:
                Rc, znm = self.remove_curvature(method='zernike', zModes = 'defocus')
                if verbose:
                    print('  Removed Z20 with amplitude A20 = {:.2f} nm'.format(znm['02'][2]))
                    print('  Equivalent Rc = {:.2f} m'.format(Rc))
            else:
                Rc, zOff, A20 = self.remove_curvature(method='sphere', w=w)
                # Adding A20 to zOff (calling it A0 now) because: Z(n=2,m=0) = A20*(r**2 - 1)
                A0 = zOff+A20
                if verbose:
                    print('  Removed Rc = {0:.2f} m'.format(Rc) )
                    print('  Equivalent Z(n=2,m=0) amplitude A20 = {0:.2f} nm'.format(A20))
    
            if verbose:
                print('  (rms, avg) = ({:.3e}, {:.3e}) nm'.format(self.rms(w),self.avg(w)))
                
                print(' Removing offset...')
            # --------------------------------------------------------
    
            zOff = self.removeOffset(w)
            if w is None:
                if verbose:
                    print('  Removed Z00 with amplitude A00 = {:.2f} nm'.format(zOff) )
            else:
                A0 = A0 + zOff
                if verbose:
                    print('  Removed z-offset (A00) = {0:.3f} nm'.format(zOff))
            if verbose:
                print('  (rms, avg) = ({:.3e}, {:.3e}) nm'.format(self.rms(w),self.avg(w)))
            
                print(' Removing tilts...')
            # --------------------------------------------------------
            if w is None:
                A1,xbeta,ybeta = self.rmTilt(method='zernike')
                if verbose:
                    for k in range(len(A1)):
                        print('  Removed Z1{:d} with amplitude A1{:d} = {:.2f} nm'.
                              format(2*k-1,2*k-1,A1[k]) )
                    print('  Equivalent tilts in radians: ')
                    print('   xbeta = {:.2e} rad'.format(xbeta))
                    print('   ybeta = {:.2e} rad'.format(ybeta))
            else:
                A1,xbeta,ybeta,zOff = self.rmTilt(method='fitSurf', w=w)
                A0 = A0 + zOff
                if verbose:
                    print('  Tilted surface removed:')
                    print('   xbeta    = {:.2e} rad'.format(xbeta))
                    print('   ybeta    = {:.2e} rad'.format(ybeta))
                    print('   z-offset = {:.2e} nm'.format(zOff))
                    print('  Equivalent Zernike amplitudes:')
                    print('   A(1,-1) = {:.2f} nm'.format(A1[0]))
                    print('   A(1, 1) = {:.2f} nm'.format(A1[1]))
            if verbose:
                print('  (rms, avg) = ({:.3e}, {:.3e}) nm'.format(self.rms(w),self.avg(w)))
    
            if w is not None:
                # Adding the total offset removed to zernikeRemoved (just for the record)
                self.zernikeRemoved = (0,0,A0)
                if verbose:
                    print(' Equivalent Z00 amplitude from accumulated z-offsets, A00 = {:.2f}'.format(A0))
                    
            if verbose:
                print(' Creating aperture map...')
                # --------------------------------------------------------
            amap = aperturemap(self.name, self.size, self.step_size,
                               self.find_radius(method='min',unit='meters'), self.center)
            if verbose:
                print('  Aperture map created.')
                                           
            if xyOffset is not None:
                if verbose:
                    print(' Offsetting mirror center in the xy-plane...')
                # --------------------------------------------------------
                self.center = (self.center[0] + xyOffset[0]/self.step_size[0],
                               self.center[1] + xyOffset[1]/self.step_size[1])
                self.xyOffset = xyOffset
                if verbose:
                    print('  New mirror center (x0, y0) = ({:.2f}, {:.2f})'.format(self.center[0],self.center[1]))
            else:
                self.xyOffset = (0,0)
                
            if verbose:
                print(' Writing maps to file...')
            # --------------------------------------------------------
            filename = self.name + '_finesse.txt'
            self.write_map(filename)
            if verbose:
                print('  Phase map written to file {:s}'.format(filename))
            filename = amap.name + '_aperture.txt'
            amap.write_map(filename)
            if verbose:
                print('  Aperture map written to file {:s}'.format(filename))
            if verbose:
                print(' Writing result information to file...')
                # --------------------------------------------------------
            filename = self.name + '_finesse_info.txt'
            self.writeMapPrepResults(filename,w)
            if verbose:
                print('  Result written to file {:s}'.format(filename))
                
                
            # Printing results to terminal
            print('--------------------------------------------------')
            print('Phase map prepared!' )
            print('Name: {:s}'.format(self.name))
            print('Type: {:s}'.format(self.type))
            print('--------------------------------------------------')
    
            if w is None:
                print('Weights:   None')
            else:
                print('Weights:   r     = {:.2f} cm'.format(w*100))
            print('Radius:    R     = {:.2f} cm'.format(self.find_radius(unit='meters')*100))
            print('Offset:    A00   = {:.2f} nm'.format(self.zernikeRemoved['00'][2]))
            print('Tilt:      A1-1  = {:.2f} nm'.format(self.zernikeRemoved['-11'][2]))
            print('           A11   = {:.2f} nm,  or'.format(self.zernikeRemoved['11'][2]))
            print('           xbeta = {:.2e} rad'.format(self.betaRemoved[0]))
            print('           ybeta = {:.2e} rad'.format(self.betaRemoved[1]))
            print('Curvature: A20   = {:.2f} nm,  or'.format(self.zernikeRemoved['02'][2]))
            print('           Rc    = {:.2f} m'.format(self.Rc))
            print('xy-offset: x0    = {:.2f} mm'.format(self.xyOffset[0]*1000))
            print('           y0    = {:.2f} mm'.format(self.xyOffset[1]*1000))
            print('Stats:     rms   = {:.3e} nm'.format(self.rms(w)))
            print('           avg   = {:.3e} nm'.format(self.avg(w)))
            print('--------------------------------------------------')
            print()
    
            return amap
        
        def writeMapPrepResults(self, filename, w):
            '''
            Writing result information to file. The same info that is written to the
            terminal.
            '''
            import time
            with open(filename,'w') as mapfile:
                mapfile.write('---------------------------------------\n')
                mapfile.write('Map:  {:s}\n'.format(self.name))
                mapfile.write('Date:  {:s}\n'.format(time.strftime("%d/%m/%Y  %H:%M:%S")))
                mapfile.write('---------------------------------------\n')
                if w is None:
                    mapfile.write('Weights:   None')
                else:
                    mapfile.write('Weights:   r     = {:.2f} cm'.format(w*100))
                mapfile.write('Diameter:  D     = {:.2f} cm\n'.format(self.find_radius(unit='meters')*200.0))
                mapfile.write('Offset:    A00   = {:.2f} nm\n'.format(self.zernikeRemoved['00'][2]))
                mapfile.write('Tilt:      A1-1  = {:.2f} nm\n'.format(self.zernikeRemoved['-11'][2]))
                mapfile.write('           A11   = {:.2f} nm,  or\n'.format(self.zernikeRemoved['11'][2]))
                mapfile.write('           xbeta = {:.2e} rad\n'.format(self.betaRemoved[0]))
                mapfile.write('           ybeta = {:.2e} rad\n'.format(self.betaRemoved[1]))
                mapfile.write('Curvature: A20   = {:.2f} nm,  or\n'.format(self.zernikeRemoved['02'][2]))
                mapfile.write('           Rc    = {:.2f} m\n'.format(self.Rc))
                mapfile.write('xy-offset: x0    = {:.2f} mm\n'.format(self.xyOffset[0]*1000))
                mapfile.write('           y0    = {:.2f} mm\n'.format(self.xyOffset[1]*1000))
                
                # mapfile.write('Offset (Z00):  {:.2f}'.format(self.zernikeRemoved['00'][2]))
              
            
        def removeOffset(self, r=None):
            '''
            Based on Simtools function 'FT_remove_offset_from_mirror_map.m' by Charlotte Bond.
            '''
            if r is None:
                offset = self.rmZernike(0,0)
            else:
                rho = self.createPolarGrid()[0]
                inside = rho<r
                inside_notNan = np.zeros(inside.shape,dtype=bool)
                inside_notNan[inside] = self.notNan[inside]
    
                offset = self.data[inside_notNan].sum()/inside_notNan.sum()
                self.data[self.notNan] = self.data[self.notNan]-offset
            
            return offset
    
        def rmZernike(self, n, m='all'):
            '''
            Removes Zernike polynomials from mirror map.
    
            Input: n, m
            n  - Radial degree of Zernike polynomial.
            m  - Azimuthal degree of Zernike polynomial. If set to 'all', all polynomials
                 of radial degree n is removed. Can also be an integer, or a list
                 specifying the Azimuthal degrees that will be removed. Must be consistent
                 with the chosen radial degree.
    
            Returns: A
            A  - List of amplitudes of the removed Zernike polynomials.
            '''
            if m=='all':
                m = range(-n,n+1,2)
                
            R = self.find_radius(unit='meters')
            A = self.zernikeConvol(n,m)
            rho, phi= self.createPolarGrid()
            rho = rho/R
    
            if isinstance(m,list):
                for k in range(len(m)):
                    Z = zernike(m[k], n, rho, phi)
                    self.data[self.notNan] = self.data[self.notNan]-A[k]*Z[self.notNan]
                    self.zernikeRemoved = (m[k], n, A[k])
            else:
                Z = zernike(m, n, rho, phi)
                self.data[self.notNan] = self.data[self.notNan]-A*Z[self.notNan]
                self.zernikeRemoved = (m, n, A)
            return A
            
    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 = []
            self.weighting = None
            
        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 type(self):
            hasR = False
            hasT = False
            
            _type = ""
            
            for m in self.__maps:
                if "reflection" in m.type: hasR = True
                
                if "transmission" in m.type: hasT = True
                
                if "both" in m.type:
                    hasR = True
                    hasT = True
            
            if hasR and not hasT: _type += "reflection "
            elif hasR and not hasT: _type += "transmission "
            elif hasR and hasT: _type += "both "
            
            return _type
            
        @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(0, self.size[0])) - self.center[0])
            
        @property
        def y(self):
            return self.step_size[1] * (np.array(range(0, self.size[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) - (np.array(self.size)-1)/2.0)
        
        @property
        def ROMWeights(self):
            return self._rom_weights
        
        def z_xy(self, wavelength=1064e-9, direction="reflection_front", 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)
                
            if self.weighting is None:
                return z_xy
            else:
                return z_xy * self.weighting
            
        def generateROMWeights(self, EIxFilename, EIyFilename=None, verbose=False, interpolate=False, newtonCotesOrder=8, nr1=1, nr2=1):
            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)
            
            w_refl_front, w_refl_back, w_tran_front, w_tran_back = (None, None, None, None)
            
            if "reflection" in self.type or "both" in self.type:
                w_refl_front = makeWeightsNew(self, EIxFilename, EIyFilename,
                                          verbose=verbose, newtonCotesOrderMapWeight=newtonCotesOrder,
                                          direction="reflection_front")
                
                w_refl_front.nr1 = nr1
                w_refl_front.nr2 = nr2
                
                w_refl_back = makeWeightsNew(self, EIxFilename, EIyFilename,
                                          verbose=verbose, newtonCotesOrderMapWeight=newtonCotesOrder,
                                          direction="reflection_back")
                
                w_refl_back.nr1 = nr1
                w_refl_back.nr2 = nr2
    
            if "transmission" in self.type or "both" in self.type:                                      
                w_tran_front = makeWeightsNew(self, EIxFilename, EIyFilename,
                                          verbose=verbose, newtonCotesOrderMapWeight=newtonCotesOrder,
                                          direction="transmission_front")
    
                w_refl_front.nr1 = nr1
                w_refl_front.nr2 = nr2
                                                
                w_tran_back  = makeWeightsNew(self, EIxFilename, EIyFilename,
                                          verbose=verbose, newtonCotesOrderMapWeight=newtonCotesOrder,
                                          direction="transmission_back")
                
                w_refl_back.nr1 = nr1
                w_refl_back.nr2 = nr2
                
            self._rom_weights = MirrorROQWeights(w_refl_front, w_refl_back, w_tran_front, w_tran_back)
            
            return self._rom_weights
    
        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, center=None):
            if center is None:
                center = (np.array(size)+1)/2.0
            surfacemap.__init__(self, name, "absorption both", size, center, step_size, 1)
            self.R = R
            
        @property
        def R(self):
            return self.__R
        
        @R.setter
        def R(self, value):
            self.__R = value
        
            xx, yy = np.meshgrid(self.x, self.y)
            
            radius = np.sqrt(xx**2 + yy**2)
            
            self.data = np.zeros(self.size[::-1])
            self.data[radius > self.R] = 1.0
            
            
    class curvedmap(surfacemap):
        
        def __init__(self, name, size, step_size, Rc):
            surfacemap.__init__(self, name, "phase reflection", size, (np.array(size)+1)/2.0, step_size, 1e-6)
            self.Rc = Rc
            
        @property
        def Rc(self):
            return self.__Rc
        
        @Rc.setter
        def Rc(self, value):
            self.__Rc = value
        
            xx, yy = np.meshgrid(self.x, self.y)
            
            Rsq = xx**2 + yy**2
            self.data = (self.Rc - math.copysign(1.0, self.Rc) * np.sqrt(self.Rc**2 - Rsq))/ self.scaling
    
    class tiltmap(surfacemap):
        """
        To create a tiltmap, plot it and write it to a file to use with Finesse:
            
            tilts = (1e-6, 1e-8) # tilt in (x, y) radians\
            dx = 1e-4
            L = 0.2
            N = L/dx
            
            tmap = tiltmap("tilt", (N, N), (dx,dx), tilts)
            tmap.plot()
            tmap.write_map("mytilt.map")
        """
        
        def __init__(self, name, size, step_size, tilt):
            surfacemap.__init__(self, name, "phase reflection", size, (np.array(size)+1)/2.0, step_size, 1e-9)
            self.tilt = tilt
            
        @property
        def tilt(self):
            return self.__tilt
        
        @tilt.setter
        def tilt(self, value):
            self.__tilt = value
            
            xx, yy = np.meshgrid(self.x, self.y)
            
            self.data = (yy * self.tilt[1] + xx * self.tilt[0])/self.scaling
            
    
    class zernikemap(surfacemap):
    	def __init__(self, name, size, step_size, radius, scaling=1e-9):
    		surfacemap.__init__(self, name, "phase reflection", size, (np.array(size)+1)/2.0, step_size, scaling)
    		self.__zernikes = {}
    		self.radius = radius
    		
    	@property
    	def radius(self): return self.__radius
    
    	@radius.setter
    	def radius(self, value, update=True):
    		self.__radius = float(value)
    		if update: self.update_data()
    
    	def setZernike(self, m, n, amplitude, update=True):
    		self.__zernikes["%i%i" % (m, n)] = (m,n,amplitude)
    		if update: self.update_data()
    
    	def update_data(self):
    		X,Y = np.meshgrid(self.x, self.y)
    		R = np.sqrt(X**2 + Y**2)
    		PHI = np.arctan2(Y, X)
    
    		data = np.zeros(np.shape(R))
    
    		for i in self.__zernikes.items():
    			data += i[1][2] * zernike(i[1][0], i[1][1], R/self.radius, PHI)
    
    		self.data = data
    	
    			
    
      
    def read_map(filename, mapFormat='finesse', scaling=1.0e-9, mapType='phase', field='both'):
        '''
        Reads surface map files and returns a surfacemap object.
        
        filename  - name of surface map file.
        mapFormat - 'finesse', 'ligo', 'zygo'. Currently only for ascii formats.
        scaling   - scaling of surface height of the mirror map [m].
        '''
        
        # Reads finesse mirror maps.
        if mapFormat.lower() == 'finesse':
            
            g = lambda x: float(x)
            with open(filename, 'r') as f:
            
                f.readline()
                name = f.readline().split(':')[1].strip()
                maptype = f.readline().split(':')[1].strip()
                size = tuple(map(g, f.readline().split(':')[1].strip().split()))
                center = tuple(map(g, f.readline().split(':')[1].strip().split()))
                step = tuple(map(g, f.readline().split(':')[1].strip().split()))
                scaling = float(f.readline().split(':')[1].strip())
            
            data = np.loadtxt(filename, dtype=np.float64,ndmin=2,comments='%')    
    
            return surfacemap(name, maptype, size, center, step, scaling, data)
            
        elif mapFormat.lower() == 'ligo' or mapFormat.lower() == 'zygo':
            '''
            Reads zygo and ligo mirror maps. Based on Charlotte Bond's Simtools
            functions 'FT_read_zygo_map.m' and 'FT_read_ligo_map.m'.
            '''
            maptype = ' '.join([mapType, field])
            if mapFormat == 'ligo':
                isLigo = True
                # Unused, but need to set it to call the function.
                isAscii = False
                # Remove '_asc.dat' for output name
                name = filename.split('_')
                name = '_'.join(name[:-1])
            else:
                isLigo = False
                tmp = filename.split('.')
                fileFormat = tmp[-1].strip()
                name = '.'.join(tmp[:-1])
                if fileFormat == 'asc':
                    isAscii = True
                else:
                    isAscii = False
                    
            # Reading the zygo/ligo surface map file. Note that the intensity data
            # iData is available here too when ascii zygo files are read, but at
            # the moment it's not used. 
            out = readZygoLigoMaps(filename, isLigo, isAscii)
            if len(out) == 3:
                hData, data, iData = out
            else:
                hData, data = out
            
            # Matrix with True where data element is NaN.
            isNan = np.isnan(data)
            # Matrix with True where data element is not NaN.
            notNan = isNan==False
            # Setting NaN-values to zero
            data[isNan] = 0
            # Scaling to the scaling chosesen as input.
            data[notNan] = (hData['hScale']/scaling)*data[notNan]
            smap = surfacemap(name, maptype, hData['size'], hData['center'],
                              hData['stepSize'], scaling, data, notNan)
            smap.recenter()
            
            return smap
    
            
        if mapFormat == 'metroPro':
            # Remove '.dat' for output name
            name = filename.split('.')[0]
            # Reading file
            data, hData = readMetroProData(filename)
            # xy-stepsize in meters
            stepSize = (hData['cameraRes'], hData['cameraRes'])
            # Temporary center
            center = ((len(data[0,:])-1)/2.0,(len(data[:,0])-1)/2.0)
            # Keeping track of nan's
            isNan = np.isnan(data)
            notNan = isNan==False
            data[isNan] = 0
            # Type of map
            maptype = ' '.join([mapType, field])
            # Unnecessary size, remove from class later
            size = data.shape[::-1]
            # Creating surfacemap-instance
            smap = surfacemap(name, maptype, size, center, stepSize, 1.0, data, notNan)
            smap.rescaleData(scaling)
            smap.recenter()
            return smap
    
    
    
    
    def readZygoLigoMaps(filename, isLigo=False, isAscii=True):
        '''
        Converts raw zygo and ligo mirror maps to the finesse
        format. Based on the matlab scripts 'FT_read_zygo_map.m' and
        'FT_read_ligo_map.m'.
        '''
        g = lambda x: float(x)
        hData = {}
        with open(filename, 'r') as f:
            # Reading header of LIGO-map
            # ------------------------------------------------------
            # Skip first two lines
            for k in range(2):
                f.readline()
            # If zygo-file, and ascii format, there is intensity
            # data. Though, the Ligo files I have seen are also
            # zygo-files, so maybe we should extract this data
            # from these too?
            line = f.readline()
            if not isLigo and isAscii:
                iCols = float(line.split()[2])
                iRows = float(line.split()[3])
    
            line = f.readline().split()
            
            if isLigo:
                y0 = float(line[0])
                x0 = float(line[1])
                rows = float(line[2])
                cols = float(line[3])
            else:
                y0 = float(line[1])
                x0 = float(line[0])
                rows = float(line[3])
                cols = float(line[2])
    
            # Skipping three lines
            for k in range(3):
                f.readline()
            line = f.readline().split()
    
            # Scaling factors
            # ----------------------------------------------
            # Interfeometric scaling factor
            S = float(line[1])
            # Wavelength
            lam = float(line[2])
            # Obliquity factor
            O = float(line[4])
            # ----------------------------------------------
            # Physical step size in metres
            if line[6] != 0:
                xstep = float(line[6])
                ystep = float(line[6])
            else:
                xstep = 1.0
                ystep = 1.0
    
            # Skipping two lines
            for k in range(2):
                f.readline()
            line = f.readline().split()
    
            # Resolution of phase data points, 1 or 0.
            phaseRes = float(line[0])
            if phaseRes == 0:
                R = 4096
            elif phaseRes == 1:
                R = 32768
    
            if not isLigo and not isAscii:
                # zygo .xyz files give phase data in microns.
                hData['hScale'] = 1.0e-6
                
            else:
                # zygo .asc and ligo-files give phase data in
                # internal units. To convert to m use hScale
                # factor.
                hData['hScale'] = S*O*lam/R
    
            # Skipping three lines
            for k in range(3):
                f.readline()
                
            if not isLigo:
                if not isAscii:
                    # For the zygo .xyz-format.
                    k = 0
                    data = np.zeros(rows*cols)
                    totRuns = cols*rows
                
                    # Read data
                    while k < cols*rows:
                        line = f.readline().split()
                        # Setting 'no data' to NaN
                        if len(line)==4:
                            data[k]=np.nan
                        # Reading value
                        elif len(line)==3:
                            data[k] = float(line[2])
                        # Otherwise stop reading
                        else:
                            k = cols*rows
                        k+=1
                else:
                    # Skipping one line
                    f.readline()
            
                    # Reading intensity data
                    iData = np.array([])
                    line = f.readline().split()
                    while line[0] != '#':
                        iData = np.append(iData, map(g,line))
                        line = f.readline().split()
                    # Reshaping intensity data
                    iData = iData.reshape(iRows, iCols).transpose()
                    iData = np.rot90(iData)
            else:
                # Skipping one line
                f.readline()
                # Skipping lines until '#' is found.
                while f.readline()[0] != '#':
                    pass
    
            if isLigo or (not isLigo and isAscii):
                # Reading phase data
                # ----------------------------------------------
                # Array with the data
                data = np.array([])
                # Reading data until next '#' is reached.
                line = f.readline().split()
                while line[0] != '#':
                    data = np.append(data, map(g,line))
                    line = f.readline().split()
                # ----------------------------------------------
    
    
        if isLigo:
            # Setting all the points outside of the mirror
            # surface to NaN. These are given a large number
            # in the file. 
            data[data == data[0]] = np.nan
    
            # Reshaping into rows and columns
            data = data.reshape(cols,rows).transpose()
            # Pretty sure that the lines below can be done in
            # less operations, but it's quick as it is.
            # ----------------------------------------------
            # Flipping right and left
            data = np.fliplr(data)
            # Rotating 90 degrees clockwise 
            data = np.rot90(data,-1)
            # Flipping right and left
            data = np.fliplr(data)
            # ----------------------------------------------
        else:
            if isAscii:
                # Setting all the points outside of the mirror
                # surface to NaN. These are given a large number
                # in the file. 
                data[data >= 2147483640] = np.nan
            # Reshaping into rows and columns.
            data = data.reshape(rows,cols).transpose()
            # Rotating to make (0,0) be in bottom left
            # corner. 
            data = np.rot90(data)
    
        hData['size'] = data.shape[::-1]
    
        # Center according to file
        hData['center'] = tuple([x0,y0])
        # xy-stepsize
        hData['stepSize'] = tuple([xstep,ystep])
    
        if not isLigo and isAscii:
            return hData, data, iData
        else:
            return hData, data
    
    
    
    def readMetroProData(filename):
        '''
        Reading the metroPro binary data files.
    
        Translated from Hiro Yamamoto's 'LoadMetroProData.m'.
        '''
        f = BinaryReader(filename)
        # Read header
        hData = readHeaderMP(f)
        if hData['format'] < 0:
            print('Error: Format unknown to readMetroProData()\nfilename: {:s}'.format(filename))
            return 0
        # Read phase map data
        # Skipping header and intensity data
        f.seek(hData['size']+hData['intNBytes'])
        # Reading data
        dat = f.read('int32',size=hData['Nx']*hData['Ny'])
        # Marking unmeasured data as NaN
        dat[dat >= hData['invalid']] = np.nan
        # Scale data to meters
        dat = dat*hData['convFactor']
        # Reshaping into Nx * Ny matrix
        dat = dat.reshape(hData['Ny'], hData['Nx'])
        # Flipping up/down, i.e., change direction of y-axis.
        dat = dat[::-1,:]
        # Auxiliary data to return
        dxy = hData['cameraRes']
    
        # Unnecessary
        #x1 = dxy * np.arange( -(len(dat[0,:])-1)/2, (len(dat[0,:])-1)/2 + 1)
        #y1 = dxy * np.arange( -(len(dat[:,0])-1)/2, (len(dat[:,1])-1)/2 + 1)
    
        return dat, hData #, x1, y1
    
    
    def readHeaderMP(f):
        '''
        Reads header of the metroPro format files.
        
        Translated from the readerHeader() function within the
        'LoadMetroProData.m' function written by Hiro Yamamoto.
        '''
        hData = {}
        hData['magicNum'] = f.read('uint32')
        hData['format'] = f.read('int16')
        hData['size'] = f.read('int32')
        # Check if the magic string and format are known ones.
        if not (hData['format'] >= 1 and hData['format']<=3 and
                hData['magicNum']-hData['format'] == int('881B036E',16)):
            hData['format'] = -1
        # Read necessary data
        hData['invalid'] = int('7FFFFFF8',16)
        f.seek(60)
        # Intensity data, which we will skip over.
        hData['intNBytes'] = f.read('int32')
        # Top-left coordinates, which are useless.
        hData['X0'] = f.read('int16')
        hData['Y0'] = f.read('int16')
        # Number of data points alng x and y
        hData['Nx'] = f.read('int16')
        hData['Ny'] = f.read('int16')
        # Total data, 4*Nx*Ny
        hData['phaNBytes'] = f.read('int32')
        f.seek(218)
        # Scale factor determined by phase resolution tag
        phaseResTag = f.read('int16')
        if phaseResTag == 0:
            phaseResVal = 4096
        elif phaseResTag == 1:
            phaseResVal = 32768
        elif phaseResTag == 2:
            phaseResVal = 131072
        else:
            phaseResVal = 0
        f.seek(164)
        intf_scale_factor = f.read('float')
        hData['waveLength'] = f.read('float')
        f.seek(176)
        obliquity_factor = f.read('float')
        # Eq. in p12-6 in MetroPro Reference Guide
        hData['convFactor'] = intf_scale_factor*obliquity_factor*\
                              hData['waveLength']/phaseResVal
        # Bin size of each measurement
        f.seek(184)
        hData['cameraRes'] = f.read('float')
        
        return hData
    
    
    
    
    
    class BinaryReaderEOFException(Exception):
        def __init__(self):
            pass
        def __str__(self):
            return 'Not enough bytes in file to satisfy read request'
    
    class BinaryReader:
        # Map well-known type names into struct format characters.
        typeNames = {
            'int8'   :'b',
            'uint8'  :'B',
            'int16'  :'h',
            'uint16' :'H',
            'int32'  :'i',
            'uint32' :'I',
            'int64'  :'q',
            'uint64' :'Q',
            'float'  :'f',
            'double' :'d',
            'char'   :'s'}
    
        def __init__(self, fileName):
            self.file = open(fileName, 'rb')
            
        def read(self, typeName, size=None):
            import struct
            
            typeFormat = BinaryReader.typeNames[typeName.lower()]
            typeFormat = '>'+typeFormat
            typeSize = struct.calcsize(typeFormat)
            if size is None:
                value = self.file.read(typeSize)
                if typeSize != len(value):
                    raise BinaryReaderEOFException
                unpacked = struct.unpack(typeFormat, value)[0]
            else:
                value = self.file.read(size*typeSize)
                if size*typeSize != len(value):
                    raise BinaryReaderEOFException
                unpacked = np.zeros(size)
                for k in range(size):
                    i = k*typeSize
                    unpacked[k] = struct.unpack(typeFormat,value[i:i+typeSize])[0]
            return unpacked
    
        def seek(self, offset, refPos=0):
            '''
            offset in bytes and refPos gives reference position, where 0
            means origin of the file, 1 uses current position and 2 uses
            the end of the file.
            '''
            self.file.seek(offset, refPos)
        
        
        def __del__(self):
            self.file.close()
    
            
    # TODO: Add options for reading virgo maps, and .xyz zygo
        # maps (need .xys file for this). Binary ligo-maps?
        # The intensity data is not used to anything here. Remove
        # or add to pykat?
        
        
        
        
    # TODO: Recreate functions from Simtools:, List taken from: ligo_maps/FT_convert_ligo_map_for_finesse.m
    # map=FT_recenter_mirror_map(map);
    # [map2,A2,Rc_out]=FT_remove_zernike_curvatures_from_map(map,Rc_in);
    # [map2,Rc_out]=FT_remove_curvature_from_mirror_map(map,Rc_in,w, display_style);
    # [map2,offset]=FT_remove_offset_from_mirror_map(map2,1e-2);
    # [map3,x_tilt,y_tilt,offset2]=FT_remove_piston_from_mirror_map(map2,w, display_style);
    # map3=FT_invert_mirror_map(map3, invert);