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