calibcant/analyze.py

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