- # units, photoSensitivity = Vphoto/(Zcant in nm),
- # so Vphoto/photoSensitivity = Zcant in nm
- # so k = J/K * K / nm^2 * (1e9nm/m)**2 = N/m
- k = kb * T_m * photoSensitivity2_m * one_o_Vphoto2_m * 1e18
-
- # propogation of errors !!!
- # first, get standard deviations
- photoSensitivity2_s = photoSensitivity2.std()
- T_s = Ts.std()
- one_o_Vphoto2_s = one_o_Vphoto2.std()
- # !!!!now, get relative errors
- photoSensitivity2_r = photoSensitivity2_s / photoSensitivity2_m
- T_r = T_s / T_m
- one_o_Vphoto2_r = one_o_Vphoto2_s / one_o_Vphoto2_m
-
- k_s = k*(photoSensitivity2_r**2 + T_r**2 + one_o_Vphoto2_r**2)**0.5
-
- return (k, k_s,
- photoSensitivity2_m, photoSensitivity2_s,
- T_m, T_s, one_o_Vphoto2_m, one_o_Vphoto2_s)
-
-def string_errors(k_m, k_s,
- photoSensitivity2_m, photoSensitivity2_s,
- T_m, T_s,
- one_o_Vphoto2_m, one_o_Vphoto2_s) :
- k_r = k_s / k_m
- photoSensitivity2_r = photoSensitivity2_s / photoSensitivity2_m
- T_r = T_s / T_m
- one_o_Vphoto2_r = one_o_Vphoto2_s / one_o_Vphoto2_m
- string = "Variable (units) : mean +/- std. dev. (relative error)\n"
- string += "Cantilever k (N/m) : %g +/- %g (%g)\n" \
- % (k_m, k_s, k_r)
- string += "photoSensitivity**2 (V/nm)**2 : %g +/- %g (%g)\n" \
- % (photoSensitivity2_m, photoSensitivity2_s, photoSensitivity2_r)
- string += "T (K) : %g +/- %g (%g)\n" \
- % (T_m, T_s, T_r)
- string += "1/Vphoto**2 (1/V)**2 : %g +/- %g (%g)" \
- % (one_o_Vphoto2_m, one_o_Vphoto2_s, one_o_Vphoto2_r)
- return string
-
-def calib_plot(bumps, Ts, vibs) :
- from pylab import figure, subplot, plot, title, show
- figure()
- subplot(311)
- plot(bumps, 'g.-')
- title('Photodiode sensitivity (V/nm)')
- subplot(312)
- plot(Ts, 'r.-')
- title('Temperature (K)')
- subplot(313)
- plot(vibs, 'b.-')
- title('Thermal deflection variance (Volts**2)')
- show()
-
-# commandline interface functions
-import scipy.io, sys
-
-def array_from_string(string):
- ret = []
- for num in string.split(',') :
- ret.append(float(num))
- assert len(ret) > 0
- return numpy.array(ret)
-
-def read_data(ifile):
- "ifile can be a filename string or open (seekable) file object"
- unlabeled_data=scipy.io.read_array(file)
- return unlabeled_data
-
-def get_array(string, filename, name) :
- "get an array from supplied command line options"
- if string != None :
- array = array_from_string(string)
- elif filename != None :
- array = read_data(filename)
- else :
- raise Exception, "no %s information given" % (name)
- return array
-
-if __name__ == '__main__' :
- # command line interface
- from optparse import OptionParser
-
- usage_string = ('%prog <bumps> <temps> <vibs>\n'
- '2008, W. Trevor King.\n'
- '\n'
- 'Takes arrays of Vphotodiode sensitivity (V/nm), Temperature (K), \n'
- 'and Vibration variance (V**2) as comma seperated lists.\n'
- 'Returns the cantilever spring constant (pN/nm).\n'
- 'for example:\n'
- ' $ %prog -b 0.02,0.03,0.025 -t 298.2,300.1 -v 6e-9,5.5e-9\n'
- )
- parser = OptionParser(usage=usage_string, version='%prog '+common.VERSION)
- parser.add_option('-b', '--bump-string', dest='bump_string',
- help='comma seperated photodiode sensitivities (V/nm)',
- type='string', metavar='BUMPS')
- parser.add_option('-t', '--temp-string', dest='temp_string',
- help='comma seperated temperatures (K)',
- type='string', metavar='TEMPS')
- parser.add_option('-v', '--vib-string', dest='vib_string',
- help='comma seperated vibration variances (V**2)',
- type='string', metavar='VIBS')
- parser.add_option('-B', '--bump-file', dest='bump_file',
- help='comma seperated photodiode sensitivities (V/nm)',
- type='string', metavar='BUMPFILE')
- parser.add_option('-T', '--temp-file', dest='temp_file',
- help='comma seperated temperatures (K)',
- type='string', metavar='TEMPFILE')
- parser.add_option('-V', '--vib-file', dest='vib_file',
- help='comma seperated vibration variances (V**2)',
- type='string', metavar='VIBFILE')
- parser.add_option('-C', '--celsius', dest='celsius',
- help='Use Celsius input temperatures instead of Kelvin (defaul %default)\n',
- action='store_true', default=False)
- parser.add_option('-o', '--output-file', dest='ofilename',
- help='write output to FILE (default stdout)',
- type='string', metavar='FILE')
- parser.add_option('-p', '--plot-inputs', dest='plot',
- help='plot the input arrays to check their distribution',
- action='store_true', default=False)
- parser.add_option('', '--verbose', dest='verbose', action='store_true',
- help='print lots of debugging information',
- default=False)
-
- options,args = parser.parse_args()
- parser.destroy()
-
- bumps = get_array(options.bump_string, options.bump_file, "bump")
- temps = get_array(options.temp_string, options.temp_file, "temp")
- vibs = get_array(options.vib_string, options.vib_file, "vib")
-
- if options.plot :
- calib_plot(bumps, temps, vibs)
-
- if options.celsius :
- for i in range(len(temps)) :
- temps[i] = T_analyze.C_to_K(temps[i])
-
- km,ks,ps2m,ps2s,Tm,Ts,ooVp2m,ooVp2s = calib_analyze(bumps, temps, vibs)
-
- if options.verbose :
- print >> sys.stderr, \
- string_errors(km,ks,ps2m,ps2s,Tm,Ts,ooVp2m,ooVp2s)
-
- if options.ofilename != None :
- print >> file(options.ofilename, 'w'), km
- else :
- print km
+ # units, photo_sensitivity = Vphoto/(Zcant in m),
+ # so Vphoto/photo_sensitivity = Zcant in m
+ # so k = J/K * K / m^2 = J / m^2 = N/m
+ k = _kB * T_m * ps_m**2 / v2_m
+
+ # propogation of errors
+ # dk/dT = k/T
+ dk_T = k/T_m * T_s
+ # dk/dps = 2k/ps
+ dk_ps = 2*k/ps_m * ps_s
+ # dk/dv2 = -k/v2
+ dk_v = -k/v2_m * v2_s
+
+ k_s = _numpy.sqrt(dk_T**2 + dk_ps**2 + dk_v**2)
+
+ _LOG.info('variable (units) : '
+ 'mean +/- std. dev. (relative error)')
+ _LOG.info('cantilever k (N/m) : %g +/- %g (%g)' % (k, k_s, k_s/k))
+ _LOG.info('photo sensitivity (V/m) : %g +/- %g (%g)'
+ % (ps_m, ps_s, ps_s/ps_m))
+ _LOG.info('T (K) : %g +/- %g (%g)'
+ % (T_m, T_s, T_s/T_m))
+ _LOG.info('vibration variance (V^2) : %g +/- %g (%g)'
+ % (v2_m, v2_s, v2_s/v2_m))
+
+ if _package_config['matplotlib']:
+ plot(bumps, temperatures, vibrations)
+
+ return (k, k_s)
+
+
+def plot(bumps, temperatures, vibrations):
+ if not _matplotlib:
+ raise _matplotlib_import_error
+ figure = _matplotlib_pyplot.figure()
+
+ bump_axes = figure.add_subplot(3, 1, 1)
+ T_axes = figure.add_subplot(3, 1, 2)
+ vib_axes = figure.add_subplot(3, 1, 3)
+
+ timestamp = _time.strftime('%H%M%S')
+ bump_axes.set_title('cantilever calibration %s' % timestamp)
+
+ bump_axes.autoscale(tight=True)
+ bump_axes.plot(bumps, 'g.-')
+ bump_axes.set_ylabel('photodiode sensitivity (V/m)')
+ T_axes.autoscale(tight=True)
+ T_axes.plot(temperatures, 'r.-')
+ T_axes.set_ylabel('temperature (K)')
+ vib_axes.autoscale(tight=True)
+ vib_axes.plot(vibrations, 'b.-')
+ vib_axes.set_ylabel('thermal deflection variance (V^2)')
+
+ if hasattr(figure, 'show'):
+ figure.show()
+ return figure
+_plot = plot # alternative name for use inside analyze_all()
+
+def save_results(filename=None, group='/', bump=None,
+ temperature=None, vibration=None, spring_constant=None,
+ spring_constant_deviation=None):
+ specs = [
+ _SaveSpec(item=bump, relpath='raw/photodiode-sensitivity',
+ array=True, units='V/m'),
+ _SaveSpec(item=temperature, relpath='raw/temperature',
+ array=True, units='K'),
+ _SaveSpec(item=vibration, relpath='raw/vibration',
+ array=True, units='V^2/Hz'),
+ _SaveSpec(item=spring_constant, relpath='processed/spring-constant',
+ units='N/m', deviation=spring_constant_deviation),
+ ]
+ _save(filename=filename, group=group, specs=specs)
+
+def analyze_all(config, data, raw_data, maximum_relative_error=1e-5,
+ filename=None, group=None, plot=False, dry_run=False):
+ "(Re)analyze (and possibly plot) all data from a `calib()` run."
+ if not data.get('bump', None):
+ data['bump'] = _numpy.zeros((config['num-bumps'],), dtype=float)
+ if not data.get('temperature', None):
+ data['temperature'] = _numpy.zeros(
+ (config['num-temperatures'],), dtype=float)
+ if not data.get('vibrations', None):
+ data['vibration'] = _numpy.zeros(
+ (config['num-vibrations'],), dtype=float)
+ if 'raw' not in data:
+ data['raw'] = {}
+ if 'bump' not in data['raw']:
+ data['raw']['bump'] = _numpy.zeros((config['num-bumps'],), dtype=float)
+ if 'temperature' not in data['raw']:
+ data['raw']['temperature'] = _numpy.zeros(
+ (config['num-temperatures'],), dtype=float)
+ if 'vibration' not in data['raw']:
+ data['raw']['vibration'] = _numpy.zeros(
+ (config['num-vibrations'],), dtype=float)
+ axis_config = config['afm']['piezo'].select_config(
+ setting_name='axes',
+ attribute_value=config['afm']['main-axis'],
+ get_attribute=_get_axis_name)
+ input_config = config['afm']['piezo'].select_config(
+ setting_name='inputs', attribute_value='deflection')
+ calibration_group = None
+ if not isinstance(group, _h5py.Group) and not dry_run:
+ f = _h5py.File(filename, mode='a')
+ group = _h5_create_group(f, group)
+ else:
+ f = None
+ try:
+ bumps_changed = len(data['raw']['bump']) != len(data['bump'])
+ for i,bump in enumerate(raw_data.get('bump', [])): # compare values
+ data['bump'][i],changed = check_bump(
+ index=i, bump=bump, config=config, z_axis_config=axis_config,
+ deflection_channel_config=input_config, plot=plot,
+ maximum_relative_error=maximum_relative_error)
+ if changed and not dry_run:
+ bumps_changed = True
+ bump_group = _h5_create_group(group, 'bump/{}'.format(i))
+ _bump_save(group=bump_group, processed=data['bump'][i])
+ temperatures_changed = len(data['raw']['temperature']) != len(
+ data['temperature'])
+ for i,temperature in enumerate(raw_data.get('temperature', [])):
+ data['temperature'][i],changed = check_temperature(
+ index=i, temperature=temperature, config=config,
+ maximum_relative_error=maximum_relative_error)
+ if changed and not dry_run:
+ temperatures_changed = True
+ temperature_group = _h5_create_group(
+ group, 'temperature/{}'.format(i))
+ _temperature_save(
+ group=temperature_group, processed=data['temperature'][i])
+ vibrations_changed = len(data['raw']['vibration']) != len(
+ data['vibration'])
+ for i,vibration in enumerate(raw_data.get('vibration', [])):
+ data['vibration'][i],changed = check_vibration(
+ index=i, vibration=vibration, config=config,
+ deflection_channel_config=input_config, plot=plot,
+ maximum_relative_error=maximum_relative_error)
+ if changed and not dry_run:
+ vibrations_changed = True
+ vibration_group = _h5_create_group(
+ group, 'vibration/{}'.format(i))
+ _vibration_save(
+ group=vibration_group, processed=data['vibration'][i])
+ if (bumps_changed or temperatures_changed or vibrations_changed
+ ) and not dry_run:
+ calibration_group = _h5_create_group(group, 'calibration')
+ if bumps_changed:
+ save_results(
+ group=calibration_group, bump=data['bump'])
+ if temperatures_changed:
+ save_results(
+ group=calibration_group, temperature=data['temperature'])
+ if vibrations_changed:
+ save_results(
+ group=calibration_group, vibration=data['vibration'])
+ if len(raw_data.get('bump', [])) != len(data['bump']):
+ raise ValueError(
+ 'not enough raw bump data: {} of {}'.format(
+ len(raw_data.get('bump', [])), len(data['bump'])))
+ if len(raw_data.get('temperature', [])) != len(data['temperature']):
+ raise ValueError(
+ 'not enough raw temperature data: {} of {}'.format(
+ len(raw_data.get('temperature', [])),
+ len(data['temperature'])))
+ if len(raw_data['vibration']) != len(data['vibration']):
+ raise ValueError(
+ 'not enough raw vibration data: {} of {}'.format(
+ len(raw_data.get('vibration', [])),
+ len(data['vibration'])))
+ k,k_s,changed = check_calibration(
+ k=data.get('processed', {}).get('spring_constant', None),
+ k_s=data.get('processed', {}).get(
+ 'spring_constant_deviation', None),
+ bumps=data['bump'],
+ temperatures=data['temperature'], vibrations=data['vibration'],
+ maximum_relative_error=maximum_relative_error)
+ if changed and not dry_run:
+ if calibration_group is None:
+ calibration_group = _h5_create_group(group, 'calibration')
+ save_results(
+ group=calibration_group,
+ spring_constant=k, spring_constant_deviation=k_s)
+ finally:
+ if f:
+ f.close()
+ if plot:
+ _plot(bumps=data['bump'],
+ temperatures=data['temperature'],
+ vibrations=data['vibration'])
+ return (k, k_s)
+
+def check_bump(index, bump, config=None, maximum_relative_error=0, **kwargs):
+ changed = False
+ try:
+ bump_config = bump['config']['bump']
+ except KeyError:
+ bump_config = config['bump']
+ sensitivity = _bump_analyze(
+ config=bump_config, data=bump['raw'], **kwargs)
+ if bump.get('processed', None) is None:
+ changed = True
+ _LOG.warn('new analysis for bump {}: {}'.format(index, sensitivity))
+ else:
+ rel_error = abs(sensitivity - bump['processed'])/bump['processed']
+ if rel_error > maximum_relative_error:
+ changed = True
+ _LOG.warn(("new analysis doesn't match for bump {}: {} -> {} "
+ "(difference: {}, relative error: {})").format(
+ index, bump['processed'], sensitivity,
+ sensitivity-bump['processed'], rel_error))
+ return (sensitivity, changed)
+
+def check_temperature(index, temperature, config=None,
+ maximum_relative_error=0, **kwargs):
+ changed = False
+ try:
+ temp_config = temperature['config']['temperature']
+ except KeyError:
+ temp_config = config['temperature']
+ temp = _temperature_analyze(
+ config=temp_config,
+ temperature=temperature['raw'], **kwargs)
+ if temperature.get('processed', None) is None:
+ changed = True
+ _LOG.warn('new analysis for temperature {}: {}'.format(index, temp))
+ else:
+ rel_error = abs(temp - temperature['processed']
+ )/temperature['processed']
+ if rel_error > maximum_relative_error:
+ changed = True
+ _LOG.warn(("new analysis doesn't match for temperature "
+ "{} -> {} (difference: {}, relative error: {})"
+ ).format(
+ index, temperature['processed'], temp,
+ temp-temperature['processed'], rel_error))
+ return (temp, changed)
+
+def check_vibration(index, vibration, config=None, maximum_relative_error=0,
+ **kwargs):
+ changed = False
+ try:
+ vib_config = vibration['config']['vibration']
+ except KeyError:
+ vib_config = config['vibration']
+ variance = _vibration_analyze(
+ config=vib_config, deflection=vibration['raw'], **kwargs)
+ if vibration.get('processed', None) is None:
+ changed = True
+ _LOG.warn('new analysis for vibration {}: {}'.format(
+ index, variance))
+ else:
+ rel_error = abs(variance-vibration['processed'])/vibration['processed']
+ if rel_error > maximum_relative_error:
+ changed = True
+ _LOG.warn(("new analysis doesn't match for vibration {}: {} != {} "
+ "(difference: {}, relative error: {})").format(
+ index, variance, vibration['processed'],
+ variance-vibration['processed'], rel_error))
+ return (variance, changed)
+
+def check_calibration(k, k_s, maximum_relative_error, **kwargs):
+ changed = False
+ new_k,new_k_s = analyze(**kwargs)
+ if k is None:
+ changed = True
+ _LOG.warn('new analysis for the spring constant: {}'.format(new_k))
+ else:
+ rel_error = abs(new_k-k)/k
+ if rel_error > maximum_relative_error:
+ changed = True
+ _LOG.warn(("new analysis doesn't match for the spring constant: "
+ "{} != {} (difference: {}, relative error: {})").format(
+ new_k, k, new_k-k, rel_error))
+ if k_s is None:
+ changed = True
+ _LOG.warn('new analysis for the spring constant deviation: {}'.format(
+ new_k_s))
+ else:
+ rel_error = abs(new_k_s-k_s)/k_s
+ if rel_error > maximum_relative_error:
+ changed = True
+ _LOG.warn(
+ ("new analysis doesn't match for the spring constant deviation"
+ ": {} != {} (difference: {}, relative error: {})").format(
+ new_k_s, k_s, new_k_s-k_s, rel_error))
+ return (new_k, new_k_s, changed)