calibcant/analyze.py

   1 # calibcant - tools for thermally calibrating AFM cantilevers
   2 #
   3 # Copyright (C) 2008-2011 W. Trevor King <wking@drexel.edu>
   4 #
   5 # This file is part of calibcant.
   6 #
   7 # calibcant is free software: you can redistribute it and/or
   8 # modify it under the terms of the GNU Lesser General Public
   9 # License as published by the Free Software Foundation, either
  10 # version 3 of the License, or (at your option) any later version.
  11 #
  12 # calibcant is distributed in the hope that it will be useful,
  13 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  14 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  15 # GNU Lesser General Public License for more details.
  16 #
  17 # You should have received a copy of the GNU Lesser General Public
  18 # License along with calibcant.  If not, see
  19 # <http://www.gnu.org/licenses/>.
  20
  21 """
  22 Separate the more general calib_analyze() from the other calib_*()
  23 functions in calibcant.  Also provide a command line interface
  24 for analyzing data acquired through other workflows.
  25
  26 The relevent physical quantities are :
  27  Vzp_out  Output z-piezo voltage (what we generate)
  28  Vzp      Applied z-piezo voltage (after external ZPGAIN)
  29  Zp       The z-piezo position
  30  Zcant    The cantilever vertical deflection
  31  Vphoto   The photodiode vertical deflection voltage (what we measure)
  32  Fcant    The force on the cantilever
  33  T        The temperature of the cantilever and surrounding solution
  34           (another thing we measure)
  35  k_b      Boltzmann's constant
  36
  37 Which are related by the parameters :
  38  zpGain           Vzp_out / Vzp
  39  zpSensitivity    Zp / Vzp
  40  photoSensitivity Vphoto / Zcant
  41  k_cant           Fcant / Zcant
  42 """
  43
  44 import numpy
  45
  46 import data_logger
  47 from splittable_kwargs import splittableKwargsFunction, \
  48     make_splittable_kwargs_function
  49
  50 from . import common
  51 from . import config
  52 from . import T_analyze
  53
  54
  55 kb = 1.3806504e-23 # Boltzmann's constant in J/K
  56
  57 #@splittableKwargsFunction((calib_plot, 'bumps', 'Ts', 'vibs'))
  58 # Some of the child functions aren't yet defined, so postpone
  59 # make-splittable until later in the module.
  60 def calib_analyze(bumps, Ts, vibs, **kwargs) :
  61     """
  62     Analyze data from get_calibration_data()
  63     return (k, k_s, !!!a2_r, T_r, one_o_Vp2_r)
  64     Inputs (all are arrays of recorded data) :
  65      bumps measured (V_photodiode / nm_tip) proportionality constant
  66      Ts    measured temperature (K)
  67      vibs  measured V_photodiode variance in free solution (V**2)
  68     Outputs :
  69      k    cantilever spring constant (in N/m, or equivalently nN/nm)
  70      k_s  standard deviation in our estimate of k
  71      !!!a2_r relative error in a**2
  72      !!!T_r  relative error in T
  73      !!!one_o_Vp2_r relative error in 1/Vphotodiode_variance
  74     Notes :
  75      We're assuming vib is mostly from thermal cantilever vibrations
  76     (and then only from vibrations in the single vertical degree of freedom),
  77     and not from other noise sources.
  78      The various relative errors are returned to help you gauge the
  79     largest source of random error in your measurement of k.
  80     If one of them is small, don't bother repeating that measurment too often.
  81     If one is large, try repeating that measurement more.
  82     Remember that you need enough samples to have a valid error estimate in
  83     the first place, and that none of this addresses any systematic errors.
  84     """
  85     calib_plot_kwargs, = calib_analyze._splitargs(calib_analyze, kwargs)
  86     photoSensitivity2 = bumps**2
  87     one_o_Vphoto2 = 1/vibs
  88
  89     photoSensitivity2_m = photoSensitivity2.mean()
  90     T_m = Ts.mean()
  91     one_o_Vphoto2_m = one_o_Vphoto2.mean()
  92     # Vphoto / photoSensitivity = x
  93     # k = kb T / <x**2> = kb T photoSensitiviy**2 * (1e9nm/m)**2 /
  94     #                                                       <Vphoto_std**2>
  95     #
  96     # units,  photoSensitivity =  Vphoto/(Zcant in nm),
  97     # so Vphoto/photoSensitivity = Zcant in nm
  98     # so k = J/K * K / nm^2 * (1e9nm/m)**2 = N/m
  99     k  = kb * T_m * photoSensitivity2_m * one_o_Vphoto2_m * 1e18
 100
 101     # propogation of errors !!!
 102     # first, get standard deviations
 103     photoSensitivity2_s = photoSensitivity2.std()
 104     T_s = Ts.std()
 105     one_o_Vphoto2_s = one_o_Vphoto2.std()
 106     # !!!!now, get relative errors
 107     photoSensitivity2_r = photoSensitivity2_s / photoSensitivity2_m
 108     T_r = T_s / T_m
 109     one_o_Vphoto2_r = one_o_Vphoto2_s / one_o_Vphoto2_m
 110
 111     k_s = k*(photoSensitivity2_r**2 + T_r**2 + one_o_Vphoto2_r**2)**0.5
 112
 113     calib_plot(bumps, Ts, vibs, **calib_plot_kwargs)
 114
 115     return (k, k_s,
 116             photoSensitivity2_m, photoSensitivity2_s,
 117             T_m, T_s, one_o_Vphoto2_m, one_o_Vphoto2_s)
 118
 119 @splittableKwargsFunction()
 120 def string_errors(k_m, k_s,
 121                   photoSensitivity2_m, photoSensitivity2_s,
 122                   T_m, T_s,
 123                   one_o_Vphoto2_m, one_o_Vphoto2_s) :
 124     k_r = k_s / k_m
 125     photoSensitivity2_r = photoSensitivity2_s / photoSensitivity2_m
 126     T_r = T_s / T_m
 127     one_o_Vphoto2_r = one_o_Vphoto2_s / one_o_Vphoto2_m
 128     string  = "Variable (units)              : mean +/- std. dev. (relative error)\n"
 129     string += "Cantilever k (N/m)            : %g +/- %g (%g)\n" \
 130               % (k_m, k_s, k_r)
 131     string += "photoSensitivity**2 (V/nm)**2 : %g +/- %g (%g)\n" \
 132               % (photoSensitivity2_m, photoSensitivity2_s, photoSensitivity2_r)
 133     string += "T (K)                         : %g +/- %g (%g)\n" \
 134               % (T_m, T_s, T_r)
 135     string += "1/Vphoto**2 (1/V)**2          : %g +/- %g (%g)" \
 136               % (one_o_Vphoto2_m, one_o_Vphoto2_s, one_o_Vphoto2_r)
 137     return string
 138
 139 @splittableKwargsFunction()
 140 def calib_save(bumps, Ts, vibs, log_dir=None) :
 141     """
 142     Save a dictonary with the bump, T, and vib data.
 143     """
 144     if log_dir != None :
 145         data = {'bump':bumps, 'T':Ts, 'vib':vibs}
 146         log = data_logger.data_log(log_dir, noclobber_logsubdir=False,
 147                                    log_name="calib")
 148         log.write_dict_of_arrays(data)
 149
 150 def calib_load(datafile) :
 151     "Load the dictionary data, using data_logger.date_load()"
 152     dl = data_logger.data_load()
 153     data = dl.read_dict_of_arrays(path)
 154     return (data['bump'], data['T'], data['vib'])
 155
 156 def calib_save_analysis(k, k_s,
 157                         photoSensitivity2_m, photoSensitivity2_s,
 158                         T_m, T_s, one_o_Vphoto2_m, one_o_Vphoto2_s,
 159                         log_dir=None) :
 160     if log_dir != None :
 161         log = data_logger.data_log(log_dir, noclobber_logsubdir=False,
 162                                    log_name="calib_analysis_text")
 163         log.write_binary(string_errors(k, k_s,
 164                                        photoSensitivity2_m, photoSensitivity2_s,
 165                                        T_m, T_s,
 166                                        one_o_Vphoto2_m, one_o_Vphoto2_s)
 167                          +'\n')
 168
 169 @splittableKwargsFunction()
 170 def calib_plot(bumps, Ts, vibs, plotVerbose=False) :
 171     if plotVerbose or config.PYLAB_VERBOSE :
 172         common._import_pylab()
 173         common._pylab.figure(config.BASE_FIGNUM+4)
 174         common._pylab.subplot(311)
 175         common._pylab.plot(bumps, 'g.-')
 176         common._pylab.title('Photodiode sensitivity (V/nm)')
 177         common._pylab.subplot(312)
 178         common._pylab.plot(Ts, 'r.-')
 179         common._pylab.title('Temperature (K)')
 180         common._pylab.subplot(313)
 181         common._pylab.plot(vibs, 'b.-')
 182         common._pylab.title('Thermal deflection variance (Volts**2)')
 183         common._flush_plot()
 184
 185 make_splittable_kwargs_function(calib_analyze,
 186                                 (calib_plot, 'bumps', 'Ts', 'vibs'))
 187
 188 @splittableKwargsFunction((calib_analyze, 'bumps', 'Ts', 'vibs'))
 189 def calib_load_analyze_tweaks(bump_tweaks, vib_tweaks, T_tweaks=None) :
 190     raise NotImplementedError
 191     a = read_tweaked_bumps(bump_tweaks)
 192     vib = V_photo_variance_from_file(vib_tweaks)
 193     if T_tweaks == None:
 194         pass
 195     return analyze_calibration_data(a, T, vib, log_dir=log_dir)
 196
 197 # commandline interface functions
 198 import scipy.io, sys
 199
 200 def array_from_string(string):
 201     ret = []
 202     for num in string.split(',') :
 203         ret.append(float(num))
 204     assert len(ret) > 0
 205     return numpy.array(ret)
 206
 207 def read_data(ifile):
 208     "ifile can be a filename string or open (seekable) file object"
 209     unlabeled_data=scipy.io.read_array(file)
 210     return unlabeled_data
 211
 212 def get_array(string, filename, name) :
 213     "get an array from supplied command line options"
 214     if string != None :
 215         array = array_from_string(string)
 216     elif filename != None :
 217         array = read_data(filename)
 218     else :
 219         raise Exception, "no %s information given" % (name)
 220     return array
 221
 222 if __name__ == '__main__' :
 223     # command line interface
 224     from optparse import OptionParser
 225
 226     usage_string = ('%prog <bumps> <temps> <vibs>\n'
 227                     '2008, W. Trevor King.\n'
 228                     '\n'
 229                     'Takes arrays of Vphotodiode sensitivity (V/nm), Temperature (K), \n'
 230                     'and Vibration variance (V**2) as comma seperated lists.\n'
 231                     'Returns the cantilever spring constant (pN/nm).\n'
 232                     'for example:\n'
 233                     ' $ %prog -b 0.02,0.03,0.025 -t 298.2,300.1 -v 6e-9,5.5e-9\n'
 234                     )
 235     parser = OptionParser(usage=usage_string, version='%prog '+common.VERSION)
 236     parser.add_option('-b', '--bump-string', dest='bump_string',
 237                       help='comma seperated photodiode sensitivities (V/nm)',
 238                       type='string', metavar='BUMPS')
 239     parser.add_option('-t', '--temp-string', dest='temp_string',
 240                       help='comma seperated temperatures (K)',
 241                       type='string', metavar='TEMPS')
 242     parser.add_option('-v', '--vib-string', dest='vib_string',
 243                       help='comma seperated vibration variances (V**2)',
 244                       type='string', metavar='VIBS')
 245     parser.add_option('-B', '--bump-file', dest='bump_file',
 246                       help='comma seperated photodiode sensitivities (V/nm)',
 247                       type='string', metavar='BUMPFILE')
 248     parser.add_option('-T', '--temp-file', dest='temp_file',
 249                       help='comma seperated temperatures (K)',
 250                       type='string', metavar='TEMPFILE')
 251     parser.add_option('-V', '--vib-file', dest='vib_file',
 252                       help='comma seperated vibration variances (V**2)',
 253                       type='string', metavar='VIBFILE')
 254     parser.add_option('-C', '--celsius', dest='celsius',
 255                       help='Use Celsius input temperatures instead of Kelvin (defaul %default)\n',
 256                       action='store_true', default=False)
 257     parser.add_option('-o', '--output-file', dest='ofilename',
 258                       help='write output to FILE (default stdout)',
 259                       type='string', metavar='FILE')
 260     parser.add_option('-p', '--plot-inputs', dest='plot',
 261                       help='plot the input arrays to check their distribution',
 262                       action='store_true', default=False)
 263     parser.add_option('', '--verbose', dest='verbose', action='store_true',
 264                       help='print lots of debugging information',
 265                       default=False)
 266
 267     options,args = parser.parse_args()
 268     parser.destroy()
 269
 270     config.TEXT_VERBOSE = options.verbose
 271     config.PYLAB_INTERACTIVE = False
 272     config.PYLAB_VERBOSE = options.plot
 273
 274     bumps = get_array(options.bump_string, options.bump_file, "bump")
 275     temps = get_array(options.temp_string, options.temp_file, "temp")
 276     vibs = get_array(options.vib_string, options.vib_file, "vib")
 277
 278     #if options.plot :
 279     #    calib_plot(bumps, temps, vibs)
 280
 281     if options.celsius :
 282         for i in range(len(temps)) :
 283             temps[i] = T_analyze.C_to_K(temps[i])
 284
 285     km,ks,ps2m,ps2s,Tm,Ts,ooVp2m,ooVp2s = \
 286         calib_analyze(bumps, temps, vibs, plotVerbose=options.plot)
 287
 288     if options.verbose :
 289         print >> sys.stderr, \
 290             string_errors(km,ks,ps2m,ps2s,Tm,Ts,ooVp2m,ooVp2s)
 291
 292     if options.ofilename != None :
 293         print >> file(options.ofilename, 'w'), km
 294     else :
 295         print km