def _load_cocktail_computed(self, an, new_list):
        if new_list:
            c = an.arar_constants
            ratios = []
            refs = {'40Ar/38Ar': nominal_value(c.atm4038),
                    '40Ar/36Ar': nominal_value(c.atm4036)}

            if self.recall_options:
                names = [r.tagname for r in self.recall_options.cocktail_options.ratios]
                ratios = [(name, name, refs.get(name, 1)) for name in names]

            cv = self._make_ratios(ratios)

            an.calculate_age()
            cv.append(ComputedValue(name='F', tag='uf',
                                    uvalue=an.uF))

            cv.append(ComputedValue(name='40Ar*', tag='rad40_percent',
                                    uvalue=an.rad40_percent))

            cv.append(ComputedValue(name='Age', tag='uage',
                                    uvalue=an.uage))

            self.computed_values = cv
            self._update_ratios()
        else:
            self._update_ratios()

def calculate_atmospheric(a38, a36, k38, ca38, ca36, decay_time,
                          production_ratios=None,
                          arar_constants=None):
    """
        McDougall and Harrison
        Roddick 1983
        Foland 1993

        iteratively calculate atm36
    """
    if production_ratios is None:
        production_ratios = {}

    if arar_constants is None:
        arar_constants = ArArConstants()

    pr = production_ratios
    m = pr.get('Cl3638', 0) * nominal_value(arar_constants.lambda_Cl36) * decay_time
    atm3836 = nominal_value(arar_constants.atm3836)

    # lcl36 = arar_constants.lambda_Cl36.nominal_value
    # atm3836 = arar_constants.atm3836.nominal_value

    # m = pr.get('Cl3638', 0) * lcl36 * decay_time
    # atm36 = ufloat(0, 0, tag='atm36')
    atm36 = 0
    for _ in range(5):
        ar38atm = atm3836 * atm36
        cl38 = a38 - ar38atm - k38 - ca38
        cl36 = cl38 * m
        atm36 = a36 - ca36 - cl36
    return atm36, cl36, cl38

def test_counter():
    if (sys.version_info < (2, 6, 0)):
        from nose.plugins.skip import SkipTest
        raise SkipTest    
    mcu = Arduino()
    reg = mcu.registers.proxy
    mcu.pins.reset()
    p = mcu.pin(5)
    p.write_mode(OUTPUT)
    p.pwm.write_value(128)
    print 'frequencies_available:', p.pwm.frequencies_available
    for fset in p.pwm.frequencies_available:
        p.pwm.frequency = fset
        assert abs(p.pwm.frequency - fset) <= 1
        print '---------------------------'
        print 'fset=', fset
        print '---------------------------'
        for ms in [10, 20, 50, 100, 200, 500, 1000]:
            for _ in range(1):
                t = ms / 1000.0
                with mcu.counter:
                    mcu.counter.run(t)
                    f = mcu.counter.frequency
                    t = mcu.counter.gate_time
                    err = f - fset

                    print 't=%s  f=%s ' % (t, f)
                    ok_(abs(nominal_value(err)) <= 0.1+std_dev(err), 
                        (abs(nominal_value(err)),std_dev(err)))

    def _plot_ratio(self, po, i):
        xs = [nominal_value(ai) for ai in self._unpack_attr(po.xtitle)]
        ys = [nominal_value(ai) for ai in self._unpack_attr(po.ytitle)]

        self.graph.new_series(x=array(xs), y=array(ys),
                              add_inspector=False,
                              marker=po.marker,
                              marker_size=po.marker_size)

    def _load_unknown_computed(self, an, new_list):
        attrs = (('Age', 'uage_w_j_err'),
                 # ('Age', 'age', None, None, 'age_err'),
                 ('w/o J', 'wo_j', '', 'uage', 'age_err_wo_j'),
                 ('K/Ca', 'kca'),
                 ('K/Cl', 'kcl'),
                 ('40Ar*', 'rad40_percent'),
                 ('F', 'uF'),
                 ('w/o Irrad', 'wo_irrad', '', 'uF', 'F_err_wo_irrad'))

        if new_list:
            def comp_factory(n, a, value=None, value_tag=None, error_tag=None):
                if value is None:
                    value = getattr(an, a)

                display_value = True
                if value_tag:
                    value = getattr(an, value_tag)
                    display_value = False

                if error_tag:
                    e = getattr(an, error_tag)
                else:
                    e = std_dev(value)

                return ComputedValue(name=n,
                                     tag=a,
                                     value=nominal_value(value) or 0,
                                     value_tag=value_tag or '',
                                     display_value=display_value,
                                     error=e or 0)

            cv = [comp_factory(*args)
                  for args in attrs]

            self.computed_values = cv
        else:
            age = an.uage
            nage, sage = nominal_value(age), std_dev(age)
            try:
                self.summary_str = u'Age={} {}{}({}%)'.format(floatfmt(nage), PLUSMINUS,
                                                              floatfmt(sage), format_percent_error(nage, sage))
            except:
                pass

            for ci in self.computed_values:
                attr = ci.tag
                if attr == 'wo_j':
                    ci.error = an.age_err_wo_j or 0
                    ci.value = nominal_value(getattr(an, ci.value_tag))
                elif attr == 'wo_irrad':
                    ci.error = an.F_err_wo_irrad or 0
                    ci.value = nominal_value(getattr(an, ci.value_tag))
                else:
                    v = getattr(an, attr)
                    if v is not None:
                        ci.value = nominal_value(v)
                        ci.error = std_dev(v)

    def _build_integrated_age_label(self, tga, n):
        txt = 'NaN'
        if not isnan(nominal_value(tga)):
            age, error = nominal_value(tga.nominal_value), std_dev(tga)

            error *= self.options.nsigma
            txt = self._build_label_text(age, error, n, sig_figs=self.options.integrated_sig_figs)

        return u'Total Integrated Age= {}'.format(txt)

 def calculate_decay_factors(self):
     arc = self.arar_constants
     # only calculate decayfactors once
     if not self.ar39decayfactor:
         dc37 = nominal_value(arc.lambda_Ar37)
         dc39 = nominal_value(arc.lambda_Ar39)
         a37df, a39df = calculate_arar_decay_factors(dc37, dc39, self.chron_segments)
         self.ar37decayfactor = a37df
         self.ar39decayfactor = a39df

    def get_baseline_corrected_value(self, include_baseline_error=None):
        if include_baseline_error is None:
            include_baseline_error = self.include_baseline_error

        b = self.baseline.uvalue
        if not include_baseline_error:
            b = nominal_value(b)
            nv = self.uvalue - b
            return ufloat(nominal_value(nv), std_dev(nv), tag=self.name)
        else:
            return self.uvalue - b

    def _value_string(self, t):
        if t == 'uF':
            a, e = self.F, self.F_err
        elif t == 'uage':
            a, e = nominal_value(self.uage), std_dev(self.uage)
        else:
            v = self.get_value(t)
            if isinstance(v, Isotope):
                v = v.get_intensity()
            a, e = nominal_value(v), std_dev(v)

        return a, e

    def _make_intermediate_summary(self, sh, ag, cols, label):
        row = self._current_row

        age_idx = next((i for i, c in enumerate(cols) if c.label == 'Age'), 0)
        cum_idx = next((i for i, c in enumerate(cols) if c.attr == 'cumulative_ar39'), 0)

        fmt = self._get_number_format('summary_age')
        kcafmt = self._get_number_format('summary_kca')

        fmt.set_bottom(1)
        kcafmt.set_bottom(1)

        fmt2 = self._workbook.add_format({'bottom': 1, 'bold': True})
        border = self._workbook.add_format({'bottom': 1})

        for i in range(age_idx + 1):
            sh.write_blank(row, i, '', fmt)

        startcol = 1
        sh.write(row, startcol, '{:02n}'.format(ag.aliquot), fmt2)
        sh.write_rich_string(row, startcol + 1, label, fmt2)
        cols[startcol + 1].calculate_width(label)

        age = ag.uage
        tn = ag.total_n
        if label == 'plateau':
            if not ag.plateau_steps:
                age = None
            else:
                txt = 'n={}/{} steps={}'.format(ag.nsteps, tn, ag.plateau_steps_str)
                sh.write(row, startcol + 2, txt, border)
                sh.write(row, cum_idx + 1, format_mswd(ag.get_plateau_mswd_tuple()), border)

        else:
            txt = 'n={}/{}'.format(ag.nanalyses, tn)
            sh.write(row, startcol + 2, txt, border)
            sh.write(row, cum_idx + 1, format_mswd(ag.get_mswd_tuple()), border)

        if age is not None:
            sh.write_number(row, age_idx, nominal_value(age), fmt)
            sh.write_number(row, age_idx + 1, std_dev(age), fmt)
        else:
            sh.write(row, age_idx, 'No plateau', border)

        sh.write_number(row, age_idx + 2, nominal_value(ag.kca), kcafmt)
        sh.write_number(row, age_idx + 3, std_dev(ag.kca), kcafmt)

        if label == 'plateau':
            sh.write_number(row, cum_idx, ag.plateau_total_ar39(), fmt)
        else:
            sh.write_number(row, cum_idx, ag.valid_total_ar39(), fmt)
        self._current_row += 1

    def __init__(self, ratios, name, *args, **kw):
        super(Result, self).__init__(*args, **kw)
        vs = array([nominal_value(ri) for ri in ratios])
        es = array([std_dev(ri) for ri in ratios])

        self.name = name
        m = ratios.mean()
        self.value = nominal_value(m)
        self.error = std_dev(m)
        wm, we = calculate_weighted_mean(vs, es)
        self.wm_value = wm
        self.wm_error = we
        self.mswd = calculate_mswd(vs, es, wm=wm)

    def _get_value(self, attr, n=3, **kw):
        v = ""
        item = self.item
        if hasattr(item, attr):
            v = getattr(self.item, attr)
            if v:
                v = floatfmt(nominal_value(v), n=n, **kw)
        elif hasattr(item, "isotopes"):
            if attr in item.isotopes:
                v = item.isotopes[attr].get_intensity()
                v = floatfmt(nominal_value(v), n=n, **kw)

        return v

    def _air_ratio(self):
        a4038 = self.isotope_group.get_ratio('Ar40/Ar38', non_ic_corr=True)
        a4036 = self.isotope_group.get_ratio('Ar40/Ar36', non_ic_corr=True)
        # e4038 = uformat_percent_error(a4038, include_percent_sign=True)
        # e4036 = uformat_percent_error(a4036, include_percent_sign=True)

        lines = [self._make_header('Ratios'),
                 'Ar40/Ar36= {} {}'.format(floatfmt(nominal_value(a4036)), errorfmt(nominal_value(a4036),
                                                                                    std_dev(a4036))),
                 'Ar40/Ar38= {} {}'.format(floatfmt(nominal_value(a4038)), errorfmt(nominal_value(a4038),
                                                                                    std_dev(a4038)))]

        return self._make_lines(lines)

	def __residenceTime(self):
		"""
			Calculate the residence time of a single step event. 
		"""
		# if previous errors were detected, the 
		# event is already rejected, don't process it 
		# any further
		if self.mdProcessingStatus != 'normal':
			return

		# set numpy warning handling. 
		# raise divide by zero errors so we 
		# can catch them
		np.seterr(divide='raise')

		ocmu=np.abs(uncertainties.nominal_value(self.mdOpenChCurrent))
		ocsd=np.abs(uncertainties.std_dev(self.mdOpenChCurrent))
		
		bcmu=np.abs(uncertainties.nominal_value(self.mdBlockedCurrent))
		bcsd=np.abs(uncertainties.std_dev(self.mdBlockedCurrent))

		# refine the start estimate
		idx=self.eStartEstimate

		try:
			while np.abs((np.abs(self.eventData[idx])-ocmu)/ocsd) > 5.0:
				idx-=1
			
			# Set the start point
			self.mdEventStart=idx+1

			# Next move the count forward so we are in the blocked channel region of the pulse
			while np.abs((np.abs(self.eventData[idx])-bcmu)/bcsd) > 0.5:
				idx+=1

			# Search for the event end. 7*sigma allows us to prevent false 
			# positives
			while np.abs((np.abs(self.eventData[idx])-bcmu)/bcsd) < 7.0:
				idx+=1

			# Finally backtrack to find the true event end
			while np.abs((np.abs(self.eventData[idx])-bcmu)/bcsd) > 5.0:
				idx-=1
		except ( IndexError, FloatingPointError ):
			self.rejectEvent('eResTime')
			return

		self.mdEventEnd=idx-1

		# residence time in ms
		self.mdResTime=1000.*((self.mdEventEnd-self.mdEventStart)/float(self.Fs))

    def __init__(self,centre, width, area, eta=0.25, asymmetry=0,background=0,**kwargs):

        #bypass __getattr__ - do this before anything else is created        
        object.__setattr__(self, '_pset', Parameters())

        self._pset['centre'] = Parameter(name='centre', value=nominal_value(centre), 
                                               vary=True, min=centre-2*width, max=centre+2*width, expr=None)
        self._pset['width'] = Parameter(name='width', value=nominal_value(width), 
                                              vary=True, 
                                              min=nominal_value(width)/3, 
                                              max=3*nominal_value(width), expr=None)
        self._pset['area'] = Parameter(name='area', value=nominal_value(area), 
                                             vary=True, 
                                             min=0, 
                                             max=None, expr=None)

        self._pset['eta'] = Parameter(name='eta', value=nominal_value(eta), 
                                            vary=True, min=-0.001, max=1.001, expr=None)
        
        self._pset['asymmetry'] = Parameter(name='asymmetry', value=nominal_value(asymmetry),
                                   vary=False, min=None, max=None, expr=None)
        
        #not currently used as a refinable parameter                           
        self._pset['background'] = Parameter(name='background', value=nominal_value(background),
                                   vary=False, min=0, max=None, expr=None)                           
       
        self.ymin_on_ymax = kwargs.get('ymin_on_ymax',0.001) 
        self.hkl = kwargs.get('hkl',None) 
        self.name = kwargs.get('name',None)  
        
        shape = kwargs.get('shape','pv')
        self.shape = shape.lower()        
    
        self.backup()

    def load_measurement(self, an, ar):

        # j = self._get_j(an)
        j = ar.j
        jf = 'NaN'
        if j is not None:
            jj = floatfmt(nominal_value(j), n=7, s=5)
            pe = format_percent_error(nominal_value(j), std_dev(j), include_percent_sign=True)
            jf = u'{} \u00b1{:0.2e}({})'.format(jj, std_dev(j), pe)

        a39 = ar.ar39decayfactor
        a37 = ar.ar37decayfactor
        ms = [MeasurementValue(name='Branch',
                               value=an.branch),
              MeasurementValue(name='DAQ Version',
                               value=an.collection_version),
              MeasurementValue(name='UUID',
                               value=an.uuid),
              MeasurementValue(name='RepositoryID',
                               value=an.repository_identifier),
              MeasurementValue(name='Spectrometer',
                               value=an.mass_spectrometer),
              MeasurementValue(name='Run Date',
                               value=an.rundate.strftime('%Y-%m-%d %H:%M:%S')),
              MeasurementValue(name='Irradiation',
                               value=self._get_irradiation(an)),
              MeasurementValue(name='J',
                               value=jf),
              MeasurementValue(name='Position Error',
                               value=floatfmt(an.position_jerr, use_scientific=True)),
              MeasurementValue(name='Lambda K',
                               value=nominal_value(ar.arar_constants.lambda_k),
                               units='1/a'),
              MeasurementValue(name='Project',
                               value=an.project),
              MeasurementValue(name='Sample',
                               value=an.sample),
              MeasurementValue(name='Material',
                               value=an.material),
              MeasurementValue(name='Comment',
                               value=an.comment),
              MeasurementValue(name='Ar39Decay',
                               value=floatfmt(a39)),
              MeasurementValue(name='Ar37Decay',
                               value=floatfmt(a37)),
              MeasurementValue(name='Sens.',
                               value=floatfmt(an.sensitivity, use_scientific=True),
                               units=an.sensitivity_units)]

        self.measurement_values = ms

    def _load_unknown_computed(self, an, new_list):
        attrs = (('Age', 'uage'),
                 # ('Age', 'age', None, None, 'age_err'),
                 ('w/o J', 'wo_j', '', 'uage', 'age_err_wo_j'),
                 ('K/Ca', 'kca'),
                 ('K/Cl', 'kcl'),
                 ('40Ar*', 'rad40_percent'),
                 ('F', 'uF'),
                 ('w/o Irrad', 'wo_irrad', '', 'uF', 'F_err_wo_irrad'))

        if new_list:
            def comp_factory(n, a, value=None, value_tag=None, error_tag=None):
                if value is None:
                    value = getattr(an, a)

                display_value = True
                if value_tag:
                    value = getattr(an, value_tag)
                    display_value = False

                if error_tag:
                    e = getattr(an, error_tag)
                else:
                    e = std_dev(value)

                return ComputedValue(name=n,
                                     tag=a,
                                     value=nominal_value(value) or 0,
                                     value_tag=value_tag,
                                     display_value=display_value,
                                     error=e or 0)

            cv = [comp_factory(*args)
                  for args in attrs]

            self.computed_values = cv
        else:
            for ci in self.computed_values:
                attr = ci.tag
                if attr == 'wo_j':
                    ci.error = an.age_err_wo_j
                    ci.value = nominal_value(getattr(an, ci.value_tag))
                elif attr == 'wo_irrad':
                    ci.error = an.F_err_wo_irrad
                    ci.value = nominal_value(getattr(an, ci.value_tag))
                else:
                    v = getattr(an, attr)
                    if v is not None:
                        ci.value = nominal_value(v)
                        ci.error = std_dev(v)

    def _update_ratios(self, an):
        for ci in self.computed_values:
            nd = ci.detectors
            n, d = nd.split('/')
            niso, diso = self._get_isotope(n), self._get_isotope(d)

            noncorrected = self._get_non_corrected_ratio(niso, diso)
            corrected, ic = self._get_corrected_ratio(niso, diso)

            ci.trait_set(value=floatfmt(nominal_value(corrected)),
                         error=floatfmt(std_dev(corrected)),
                         noncorrected_value=nominal_value(noncorrected),
                         noncorrected_error=std_dev(noncorrected),
                         ic_factor=nominal_value(ic))

    def _monte_carlo_error_propagation(self, vr, m):
        lambda_total = self._lambda_t
        el = self._lambda_ec
        f = self._f

        vel = nominal_value(el) + std_dev(el) * randn(self._n)
        vt = nominal_value(lambda_total) + std_dev(lambda_total) * randn(self._n)
        vf = nominal_value(f) + std_dev(f) * randn(self._n)

        vt_mc = ones(1, m) * vt
        vf_mc = ones(1, m) * vf
        vel_mc = ones(1, m) * vel

        t_mc = log(vt_mc / vel_mc * vf_mc * vr + 1) / vt_mc
        return mean(t_mc), std(t_mc)

    def get_mean_raw(self, tau=None):
        vs = []
        corrfunc = self._deadtime_correct
        for r in six.itervalues(self._cp):
            n = int(r['NShots'])
            nv = ufloat(float(r['Ar40']), float(r['Ar40err'])) * 6240
            dv = ufloat(float(r['Ar36']), float(r['Ar36err'])) * 6240
            if tau:
                dv = corrfunc(dv, tau * 1e-9)

            vs.append((n, nv / dv))

        key = lambda x: x[0]
        vs = sorted(vs, key=key)

        mxs = []
        mys = []
        mes = []
        for n, gi in groupby(vs, key=key):
            mxs.append(n)
            ys, es = list(zip(*[(nominal_value(xi[1]), std_dev(xi[1])) for xi in gi]))

            wm, werr = calculate_weighted_mean(ys, es)
            mys.append(wm)
            mes.append(werr)

        return mxs, mys, mes