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