diff --git a/pykat/optics/maps.py b/pykat/optics/maps.py
index c792d0a29574f9ac42aa9e1fa67e2c36fe77fc8c..2afb86b0550f6aaf10606b30db7689f16daecd86 100644
--- a/pykat/optics/maps.py
+++ b/pykat/optics/maps.py
@@ -283,37 +283,58 @@ class surfacemap(object):
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]
-
+
+ # min and max of z-values
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 = 100*self.x
+ yRange = 100*self.y
+
+ # This line is added by DT to be able to plot
+ # rectangular matrices. Effectively, I swapped the
+ # x/y-axes. Preferrably, this should be corrected above
+ # instead, but since I'm not completely sure of how the
+ # coordinate system of these maps look I'll wait with
+ # that. Here, I assume that row 0 of the matrix should
+ # be plotted with y = Y[0], and that column 0 should be
+ # plotted with x = X[0]. To be fully correct, I should
+ # add one column and one row so that each matrix value
+ # is plotted within the correct rectangle.
+ # ------------------------------------------------------
+ xRange, yRange = np.meshgrid(yRange,xRange)
+ # ------------------------------------------------------
+
fig = pylab.figure()
- axes = pylab.pcolormesh(xrange, yrange, self.data, vmin=zmin, vmax=zmax)
+ axes = pylab.pcolormesh(xRange, yRange, self.data,
+ vmin=zmin, vmax=zmax)
+
pylab.xlabel('x [cm]')
pylab.ylabel('y [cm]')
@@ -671,24 +692,186 @@ class zernikemap(surfacemap):
-def read_map(filename):
- with open(filename, 'r') as f:
-
- f.readline()
- name = f.readline().split(':')[1].strip()
- maptype = f.readline().split(':')[1].strip()
- size = tuple(map(lambda x: float(x), f.readline().split(':')[1].strip().split()))
- center = tuple(map(lambda x: float(x), f.readline().split(':')[1].strip().split()))
- step = tuple(map(lambda x: float(x), f.readline().split(':')[1].strip().split()))
- scaling = float(f.readline().split(':')[1].strip())
-
+def read_map(filename, mapFormat='finesse'):
+ # Function turning input x into float.
+ g = lambda x: float(x)
+ if mapFormat == 'finesse':
+
+ 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='%')
+
+ # Converts raw ligo mirror map to the finesse format. Based
+ # on translation of the matlab script 'FT_read_ligo_map.m'.
+ #
+ elif mapFormat == 'ligo':
+ # Unknown usage
+ baseid = 'read_ligo_map'
+ # Remove '_asc.dat' for output name
+ name = filename.split('_')
+ name = '_'.join(name[:-1])
+
+ # Unknowns (why are these values hard coded here?)
+ # ------------------------------------------------------
+ # Standard maps have type 'phase' (they store surface
+ # heights)
+ maptype = 0
+ # Both (reflected and transmitted) light fields are
+ # affected
+ field = 0
+ # Measurements in nanometers
+ scaling = 1.0e-9
+ # ------------------------------------------------------
+
+ # Reading header of LIGO-map (Zygo file? Says Zygo in
+ # header...)
+ # ------------------------------------------------------
+ with open(filename, 'r') as f:
+ # Skip first three lines
+ for k in range(3):
+ f.readline()
+ line = f.readline().split()
+
+ # Unknown
+ # ----------------------------------------------
+ y0 = float(line[0])
+ x0 = float(line[1])
+ rows = float(line[2])
+ cols = float(line[3])
+ # ----------------------------------------------
+
+ # Skipping three lines
+ for k in range(3):
+ f.readline()
+ line = f.readline().split()
+
+ # Unknown (Scaling factors)
+ # ----------------------------------------------
+ # Interfeometric scaling factor (?)
+ S = float(line[1])
+ # wavelength (of what?)
+ lam = float(line[2])
+ # Obliquity factor (?)
+ O = float(line[4])
+ # Guessing physical step size
+ xstep = float(line[6])
+ ystep = float(line[6])
+ # ----------------------------------------------
+
+ # Skipping two lines
+ for k in range(2):
+ f.readline()
+ line = f.readline().split()
+
+ # Unknown
+ # Resolution of phase data points, 1 or 0.
+ phaseRes = float(line[0])
+ if phaseRes == 0:
+ R = 4096
+ elif phaseRes == 1:
+ R = 32768
+ else:
+ print('Error, invalid phaseRes')
+
+ # Unknown
+ hScale = S*O*lam/R
+
+ # Skipping four lines
+ for k in range(4):
+ f.readline()
+ # Skipping lines until '#' is found.
+ while f.readline()[0] != '#':
+ pass
+
+ # Reading the data
+ # ----------------------------------------------
+ # Function converting string to float
+ g = lambda x: float(x)
+ # 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()
+ # ----------------------------------------------
+
+ # 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
+ # more efficient, 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)
+ # ----------------------------------------------
+
+ # Scaling to nanometer (change this to a user
+ # defined value?) Still don't know where
+ # 'hScale' really comes from.
+ data = (hScale/scaling)*data
+ size = data.shape
+
+ if maptype == 0:
+ mType = 'phase'
+ else:
+ mType = 'Unknown'
+ if field == 0:
+ fType = 'both'
+ else:
+ fType = 'unknown'
+
+ maptype = ' '.join([mType, fType])
+
+ print(maptype)
+ # Wrong! fix by creating recenter method.
+ center = tuple([x0,y0])
+ step = tuple([xstep,ystep])
+
+ # Simple re-centering of mirror, translated from
+ # 'FT_recenter_mirror_map.m'
+ # -------------------------------------------------
+ # Matrix with ones where data element is not NaN.
+ isNan = np.isnan(data)
+ notNan = isNan==False
+ # Row and column indices with non-NaN elements
+ rIndx, cIndx = notNan.nonzero()
+ # Finding centres
+ x0 = float(cIndx.sum())/len(cIndx)
+ y0 = float(rIndx.sum())/len(rIndx)
+ center = tuple([x0,y0])
+ # -------------------------------------------------
+ # Changing NaN to zeros. Just to be able to plot the
+ # map with surfacemap.plot().
+ data[isNan] = 0
- data = np.loadtxt(filename, dtype=np.float64,ndmin=2,comments='%')
- return surfacemap(name,maptype,size,center,step,scaling,data)
+ # TODO: Add options for reading zygo and virgo maps too.
+
+ return surfacemap(name, maptype, size, center, step,
+ scaling, data)
+
+
# 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);
@@ -696,3 +879,6 @@ def read_map(filename):
# [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);
+
+# Understand the internal coordinate system of the
+# maps/matrices.