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(string_errors(k, k_s,
 167                                        photoSensitivity2_m, photoSensitivity2_s,
 168                                        T_m, T_s,
 169                                        one_o_Vphoto2_m, one_o_Vphoto2_s)
 170                          +'\n')
 171
 172 @splittableKwargsFunction()
 173 def calib_plot(bumps, Ts, vibs, plotVerbose=False) :
 174     if plotVerbose or config.PYLAB_VERBOSE :
 175         common._import_pylab()
 176         common._pylab.figure(config.BASE_FIGNUM+4)
 177         common._pylab.subplot(311)
 178         common._pylab.plot(bumps, 'g.-')
 179         common._pylab.title('Photodiode sensitivity (V/nm)')
 180         common._pylab.subplot(312)
 181         common._pylab.plot(Ts, 'r.-')
 182         common._pylab.title('Temperature (K)')
 183         common._pylab.subplot(313)
 184         common._pylab.plot(vibs, 'b.-')
 185         common._pylab.title('Thermal deflection variance (Volts**2)')
 186         common._flush_plot()
 187
 188 make_splittable_kwargs_function(calib_analyze,
 189                                 (calib_plot, 'bumps', 'Ts', 'vibs'))
 190
 191 @splittableKwargsFunction((calib_analyze, 'bumps', 'Ts', 'vibs'))
 192 def calib_load_analyze_tweaks(bump_tweaks, vib_tweaks, T_tweaks=None) :
 193     raise NotImplementedError
 194     a = read_tweaked_bumps(bump_tweaks)
 195     vib = V_photo_variance_from_file(vib_tweaks)
 196     if T_tweaks == None:
 197         pass
 198     return analyze_calibration_data(a, T, vib, log_dir=log_dir)
 199
 200 # commandline interface functions
 201 import scipy.io, sys
 202
 203 def array_from_string(string):
 204     ret = []
 205     for num in string.split(',') :
 206         ret.append(float(num))
 207     assert len(ret) > 0
 208     return numpy.array(ret)
 209
 210 def read_data(ifile):
 211     "ifile can be a filename string or open (seekable) file object"
 212     unlabeled_data=scipy.io.read_array(file)
 213     return unlabeled_data
 214
 215 def get_array(string, filename, name) :
 216     "get an array from supplied command line options"
 217     if string != None :
 218         array = array_from_string(string)
 219     elif filename != None :
 220         array = read_data(filename)
 221     else :
 222         raise Exception, "no %s information given" % (name)
 223     return array
 224
 225 if __name__ == '__main__' :
 226     # command line interface
 227     from optparse import OptionParser
 228
 229     usage_string = ('%prog <bumps> <temps> <vibs>\n'
 230                     '2008, W. Trevor King.\n'
 231                     '\n'
 232                     'Takes arrays of Vphotodiode sensitivity (V/nm), Temperature (K), \n'
 233                     'and Vibration variance (V**2) as comma seperated lists.\n'
 234                     'Returns the cantilever spring constant (pN/nm).\n'
 235                     'for example:\n'
 236                     ' $ %prog -b 0.02,0.03,0.025 -t 298.2,300.1 -v 6e-9,5.5e-9\n'
 237                     )
 238     parser = OptionParser(usage=usage_string, version='%prog '+common.VERSION)
 239     parser.add_option('-b', '--bump-string', dest='bump_string',
 240                       help='comma seperated photodiode sensitivities (V/nm)',
 241                       type='string', metavar='BUMPS')
 242     parser.add_option('-t', '--temp-string', dest='temp_string',
 243                       help='comma seperated temperatures (K)',
 244                       type='string', metavar='TEMPS')
 245     parser.add_option('-v', '--vib-string', dest='vib_string',
 246                       help='comma seperated vibration variances (V**2)',
 247                       type='string', metavar='VIBS')
 248     parser.add_option('-B', '--bump-file', dest='bump_file',
 249                       help='comma seperated photodiode sensitivities (V/nm)',
 250                       type='string', metavar='BUMPFILE')
 251     parser.add_option('-T', '--temp-file', dest='temp_file',
 252                       help='comma seperated temperatures (K)',
 253                       type='string', metavar='TEMPFILE')
 254     parser.add_option('-V', '--vib-file', dest='vib_file',
 255                       help='comma seperated vibration variances (V**2)',
 256                       type='string', metavar='VIBFILE')
 257     parser.add_option('-C', '--celsius', dest='celsius',
 258                       help='Use Celsius input temperatures instead of Kelvin (defaul %default)\n',
 259                       action='store_true', default=False)
 260     parser.add_option('-o', '--output-file', dest='ofilename',
 261                       help='write output to FILE (default stdout)',
 262                       type='string', metavar='FILE')
 263     parser.add_option('-p', '--plot-inputs', dest='plot',
 264                       help='plot the input arrays to check their distribution',
 265                       action='store_true', default=False)
 266     parser.add_option('', '--verbose', dest='verbose', action='store_true',
 267                       help='print lots of debugging information',
 268                       default=False)
 269
 270     options,args = parser.parse_args()
 271     parser.destroy()
 272
 273     config.TEXT_VERBOSE = options.verbose
 274     config.PYLAB_INTERACTIVE = False
 275     config.PYLAB_VERBOSE = options.plot
 276
 277     bumps = get_array(options.bump_string, options.bump_file, "bump")
 278     temps = get_array(options.temp_string, options.temp_file, "temp")
 279     vibs = get_array(options.vib_string, options.vib_file, "vib")
 280
 281     #if options.plot :
 282     #    calib_plot(bumps, temps, vibs)
 283
 284     if options.celsius :
 285         for i in range(len(temps)) :
 286             temps[i] = T_analyze.C_to_K(temps[i])
 287
 288     km,ks,ps2m,ps2s,Tm,Ts,ooVp2m,ooVp2s = \
 289         calib_analyze(bumps, temps, vibs, plotVerbose=options.plot)
 290
 291     if options.verbose :
 292         print >> sys.stderr, \
 293             string_errors(km,ks,ps2m,ps2s,Tm,Ts,ooVp2m,ooVp2s)
 294
 295     if options.ofilename != None :
 296         print >> file(options.ofilename, 'w'), km
 297     else :
 298         print km