Began versioning
authorW. Trevor King <wking@drexel.edu>
Wed, 22 Oct 2008 19:05:38 +0000 (15:05 -0400)
committerW. Trevor King <wking@drexel.edu>
Wed, 22 Oct 2008 19:05:38 +0000 (15:05 -0400)
README [new file with mode: 0644]
gnuplot_getSurfPos.py [new file with mode: 0755]
test_vib_analyze.py [new file with mode: 0644]
x_piezo.py [new file with mode: 0644]
z_piezo.py [new file with mode: 0644]
z_piezo_utils.py [new file with mode: 0644]

diff --git a/README b/README
new file mode 100644 (file)
index 0000000..7ce2cc9
--- /dev/null
+++ b/README
@@ -0,0 +1,4 @@
+python calibrate_raw.py ../comedi_simult_AIO/20080128/far_45.raw
+
+veeco_header is the beginning of a saved Nanoscope image file.
+It lists the parameters used by the Nanoscope during imaging.
diff --git a/gnuplot_getSurfPos.py b/gnuplot_getSurfPos.py
new file mode 100755 (executable)
index 0000000..7983b92
--- /dev/null
@@ -0,0 +1,44 @@
+#!/usr/bin/python
+#
+# Usage: gnu_z_peizo_getSurfPos <filename>
+#
+# Where <filename> contains the data in ascii pairs
+# followed by a blank line
+# followed by lines: 
+#  #Guess params: [A, B, C, D]
+#  #Fit params: [A, B, C, D]
+
+import os
+import sys
+
+def get_params_from_line (line) :
+    "'#BlaBlaBla: [123, 456, 7.89, -10]' -> (123, 456, 7.89, -10)"
+    array_and_tailing_bracket = line.split('[')[1]
+    array = array_and_tailing_bracket.split(']')[0]
+    nums = array.split(',')
+    return (nums[0], nums[1], nums[2], nums[3])
+
+def get_line_starting_with(filename, start_string) :
+    ln = len(start_string)
+    for line in file(filename).readlines() :
+        if line[0:ln] == start_string :
+            return line
+
+filename = sys.argv[1]
+x_data_col = 1 # gnuplot columns start with 1
+y_data_col = 3
+
+start_line = get_line_starting_with(filename, "#Guess params:")
+s0,s1,s2,s3 = get_params_from_line(start_line)
+fit_line = get_line_starting_with(filename, "#Fit params:")
+f0,f1,f2,f3 = get_params_from_line(fit_line)
+
+print ""
+print "start(x) =  x < %s  ?  (%s)+(%s)*x  :  (%s)+(%s)*(%s)+(%s)*(x-(%s))" \
+                  % (s2,    s0,s1,      s0,s1,s2,s3,   s2)
+print "set arrow from %s, graph 0 to %s, graph 1 nohead" % (s2, s2)
+print "fit(x)   =  x < %s  ?  (%s)+(%s)*x  :  (%s)+(%s)*(%s)+(%s)*(x-(%s))" \
+                  % (f2,    f0,f1,      f0,f1,f2,f3,   f2)
+print "set arrow from %s, graph 0 to %s, graph 1 nohead" % (f2, f2)
+print "plot '%s' using %d:%d, start(x), fit(x)" % (filename, x_data_col, y_data_col)
+
diff --git a/test_vib_analyze.py b/test_vib_analyze.py
new file mode 100644 (file)
index 0000000..b9c7a20
--- /dev/null
@@ -0,0 +1,8 @@
+import sys
+import calibrate_cantilever
+
+tweak_file = sys.argv[1]
+
+calibrate_cantilever.PYLAB_INTERACTIVE = False
+calibrate_cantilever.vib_load_analyze_tweaked(tweak_file, plotVerbose=True)
+calibrate_cantilever._final_flush_plot()
diff --git a/x_piezo.py b/x_piezo.py
new file mode 100644 (file)
index 0000000..6e87591
--- /dev/null
@@ -0,0 +1,90 @@
+import comedi_single_aio
+#import sngAO
+
+class x_piezoError (Exception) :
+    pass
+
+class x_piezo :
+    def __init__(self, xp_chan=1) :
+        self.verbose = False
+        self.xpSensitivity = 34.0 # nm/Volt
+        self.gain = 20 # Vpiezo / Voutput
+
+        self.AO = comedi_single_aio.ao_obj(chan=(xp_chan,))
+        self.AO.close()
+       #self.AO = sngAO.sngAOobj([xp_chan], -10, 10)
+
+       self.curOut = [self.nm2out(0)]*1
+       self.xpMax = self.nm2out(2000) # limits at 2 microns
+       self.xpMin = self.nm2out(-2000) # and -2 microns
+
+       self.wanderStep = self.nm2out(5) - self.nm2out(0)
+       # define default value changed callbacks
+       self.setExternalXPiezo = None
+    def curPos(self) :
+        return self.curOut[0]
+    def nm2out(self, nm) :
+        # nm / (nm/Vp * Vp/Vout) = Vout
+        # where Vout is Volts output by the DAQ card
+        # and Vp is Volts applied across the piezo
+        return self.phys_to_out(nm / (self.xpSensitivity * self.gain))
+    def out2nm(self, out) :
+        return self.out_to_phys(out) * self.xpSensitivity * self.gain
+    def out_to_phys(self, output) :
+        return self.AO.comedi_to_phys(0, output)
+    def phys_to_out(self, physical) :
+        return self.AO.phys_to_comedi(0, physical)
+    def pCurVals(self) :
+        print "X piezo output  : %6d" % self.curPos()
+    def jumpToPos(self, pos) :
+        self._check_range(pos)
+        if self.verbose :
+            print "Jump X piezo to %d (%g nm)" % (pos, self.out2nm(pos))
+       self.curOut[0] = pos
+        self.AO.open()
+       #self.AO.reserve()
+       self.AO.write(self.curOut)
+        self.AO.close()
+       #self.AO.unreserve()
+       if self.setExternalXPiezo != None :
+           self.setExternalXPiezo(self.out2nm(self.curPos()))
+       return self.curPos()
+    def wander(self) :
+               """
+       steps in one direction until it hits the boundary,
+       then reverses and walks back the other way.
+       """
+       newPos = self.curPos() + self.wanderStep
+        if self.verbose :
+            print "wander from %d to %d" % (self.curPos(), newPos)
+       try :
+           self._check_range(newPos)
+        except Exception :
+           # reached one limit, so turn around
+           self.wanderStep = -self.wanderStep
+           newPos += 2*self.wanderStep
+       return self.jumpToPos(newPos)           
+    def _check_range(self, pos) :
+        if pos > self.xpMax :
+           raise x_piezoError, "X piezo pos = %d > %d = max" % (pos, self.xpMax)
+       if pos < self.xpMin :
+           raise x_piezoError, "X piezo pos = %d < %d = min" % (pos, self.xpMin)
+
+def test() :
+    print "test x_piezo()"
+    xp = x_piezo()
+    xp.verbose = True
+    xp.jumpToPos(xp.nm2out(0))
+    xp.wander()
+    xp.wander()
+    xp.wander()
+    xp.pCurVals()
+    curPos = xp.curPos()
+    curPosPhys = xp.out_to_phys(curPos)
+    curPos2 = xp.phys_to_out(curPosPhys)
+    if curPos != curPos2 :
+        raise x_piezoError, "Conversion %d -> %g -> %d not round trip" % (curPos, curPosPhys, curPos2)
+    print "success"
+
+if __name__ == "__main__" :
+    test()
diff --git a/z_piezo.py b/z_piezo.py
new file mode 100644 (file)
index 0000000..ceb4d90
--- /dev/null
@@ -0,0 +1,397 @@
+"""The z_piezo module controls the one-dimensional displacement of a
+piezoelectric actuator, while allowing simultaneous data aquisition on
+two channels.  The module is called z_piezo, because for force spectroscopy,
+the axis that needs simultaneous aquisition is the Z axis, although this
+module could be used to provide control in any "fast" axis.
+
+There are a few globals controlling the behavior of the entire module.
+
+  USE_ABCD_DEFLECTION (default False)
+    selects between using a preprocessed vertical deflection signal
+    and using the 'raw' 4 photodiode output segments.
+  TEXT_VERBOSE  (default False)
+    print text messages to stderr displaying actions taken
+  PYLAB_INTERACTIVE_VERBOSE (default False)
+    display pylab graphs of any ramps taken
+  BASE_FIG_NUM (default 50)
+    the figure number for z_piezo pylab figures will be this number + {0,1,..}
+  LOG_RAMPS (default False)
+    save the input and output of any ramps to the LOG_DIR directory
+  LOG_DIR = '/home/wking/rsrch/data/z_piezo'
+    where the ramp logs are saved to
+  DEFAULT_ZP_CHAN (default 0)
+    output channel controlling the z piezo
+  DEFAULT_ZP_MON_CHAN (default 1)
+    input channel monitoring the z piezo signal
+  DEFAULT_DEF_CHAN (default 0, or in ABCD mode (2,3,4,5))
+    input channel(s) monitoring deflection (or some other feedback signal)
+"""
+
+USE_ABCD_DEFLECTION = False
+TEXT_VERBOSE = False
+PYLAB_INTERACTIVE_VERBOSE = False
+BASE_FIG_NUM = 50
+LOG_RAMPS = False
+LOG_DIR = '/home/wking/rsrch/data/z_piezo'
+
+DEFAULT_ZP_CHAN = 0
+DEFAULT_ZP_MON_CHAN = 1
+DEFAULT_DEF_CHAN = 0
+if USE_ABCD_DEFLECTION :
+    DEFAULT_DEF_CHAN = (2,3,4,5)
+    # The 4 photodiode segments are
+    #  A B
+    #  C D
+    # looking in along the laser
+
+import pycomedi.single_aio as single_aio
+import pycomedi.simult_aio as simult_aio
+
+from numpy import array, arange, ones, zeros, linspace, uint16, float, sin, pi
+from scipy.stats import linregress
+import data_logger
+
+if PYLAB_INTERACTIVE_VERBOSE :
+    from pylab import figure, plot, title, legend, hold, subplot
+    import time # for timestamping lines on plots
+
+class error (Exception) :
+    pass
+class outOfRange (error) :
+    pass
+
+# frontend:
+
+class z_piezo :
+    def __init__(self, zp_chan=DEFAULT_ZP_CHAN,
+                 zp_mon_chan=DEFAULT_ZP_MON_CHAN,
+                 def_chan=DEFAULT_DEF_CHAN) :
+        self.sensitivity = 7.41 # nm_surface/Volt_piezo
+        self.gain = 20 # Vpiezo / Voutput
+        self.nm2Vout = self.sensitivity * self.gain # nm_surface / V_output
+        self.chan_info = _chan_info(zp_chan, zp_mon_chan, def_chan)
+
+        self.curIn = array([0]*self.chan_info.numIn, dtype=uint16)
+        self.curOut = array([0]*self.chan_info.numOut, dtype=uint16)
+        self._jump = _z_piezo_jump(self.chan_info)
+        self._ramp = _z_piezo_ramp(self.chan_info)
+        self.zpMax = int(self.pos_nm2out(1000)) # limits at 1 micron
+        self.zpMin = int(self.pos_nm2out(-1000)) # and -1 micron
+        # define default value changed callbacks
+        self.setExternalZPiezo = None
+        self.setExternalZPiezoMon = None
+        self.setExternalDeflection = None
+    def __del__(self) :
+        pass
+    def curPos(self) :
+        return int(self.curOut[self.chan_info.zp_ind]) # cast as int so arithmetic is signed
+    def curPosMon(self) :
+        return int(self.curIn[self.chan_info.zp_mon_ind])
+    def curDef(self) :
+        if USE_ABCD_DEFLECTION :
+            A = int(self.curIn[self.chan_info.def_ind[0]])
+            B = int(self.curIn[self.chan_info.def_ind[1]])
+            C = int(self.curIn[self.chan_info.def_ind[2]])
+            D = int(self.curIn[self.chan_info.def_ind[3]])
+            df = float((A+B)-(C+D))/(A+B+C+D)
+            dfout = int(df * 2**15) + 2**15
+            if TEXT_VERBOSE :
+                print "Current deflection %d (%d, %d, %d, %d)" \
+                    % (dfout, A, B, C, D)
+            return dfout
+        else : # reading the deflection voltage directly
+            return int(self.curIn[self.chan_info.def_ind])
+    def pos_out2V(self, out) :
+        return self._jump.out_to_phys(out)
+    def pos_V2out(self, V) :
+        # Vout is Volts output by the DAQ card
+        # and Vpiezo is Volts applied across the piezo
+        return self._jump.phys_to_out(V)
+    def pos_nm2out(self, nm) :
+        # Vout is Volts output by the DAQ card
+        # and Vpiezo is Volts applied across the piezo
+        return self.pos_V2out(nm / self.nm2Vout)
+    def pos_out2nm(self, out) :
+        return self.pos_out2V(out) * self.nm2Vout
+    def posMon_nm2out(self, nm) :
+        return self._jump.phys_to_out(nm / self.nm2Vout)
+    def posMon_out2nm(self, out) :
+        return self._jump.out_to_phys(out) * self.nm2Vout
+    def def_in2V(self, input) :
+        return self._jump.out_to_phys(input)
+    def def_V2in(self, physical) :
+        print physical, type(physical) # temporary debugging printout
+        return self._jump.phys_to_out(physical)
+    def curVals(self) :
+        return {'Z piezo output': self.curPos(),
+                'Deflection input':self.curDef(),
+                'Z piezo input':self.curPosMon()}
+    def pCurVals(self, curVals=None) :
+        if curVals == None :
+            curVals = self.curVals()
+        print "Z piezo output  : %6d" % curVals["Z piezo output"]
+        print "Z piezo input   : %6d" % curVals["Z piezo input"]
+        print "Deflection input: %6d" % curVals["Deflection input"]
+    def updateInputs(self) :
+        self.curIn = self._jump.updateInputs()
+        if self.setExternalZPiezoMon != None :
+            self.setExternalZPiezoMon(self.posMon_out2nm(self.curPosMon()))
+        if self.setExternalDeflection != None :
+            self.setExternalDeflection(self.def_in2V(self.curDef()))
+    def jumpToPos(self, pos) :
+        self._check_range(pos)
+        self.curOut[self.chan_info.zp_ind] = pos
+        self.curIn = self._jump.jumpToPos(self.curOut)
+        if self.setExternalZPiezo != None :
+            self.setExternalZPiezo(self.pos_out2nm(self.curPos()))
+        if TEXT_VERBOSE :
+            print "Jumped to"
+            self.pCurVals()
+        return self.curVals()
+    def ramp(self, out, freq) :
+        for pos in out :
+            self._check_range(pos)
+        out = self._ramp.ramp(out, freq)
+        self.curOut[self.chan_info.zp_ind] = out["Z piezo output"][-1]
+        self.curIn[self.chan_info.zp_mon_ind] = out["Z piezo input"][-1]
+        if USE_ABCD_DEFLECTION :
+            for i in range(4) : # i is the photodiode element (0->A, 1->B, ...)
+                self.curIn[i] = out["Deflection segment"][i][-1]
+        else :
+            self.curIn[self.chan_info.def_ind] = out["Deflection input"][-1]
+
+        if self.setExternalZPiezo != None :
+            self.setExternalZPiezo(self.pos_out2nm(self.curPos()))
+        if self.setExternalZPiezoMon != None :
+            self.setExternalZPiezoMon(self.posMon_out2nm(self.curPosMon()))
+        if self.setExternalDeflection != None :
+            self.setExternalDeflection(self.def_in2V(self.curDef()))
+        if LOG_RAMPS :
+            log = data_logger.data_log(LOG_DIR, noclobber_logsubdir=False,
+                                       log_name="ramp")
+            if USE_ABCD_DEFLECTION :
+                v = out.pop("Deflection segment")
+                log.write_dict_of_arrays(out)
+                out["Deflection segment"] = v
+            else :
+                log.write_dict_of_arrays(out)
+        if TEXT_VERBOSE :
+            print "Ramped from (%g, %g) to (%g, %g) in %d points" \
+                % (out["Z piezo output"][0], out["Deflection input"][0],
+                   out["Z piezo output"][-1], out["Deflection input"][-1],
+                   len(out["Z piezo output"]))
+        if PYLAB_INTERACTIVE_VERBOSE :
+            figure(BASE_FIG_NUM)
+            subplot(211)
+            title('ramp (def_in vs zp_out)')
+            timestamp = time.strftime('%H%M%S')
+            plot(out["Z piezo output"], out["Deflection input"], '.', label=timestamp)
+            subplot(212)
+            title('zp_in vs zp_out')
+            plot(out["Z piezo output"], out["Z piezo input"], '.', label=timestamp)
+            legend(loc='best')
+        return out
+    def _check_range(self, pos) :
+        if pos > self.zpMax :
+            raise outOfRange, "Z piezo pos = %d > %d = max" % (pos, self.zpMax)
+        if pos < self.zpMin :
+            raise outOfRange, "Z piezo pos = %d < %d = min" % (pos, self.zpMin)
+
+# backends:
+class _z_piezo_jump :
+    def __init__(self, chan_info) :
+        self.chan_info = chan_info
+        self.reserved = False
+        if USE_COMEDI :
+            self.AI = comedi_single_aio.ai_obj(chan=self.chan_info.inChan)
+            self.AI.close()
+            self.AO = comedi_single_aio.ao_obj(chan=self.chan_info.outChan)
+            self.AO.close()
+        else :
+            self.AI = sngAI.sngAIobj(self.chan_info.inChan, -10, 10)
+            self.AO = sngAO.sngAOobj(self.chan_info.outChan, -10, 10)
+    def reserve(self) :
+        if USE_COMEDI :
+            self.AI.open()
+            self.AO.open()
+        else :
+            self.AI.reserve()
+            self.AO.reserve()
+        self.reserved = True
+    def release(self) :
+        if USE_COMEDI :
+            self.AI.close()
+            self.AO.close()
+        else :
+            self.AI.unreserve()
+            self.AO.unreserve()
+        self.reserved = False
+    def updateInputs(self) :
+        if self.reserved == False :
+            self.AI.open()
+            #self.reserve()
+            self.reserved = False
+        if USE_COMEDI :
+            In = self.AI.read()
+        else :
+            (er, In) = self.AI.read()
+        if self.reserved == False :
+            self.AI.close()
+            #self.release()
+        return In
+    def jumpToPos(self, out) :
+        if self.reserved == False :
+            self.reserve()
+            self.reserved = False # set to true inside reserve(), which we don't want
+        self.AO.write(out)
+        if self.reserved == False :
+            self.release()
+        return self.updateInputs()
+    def out_to_phys(self, output) :
+        if USE_COMEDI :
+            return self.AO.comedi_to_phys(0, output)
+    def phys_to_out(self, physical) :
+        if USE_COMEDI :
+            return self.AO.phys_to_comedi(0, physical)
+
+class _z_piezo_ramp :
+    def __init__(self, chan_info) :
+        self.chan_info = chan_info
+        self.verbose = True
+        self.attempts=0
+        self.failures=0
+    def ramp(self, outArray, freq) :
+        # allocate and run the task
+        npoints = int(len(outArray)/self.chan_info.numOut)
+        in_data = array([0]*npoints*self.chan_info.numIn, dtype=uint16)
+        if type(outArray) != type(in_data) or outArray.dtype != uint16 :
+            out_data = array([0]*npoints*self.chan_info.numOut, dtype=uint16)
+            for i in range(npoints) :
+                out_data[i] = outArray[i]
+        else :
+            out_data = outArray
+
+        correlated = False
+        ramp={}
+        while not correlated :
+            if USE_COMEDI :
+                AIO = comedi_simult_aio.aio_obj(in_chan=self.chan_info.inChan,
+                                                out_chan=self.chan_info.outChan)
+                if self.verbose :
+                    print "setup AIO"
+                AIO.setup(freq, out_data)
+                if self.verbose :
+                    print "arm AIO"
+                AIO.arm()
+                if self.verbose :
+                    print "read AIO"
+                AIO.start_read(in_data)
+                if self.verbose :
+                    print "close AIO"
+                AIO.close()
+                if self.verbose :
+                    print "finished AIO"
+            self.attempts += 1
+            ramp["Z piezo output"] = out_data[self.chan_info.zp_ind::self.chan_info.numOut]
+            ramp["Z piezo input"] = in_data[self.chan_info.zp_mon_ind::self.chan_info.numIn]
+            failed = False
+            gradient, intercept, r_value, p_value, std_err = linregress(ramp["Z piezo output"],
+                                                                        ramp["Z piezo input"])
+            rnge = ramp["Z piezo output"].max()-ramp["Z piezo output"].min()
+            if gradient < .8 and rnge > 100 :
+                print "ramp failed slope (%g, %g), try again" % (gradient, rnge)
+                failed = True
+            if not failed :
+                if USE_ABCD_DEFLECTION :
+                    ramp["Deflection segment"] = []
+                    for c_ind in self.chan_info.def_ind :
+                        ramp["Deflection segment"].append(in_data[c_ind::self.chan_info.numIn])
+                    A = ramp["Deflection segment"][0].astype(float)
+                    B = ramp["Deflection segment"][1].astype(float)
+                    C = ramp["Deflection segment"][2].astype(float)
+                    D = ramp["Deflection segment"][3].astype(float)
+                    ramp["Deflection input"] = (((A+B)-(C+D))/(A+B+C+D) * 2**15).astype(uint16)
+                    if USE_COMEDI :
+                        ramp["Deflection input"] += 2**15 # HACK, comedi uses unsigned ints...
+                else :
+                    ramp["Deflection input"] = in_data[self.chan_info.def_ind::self.chan_info.numIn]
+                correlated = True
+                if PYLAB_INTERACTIVE_VERBOSE :
+                    figure(BASE_FIG_NUM+1)
+                    timestamp = time.strftime('%H%M%S')
+                    plot(ramp["Z piezo output"], ramp["Z piezo input"], '.', label=timestamp)
+                    title('goods')
+                    legend(loc='best')
+            else :
+                self.failures += 1
+                if PYLAB_INTERACTIVE_VERBOSE :
+                    figure(BASE_FIG_NUM+2)
+                    timestamp = time.strftime('%H%M%S')
+                    plot(ramp["Z piezo output"], ramp["Z piezo input"], '.', label=timestamp)
+                    legend(loc='best')
+                    title('bads')
+
+        #AIO = simAIO.simAIOobj(npoints,
+        #                       freq,
+        #                       self.chan_info.inChan,
+        #                       self.chan_info.outChan)
+        #err = AIO.configure()
+        #err = AIO.commit()
+        #err = AIO.write(outArray)
+        #err = AIO.run()
+        #(err, inArray) = AIO.read()
+        #if err != 0 :
+        #    raise error, "Error in simAIO"
+        #if len(inArray)*self.chan_info.numOut != len(outArray)*self.chan_info.numIn :
+        #    raise error, "Input array of wrong length"
+        #
+        #ramp={}
+        #ramp["Z piezo output"] = outArray[self.chan_info.zp_ind::self.chan_info.numOut]
+        #ramp["Z piezo input"] = inArray[self.chan_info.zp_mon_ind::self.chan_info.numIn]
+        #ramp["Deflection input"] = inArray[self.chan_info.def_ind::self.chan_info.numIn]
+        
+        # task automatically freed as AIO goes out of scope
+        return ramp
+
+class _chan_info :
+    def __init__(self, zp_chan, zp_mon_chan, def_chan) :
+        self.inChan = []
+        if USE_ABCD_DEFLECTION :
+            self.def_ind = [] # index in inChan array
+            i=0
+            for chan in def_chan :
+                self.inChan.append(chan)
+                self.def_ind.append(i)
+                i += 1
+        else :
+            self.inChan.append(def_chan)
+            self.def_ind = 0 # index in inChan array
+        self.inChan.append(zp_mon_chan)
+        self.zp_mon_ind = zp_mon_chan
+        self.numIn = len(self.inChan)
+        
+        self.outChan = [zp_chan]
+        self.zp_ind = 0 # index in outChan array
+        self.numOut = len(self.outChan)
+        
+
+def test() :
+    print "test z_piezo()"
+    zp = z_piezo()
+    zp.verbose = True
+    zp.jumpToPos(zp.pos_nm2out(0))
+    if zp.verbose :
+        zp.pCurVals()
+    curPos = zp.curPos()
+    curPosNm = zp.pos_out2nm(curPos)
+    curPos2 = zp.pos_nm2out(curPosNm)
+    if curPos != curPos2 :
+        raise error, "Conversion %d -> %g -> %d not round trip" % (curPos, curPosNm, curPos2)
+    ramp_positions = linspace(zp.curPos(), zp.pos_nm2out(100), 20)
+    out = zp.ramp(ramp_positions, 10000)
+    ramp_positions = linspace(zp.curPos(), zp.pos_nm2out(0), 20)
+    out2 = zp.ramp(ramp_positions, 10000)
+
+if __name__ == "__main__" :
+    test()
diff --git a/z_piezo_utils.py b/z_piezo_utils.py
new file mode 100644 (file)
index 0000000..b5a2a54
--- /dev/null
@@ -0,0 +1,325 @@
+TEXT_VERBOSE = False
+PYLAB_INTERACTIVE_VERBOSE = False
+BASE_FIG_NUM = 60
+LOG_DATA = True
+LOG_DIR = '/home/wking/rsrch/data/z_piezo_utils'
+
+NULL_STEPS_BEFORE_SINGLE_PT_SWEEP = False
+SLEEP_DURING_SURF_POS_AQUISITION = False # doesn't help
+
+from numpy import array, arange, ones, zeros, linspace, uint16, float, sin, cos, pi
+from scipy.stats import linregress
+import curses_check_for_keypress
+import fit
+import data_logger
+import time
+
+if PYLAB_INTERACTIVE_VERBOSE :
+    from pylab import figure, plot, title, legend, hold, subplot
+        
+class error (Exception) :
+    pass
+class poorFit (error) :
+    pass
+class tooClose (error) :
+    pass
+
+def moveToPosOrDef(zpiezo, pos, defl, step, return_data=False) :
+    if return_data or LOG_DATA or PYLAB_INTERACTIVE_VERBOSE :
+        aquire_data = True
+    else :
+        aquire_data = False
+    zpiezo._check_range(pos)
+    if step == 0 :
+        raise error, "Must have non-zero step size"
+    elif step < 0 and pos > zpiezo.curPos() :
+        step = -step
+    elif step > 0 and pos < zpiezo.curPos() :
+        step = -step
+        
+    zpiezo._jump.reserve()
+    zpiezo.updateInputs()
+    if TEXT_VERBOSE :
+        print "current z_piezo position %6d, current deflection %6d <? %6d" % (zpiezo.curPos(), zpiezo.curDef(), defl)
+    if aquire_data :
+        def_array=[zpiezo.curDef()]
+        pos_array=[zpiezo.curPos()]
+    if NULL_STEPS_BEFORE_SINGLE_PT_SWEEP == True :
+        # take a few null steps to let the input deflection stabalize
+        for i in range(100) :
+            zpiezo.jumpToPos(zpiezo.curPos())
+            if aquire_data :
+                def_array.append(zpiezo.curDef())
+                pos_array.append(zpiezo.curPos())            
+    # step in until we hit our target position or exceed our target deflection
+    while zpiezo.curPos() != pos and zpiezo.curDef() < defl :
+        distTo = pos - zpiezo.curPos()
+        if abs(distTo) < abs(step) :
+            jumpTo = pos
+        else :
+            jumpTo = zpiezo.curPos() + step
+        zpiezo.jumpToPos(jumpTo)
+        if TEXT_VERBOSE :
+            print "current z_piezo position %6d, current deflection %6d >=< %6d" % (zpiezo.curPos(), zpiezo.curDef(), defl)
+        if aquire_data :
+            def_array.append(zpiezo.curDef())
+            pos_array.append(zpiezo.curPos())
+    zpiezo._jump.release()
+    
+    if zpiezo.setExternalZPiezo != None :
+        zpiezo.setExternalZPiezo(zpiezo.pos_out2nm(zpiezo.curPos()))
+    if PYLAB_INTERACTIVE_VERBOSE :
+        figure(BASE_FIG_NUM)
+        timestamp = time.strftime('%H%M%S')
+        plot(pos_array, def_array, '.', label=timestamp)
+        title('step approach')
+        legend(loc='best')
+    if LOG_DATA :
+        log = data_logger.data_log(LOG_DIR, noclobber_logsubdir=False,
+                                   log_name="moveToPosOrDef")
+        data = {"Z piezo output":array(pos_array),
+                "Deflection input":array(def_array)}
+        log.write_dict_of_arrays(data)
+    if return_data :
+        data = {"Z piezo output":array(pos_array),
+                "Deflection input":array(def_array)}
+        return (zpiezo.curVals(), data)
+    else :
+        return zpiezo.curVals()
+
+
+def wiggleForInterferenceMin(zpiezo, wig_amp=20000, wig_freq=100) :
+    "Output sin to measure interference"
+    if not PYLAB_INTERACTIVE_VERBOSE :
+        from pylab import figure, plot, title, hold, ion, ioff, draw, axes
+    npoints = 1000
+    scanfreq = wig_freq*npoints
+    out = zeros((npoints,), dtype=uint16)
+    for i in range(npoints) :
+        out[i] = int(sin(4*pi*i/npoints)*wig_amp+32000)
+        # 4 for 2 periods
+    figure(BASE_FIG_NUM+1)
+    a = axes()
+    if LOG_DATA :
+        log = data_logger.data_log(LOG_DIR, noclobber_logsubdir=False,
+                                   log_name="wiggle")
+    title('wiggle for interference')
+    hold(False)
+    p = plot(out, out, 'b.-')
+    c = curses_check_for_keypress.check_for_keypress()
+    while c.input() == None :
+        data = zpiezo.ramp(out, scanfreq)
+        p[0].set_ydata(data["Deflection input"])
+        dmin = data["Deflection input"].min()
+        dmax = data["Deflection input"].max()
+        a.set_ylim([dmin,dmax])
+        if LOG_DATA :
+            log.write_dict_of_arrays(data)
+        draw()
+    hold(True)
+
+def getSurfPosData(zpiezo, maxDefl, textVerbose=False, plotVerbose=False) :
+    "Measure the distance to the surface"
+    if plotVerbose and not PYLAB_INTERACTIVE_VERBOSE :
+        from pylab import figure, plot, title, legend, hold, subplot
+    origPos = zpiezo.curPos()
+    # fully retract the piezo
+    if TEXT_VERBOSE or textVerbose :
+        print "\tRetract the piezo"
+    out = linspace(zpiezo.curPos(), zpiezo.zpMin, 2000)
+    ret1 = zpiezo.ramp(out, 10e3)
+    del(out) # start with a clean output vector, seems to help reduce bounce?
+    zpiezo.updateInputs() # opening output seems to cause bounce?
+    # locate high deflection position
+    if TEXT_VERBOSE or textVerbose :
+        print "\tApproach until there is dangerous deflection (> %d)" % maxDefl
+    if SLEEP_DURING_SURF_POS_AQUISITION == True :
+        time.sleep(.2) # sleeping briefly seems to reduce bounce?
+    cp, mtpod = moveToPosOrDef(zpiezo, zpiezo.zpMax, maxDefl,
+                               (zpiezo.zpMax-zpiezo.zpMin)/2000,
+                               return_data=True)
+    highContactPos = zpiezo.curPos()
+    if highContactPos <= origPos :   # started out too close, we need some room to move in...
+        if TEXT_VERBOSE or textVerbose :
+            print "\tToo close, return to the original position: ", origPos
+        out = linspace(zpiezo.curPos(), origPos, 2000)
+        ret3=zpiezo.ramp(out, 10e3)
+        del(out)
+        zpiezo.updateInputs()
+        if PYLAB_INTERACTIVE_VERBOSE or plotVerbose :
+            figure(BASE_FIG_NUM+2)
+            def plot_dict(d, label) :
+                print d.keys()
+                plot(d["Z piezo output"], d["Deflection input"], '.',label=label)
+            plot_dict(ret1, 'First retract')
+            plot_dict(mtpod, 'Single step in')
+            plot_dict(ret3, 'Return to start')
+            legend(loc='best')
+            title('Too close')
+        raise tooClose, "Way too close!\n(highContactPos %d <= origPos %d)" % (highContactPos, origPos) 
+    zpiezo.updateInputs()
+    # fully retract the piezo again
+    if TEXT_VERBOSE or textVerbose :
+        print "\tRetract the piezo again"
+    if SLEEP_DURING_SURF_POS_AQUISITION == True :
+        time.sleep(.2)
+    out = linspace(highContactPos, zpiezo.zpMin, 2000)
+    ret2 = zpiezo.ramp(out, 10e3)
+    del(out)
+    zpiezo.updateInputs()
+    if zpiezo.curPos() != zpiezo.zpMin :
+        raise Exception, "bad code"
+    # scan to the high contact position
+    if TEXT_VERBOSE or textVerbose :
+        print "\tRamp in to the deflected position: ", highContactPos
+    if SLEEP_DURING_SURF_POS_AQUISITION == True :
+        time.sleep(.2)
+    out = linspace(zpiezo.curPos(), highContactPos, 4000)
+    data = zpiezo.ramp(out, 10e3)
+    del(out)
+    zpiezo.updateInputs()
+    if LOG_DATA :
+        log = data_logger.data_log(LOG_DIR, noclobber_logsubdir=False,
+                                   log_name="getSurfPos")
+        log.write_dict_of_arrays(data)
+    if SLEEP_DURING_SURF_POS_AQUISITION == True :
+        time.sleep(.2)
+    # return to the original position
+    if TEXT_VERBOSE or textVerbose :
+        print "\tReturn to the original position: ", origPos
+    out = linspace(highContactPos, origPos, 2000)
+    ret3=zpiezo.ramp(out, 10e3)
+    del(out)
+    zpiezo.updateInputs()
+    return {"ret1":ret1, "mtpod":mtpod, "ret2":ret2,
+            "approach":data, "ret3":ret3}
+
+def analyzeSurfPosData(ddict, textVerbose=False, plotVerbose=False, retAllParams=False) :
+    # ususes ddict["approach"] for analysis
+    # the others are just along to be plotted
+
+    data = ddict["approach"]
+    # analyze data, using bilinear model
+    # y = p0 + p1 x                for x <= p2
+    #   = p0 + p1 p2 + p3 (x-p2)   for x >= p2
+    if TEXT_VERBOSE or textVerbose :
+        print "\tAnalyze the data"
+    dump_before_index = 0 # 25 # HACK!!
+    # Generate a reasonable guess...
+    startPos = int(data["Z piezo output"].min())
+    finalPos = int(data["Z piezo output"].max())
+    startDef = int(data["Deflection input"].min())
+    finalDef = int(data["Deflection input"].max())
+    # startDef and startPos are probably for 2 different points
+    
+    left_offset   = startDef
+    left_slope    = 0
+    kink_position = (finalPos+startPos)/2.0
+    right_slope   = 2.0*(finalDef-startDef)/(finalPos-startPos)
+    pstart = [left_offset, left_slope, kink_position, right_slope]
+
+    offset_scale = (finalPos - startPos)/100
+    left_slope_scale = right_slope/10
+    kink_scale = (finalPos-startPos)/100
+    right_slope_scale = right_slope
+    scale = [offset_scale, left_slope_scale, kink_scale, right_slope_scale]
+    if TEXT_VERBOSE or textVerbose :
+        print "start Def %d, final def %d, start pos %d, final pos %d" % (startDef, finalDef, startPos, finalPos)
+        print "Guessed params ", pstart
+    params = fit.fit(data=[data["Z piezo output"][dump_before_index:], data["Deflection input"][dump_before_index:]],
+                     model='bilinear', start=pstart, scale=scale)
+    if TEXT_VERBOSE or textVerbose :
+        print "Best fit parameters: ", params
+
+    if PYLAB_INTERACTIVE_VERBOSE or plotVerbose :
+        if plotVerbose and not PYLAB_INTERACTIVE_VERBOSE :
+            from pylab import figure, plot, title, legend, hold, subplot
+        figure(BASE_FIG_NUM+3)
+        def plot_dict(d, label) :
+            plot(d["Z piezo output"], d["Deflection input"], '.',label=label)
+        try :
+            plot_dict(ddict['ret1'], 'First retract')
+            plot_dict(ddict['mtpod'], 'Single step in')
+            plot_dict(ddict['ret2'], 'Second retract')
+        except KeyError :
+            pass # if we don't have the surrounding data, don't worry
+        plot_dict(data, 'Main approach')
+        try :
+            plot_dict(ddict['ret3'], 'Return to start')
+        except KeyError :
+            pass # if we don't have the surrounding data, don't worry
+        def fit_fn(x, params) :
+            if x <= params[2] :
+                return params[0] + params[1]*x
+            else :
+                return params[0] + params[1]*params[2] + params[3]*(x-params[2])
+        plot([startPos, params[2], finalPos],
+             [fit_fn(startPos, params), fit_fn(params[2], params), fit_fn(finalPos, params)],
+             '-',label='fit')
+        title('surf_pos %5g %5g %5g %5g' % (params[0], params[1], params[2], params[3]))
+        legend(loc='best')
+
+    # check that the fit is reasonable
+    # params[1] is slope in non-contact region
+    # params[3] is slope in contact region
+    if abs(params[1]*10) > abs(params[3]) :
+        raise poorFit, "Slope in non-contact region, or no slope in contact"
+    if params[2] < startPos+0.02*(finalPos-startPos) : # params[2] is kink position
+        raise poorFit, "No kink (kink %g less than %g, need more space to left)" % \
+            (params[2], startPos+0.02*(finalPos-startPos))
+    if params[2] > finalPos-0.02*(finalPos-startPos) :
+        raise poorFit, "No kink (kink %g more than %g, need more space to right)" % \
+            (params[2], finalPos-0.02*(finalPos-startPos))
+    if params[3] < 0.5 * right_slope : # params[3] is slope in contact region
+        raise poorFit, "Too far"
+    if TEXT_VERBOSE or textVerbose :
+        print "\tSurface position: ", params[2]
+    if retAllParams :
+        return params
+    return params[2]
+
+def getSurfPos(zpiezo, maxDefl, textVerbose=False, plotVerbose=False) :
+    ddict = getSurfPosData(zpiezo, maxDefl, textVerbose, plotVerbose)
+    return analyzeSurfPosData(ddict, textVerbose, plotVerbose)
+
+def test_smoothness(zp, plotVerbose=True) :
+    posA = 20000
+    posB = 50000
+    setpoint = zp.def_V2in(3)
+    steps = 200
+    outfreq = 1e5
+    outarray = linspace(posB, posA, 1000)
+    indata=[]
+    outdata=[]
+    curVals = zp.jumpToPos(posA)
+    zp.pCurVals(curVals)
+    time.sleep(1) # let jitters die down
+    for i in range(10) :
+        print "ramp %d to %d" % (zp.curPos(), posB)
+        curVals, data = moveToPosOrDef(zp, posB, setpoint, step=steps,
+                                       return_data = True)
+        indata.append(data)
+        out = zp.ramp(outarray, outfreq)
+        outdata.append(out)
+    if plotVerbose :
+        from pylab import figure, plot, title, legend, hold, subplot        
+    if PYLAB_INTERACTIVE_VERBOSE or plotVerbose :
+        figure(BASE_FIG_NUM+4)
+        for i in range(10) :
+            plot(indata[i]['Z piezo output'],
+                 indata[i]['Deflection input'], '+--', label='in')
+            plot(outdata[i]['Z piezo output'],
+                 outdata[i]['Deflection input'], '.-', label='out')
+        title('test smoothness (step in, ramp out)')
+        #legend(loc='best')
+    
+def test() :
+    import z_piezo
+    zp = z_piezo.z_piezo()
+    curVals = zp.moveToPosOrDef(zp.pos_nm2out(600), defl=zp.curDef()+6000, step=(zp.pos_nm2out(10)-zp.pos_nm2out(0)))
+    if TEXT_VERBOSE :
+        zp.pCurVals(curVals)
+    pos = zp.getSurfPos(maxDefl=zp.curDef()+6000)
+    if TEXT_VERBOSE :
+        print "Surface at %g nm", pos
+    print "success"