calibcant/analyze.py

   1 #!/usr/bin/python
   2 #
   3 # calibcant - tools for thermally calibrating AFM cantilevers
   4 #
   5 # Copyright (C) 2008-2010 W. Trevor King <wking@drexel.edu>
   6 #
   7 # This file is part of CalibCant.
   8 #
   9 # CalibCant is free software: you can redistribute it and/or
  10 # modify it under the terms of the GNU Lesser General Public
  11 # License as published by the Free Software Foundation, either
  12 # version 3 of the License, or (at your option) any later version.
  13 #
  14 # CalibCant is distributed in the hope that it will be useful,
  15 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  16 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  17 # GNU Lesser General Public License for more details.
  18 #
  19 # You should have received a copy of the GNU Lesser General Public
  20 # License along with CalibCant.  If not, see
  21 # <http://www.gnu.org/licenses/>.
  22
  23 """
  24 Separate the more general calib_analyze() from the other calib_*()
  25 functions in calibcant.  Also provide a command line interface
  26 for analyzing data acquired through other workflows.
  27
  28 The relevent physical quantities are :
  29  Vzp_out  Output z-piezo voltage (what we generate)
  30  Vzp      Applied z-piezo voltage (after external ZPGAIN)
  31  Zp       The z-piezo position
  32  Zcant    The cantilever vertical deflection
  33  Vphoto   The photodiode vertical deflection voltage (what we measure)
  34  Fcant    The force on the cantilever
  35  T        The temperature of the cantilever and surrounding solution
  36           (another thing we measure)
  37  k_b      Boltzmann's constant
  38
  39 Which are related by the parameters :
  40  zpGain           Vzp_out / Vzp
  41  zpSensitivity    Zp / Vzp
  42  photoSensitivity Vphoto / Zcant
  43  k_cant           Fcant / Zcant
  44 """
  45
  46 import numpy
  47
  48 import data_logger
  49 from splittable_kwargs import splittableKwargsFunction, \
  50     make_splittable_kwargs_function
  51
  52 import .common
  53 import .config
  54 import .T_analyze
  55
  56
  57 kb = 1.3806504e-23 # Boltzmann's constant in J/K
  58
  59 #@splittableKwargsFunction((calib_plot, 'bumps', 'Ts', 'vibs'))
  60 # Some of the child functions aren't yet defined, so postpone
  61 # make-splittable until later in the module.
  62 def calib_analyze(bumps, Ts, vibs, **kwargs) :
  63     """
  64     Analyze data from get_calibration_data()
  65     return (k, k_s, !!!a2_r, T_r, one_o_Vp2_r)
  66     Inputs (all are arrays of recorded data) :
  67      bumps measured (V_photodiode / nm_tip) proportionality constant
  68      Ts    measured temperature (K)
  69      vibs  measured V_photodiode variance in free solution (V**2)
  70     Outputs :
  71      k    cantilever spring constant (in N/m, or equivalently nN/nm)
  72      k_s  standard deviation in our estimate of k
  73      !!!a2_r relative error in a**2
  74      !!!T_r  relative error in T
  75      !!!one_o_Vp2_r relative error in 1/Vphotodiode_variance
  76     Notes :
  77      We're assuming vib is mostly from thermal cantilever vibrations
  78     (and then only from vibrations in the single vertical degree of freedom),
  79     and not from other noise sources.
  80      The various relative errors are returned to help you gauge the
  81     largest source of random error in your measurement of k.
  82     If one of them is small, don't bother repeating that measurment too often.
  83     If one is large, try repeating that measurement more.
  84     Remember that you need enough samples to have a valid error estimate in
  85     the first place, and that none of this addresses any systematic errors.
  86     """
  87     calib_plot_kwargs, = calib_analyze._splitargs(calib_analyze, kwargs)
  88     photoSensitivity2 = bumps**2
  89     one_o_Vphoto2 = 1/vibs
  90
  91     photoSensitivity2_m = photoSensitivity2.mean()
  92     T_m = Ts.mean()
  93     one_o_Vphoto2_m = one_o_Vphoto2.mean()
  94     # Vphoto / photoSensitivity = x
  95     # k = kb T / <x**2> = kb T photoSensitiviy**2 * (1e9nm/m)**2 /
  96     #                                                       <Vphoto_std**2>
  97     #
  98     # units,  photoSensitivity =  Vphoto/(Zcant in nm),
  99     # so Vphoto/photoSensitivity = Zcant in nm
 100     # so k = J/K * K / nm^2 * (1e9nm/m)**2 = N/m
 101     k  = kb * T_m * photoSensitivity2_m * one_o_Vphoto2_m * 1e18
 102
 103     # propogation of errors !!!
 104     # first, get standard deviations
 105     photoSensitivity2_s = photoSensitivity2.std()
 106     T_s = Ts.std()
 107     one_o_Vphoto2_s = one_o_Vphoto2.std()
 108     # !!!!now, get relative errors
 109     photoSensitivity2_r = photoSensitivity2_s / photoSensitivity2_m
 110     T_r = T_s / T_m
 111     one_o_Vphoto2_r = one_o_Vphoto2_s / one_o_Vphoto2_m
 112
 113     k_s = k*(photoSensitivity2_r**2 + T_r**2 + one_o_Vphoto2_r**2)**0.5
 114
 115     calib_plot(bumps, Ts, vibs, **calib_plot_kwargs)
 116
 117     return (k, k_s,
 118             photoSensitivity2_m, photoSensitivity2_s,
 119             T_m, T_s, one_o_Vphoto2_m, one_o_Vphoto2_s)
 120
 121 @splittableKwargsFunction()
 122 def string_errors(k_m, k_s,
 123                   photoSensitivity2_m, photoSensitivity2_s,
 124                   T_m, T_s,
 125                   one_o_Vphoto2_m, one_o_Vphoto2_s) :
 126     k_r = k_s / k_m
 127     photoSensitivity2_r = photoSensitivity2_s / photoSensitivity2_m
 128     T_r = T_s / T_m
 129     one_o_Vphoto2_r = one_o_Vphoto2_s / one_o_Vphoto2_m
 130     string  = "Variable (units)              : mean +/- std. dev. (relative error)\n"
 131     string += "Cantilever k (N/m)            : %g +/- %g (%g)\n" \
 132               % (k_m, k_s, k_r)
 133     string += "photoSensitivity**2 (V/nm)**2 : %g +/- %g (%g)\n" \
 134               % (photoSensitivity2_m, photoSensitivity2_s, photoSensitivity2_r)
 135     string += "T (K)                         : %g +/- %g (%g)\n" \
 136               % (T_m, T_s, T_r)
 137     string += "1/Vphoto**2 (1/V)**2          : %g +/- %g (%g)" \
 138               % (one_o_Vphoto2_m, one_o_Vphoto2_s, one_o_Vphoto2_r)
 139     return string
 140
 141 @splittableKwargsFunction()
 142 def calib_save(bumps, Ts, vibs, log_dir=None) :
 143     """
 144     Save a dictonary with the bump, T, and vib data.
 145     """
 146     if log_dir != None :
 147         data = {'bump':bumps, 'T':Ts, 'vib':vibs}
 148         log = data_logger.data_log(log_dir, noclobber_logsubdir=False,
 149                                    log_name="calib")
 150         log.write_dict_of_arrays(data)
 151
 152 def calib_load(datafile) :
 153     "Load the dictionary data, using data_logger.date_load()"
 154     dl = data_logger.data_load()
 155     data = dl.read_dict_of_arrays(path)
 156     return (data['bump'], data['T'], data['vib'])
 157
 158 def calib_save_analysis(k, k_s,
 159                         photoSensitivity2_m, photoSensitivity2_s,
 160                         T_m, T_s, one_o_Vphoto2_m, one_o_Vphoto2_s,
 161                         log_dir=None) :
 162     if log_dir != None :
 163         log = data_logger.data_log(log_dir, noclobber_logsubdir=False,
 164                                    log_name="calib_analysis_text")
 165         log.write_binary(string_errors(k, k_s,
 166                                        photoSensitivity2_m, photoSensitivity2_s,
 167                                        T_m, T_s,
 168                                        one_o_Vphoto2_m, one_o_Vphoto2_s)
 169                          +'\n')
 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     config.TEXT_VERBOSE = options.verbose
 273     config.PYLAB_INTERACTIVE = False
 274     config.PYLAB_VERBOSE = options.plot
 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 = \
 288         calib_analyze(bumps, temps, vibs, plotVerbose=options.plot)
 289
 290     if options.verbose :
 291         print >> sys.stderr, \
 292             string_errors(km,ks,ps2m,ps2s,Tm,Ts,ooVp2m,ooVp2s)
 293
 294     if options.ofilename != None :
 295         print >> file(options.ofilename, 'w'), km
 296     else :
 297         print km