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