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