def calibrate(self, title, rng=None, cal=None):
        """ Returns calibrated bin centers for saved histogram with given title.
        If the cal argument is given, title must exactly match a key in the file
        self.ofile.

        "title" title of histogram. If cal is not None then title must exactly
            match a key in the hdf5 file self.ofile.
        "rng" is the index range of elements to return
        "cal" is a set of polynomial coefficients for np.polyval. A linear
            scaling by a factor s would require cal=[s,0]
        """
        rng = rng or [None, None]
        if cal is None:
            bins = self.bcent(title, rng)
            return np.polyval(self.cal, bins)
        else:
            with h5py.File(self.ofile, 'r') as ofile:
                if title not in ofile.keys():
                    raise ValueError(title + ' is not in file ' + self.ofile)
                elif title[0] == 'b':
                    bins = self.bcent(title[1:], rng)
                    return np.polyval(cal, bins)
                else:
                    axis = ofile[title][rng[0]:rng[1]]
                    return np.polyval(cal, axis)

def showExamplePolyFit(xs,ys,fitDegree1 = 1,fitDegree2 = 2):
    pylab.figure()    
    pylab.plot(xs,ys,'r.',ms=2.0,label = "measured")

    # poly fit to noise
    coeeff = numpy.polyfit(xs, ys, fitDegree1)

    # Predict the curve
    pys = numpy.polyval(numpy.poly1d(coeeff), xs)

    se = mse(ys, pys)
    r2 = rSquared(ys, pys)

    pylab.plot(xs,pys, 'g--', lw=5,label="%d degree fit, SE = %0.10f, R2 = %0.10f" %(fitDegree1,se,r2))

    # Poly fit to noise
    coeeffs = numpy.polyfit(xs, ys, fitDegree2)

    # Predict the curve
    pys = numpy.polyval(numpy.poly1d(coeeffs), xs)

    se = mse(ys, pys)
    r2 = rSquared(ys, pys)

    pylab.plot(xs,pys, 'b--', lw=5,label="%d degree fit, SE = %0.10f, R2 = %0.10f" %(fitDegree2,se,r2))

    pylab.legend()

def cd_sphere_vector(Re):
    "Computes the drag coefficient of a sphere as a function of the Reynolds number Re."
    # Curve fitted after fig . A -56 in Evett & Liu :% " Fluid Mechanics & Hydraulics ",
    # Schaum ' s Solved Problems McGraw - Hill 1989.

    from numpy import log10,array,polyval
    CD = zeros_like(Re)
   
    CD = where(Re<0,0.0,0.0)     # condition 1
    
    CD = where((Re > 0.0) & (Re <=0.5),24/Re,CD) # condition 2

    p = array([4.22,-14.05,34.87,0.658])
    CD = where((Re > 0.5) & (Re <=100.0),polyval(p,1.0/Re),CD) #condition 3

    p = array([-30.41,43.72,-17.08,2.41])
    CD = where((Re >100.0)  & (Re <=1.0e4) ,polyval(p,1.0/log10(Re)),CD) #condition 4

    p = array([-0.1584,2.031,-8.472,11.932])
    CD = where((Re > 1.0e4)  &  (Re <=3.35e5),polyval(p,log10(Re)),CD) #condition 5

    CD = where((Re > 3.35e5) & (Re <=5.0e5),91.08*(log10(Re/4.5e5))**4 + 0.0764,CD) #condition 6

    p  = array([-0.06338,1.1905,-7.332,14.93])
    CD = where((Re > 5.05e5)  &  (Re <=8.0e6),polyval(p,log10(Re)),CD) #condition 7
    
    CD = where(Re>8.0e6,0.2,CD)  # condition 8

    return CD

    def _make_axes(self):
        '''Construct axes from calibration fields in header file
        '''
        xcalib = self.header.xcalibration
        ycalib = self.header.ycalibration

        xcalib_valid = struct.unpack('?', xcalib.calib_valid)

        if xcalib_valid:
            xcalib_order, = struct.unpack('>B', xcalib.polynom_order) # polynomial order
            px = xcalib.polynom_coeff[:xcalib_order+1]
            px = np.array(px[::-1]) # reverse coefficients to use numpy polyval
            pixels = np.arange(1, self.header.xdim + 1)
            px = np.polyval(px, pixels)
        else:
            px = np.arange(1, self.header.xdim + 1)

        ycalib_valid = struct.unpack('?', ycalib.calib_valid)

        if ycalib_valid:
            ycalib_order, = struct.unpack('>B', ycalib.polynom_order) # polynomial order
            py = ycalib.polynom_coeff[:ycalib_order+1]
            py = np.array(py[::-1]) # reverse coefficients to use numpy polyval
            pixels = np.arange(1, self.header.ydim + 1)
            py = np.polyval(py, pixels)
        else:
            py = np.arange(1, self.header.ydim + 1)

        self._xaxis = px
        self._yaxis = py

        return px, py

 def sampleR3(averagedist,boxdims):
     """low-discrepancy sampling using primes.
     The samples are evenly distributed with an average distance of averagedist inside the box with dimensions boxdims.
     Algorithim from "Geometric Discrepancy: An Illustrated Guide" by Jiri Matousek"""
     minaxis = numpy.argmin(boxdims)
     maxaxis = numpy.argmax(boxdims)
     meddimdist = numpy.sort(boxdims)[1]
     # convert average distance to number of samples.... do simple 3rd degree polynomial fitting...
     x = meddimdist/averagedist
     if x < 25.6:
         N = int(numpy.polyval([ -3.50181522e-01,   2.70202333e+01,  -3.10449514e+02, 1.07887093e+03],x))
     elif x < 36.8:
         N = int(numpy.polyval([  4.39770585e-03,   1.10961031e+01,  -1.40066591e+02, 1.24563464e+03],x))
     else:
         N = int(numpy.polyval([5.60147111e-01,  -8.77459988e+01,   7.34286834e+03, -1.67779452e+05],x))
     pts = numpy.zeros((N,3))
     pts[:,0] = numpy.linspace(0.0,meddimdist,N)
     pts[:,1] = meddimdist*numpy.mod(0.5+0.5*numpy.sqrt(numpy.arange(0,5.0*N,5.0)),1.0)
     pts[:,2] = meddimdist*numpy.mod(0.5+3*numpy.sqrt(numpy.arange(0,13.0*N,13.0)),1.0)
     if boxdims[minaxis] < meddimdist:
         pts = pts[pts[:,minaxis]<=boxdims[minaxis],:]
     if boxdims[maxaxis] > meddimdist:
         # have to copy across the max dimension
         numfullcopies = numpy.floor(boxdims[maxaxis]/meddimdist)
         oldpts = pts
         pts = numpy.array(oldpts)
         for i in range(int(numfullcopies)-1):
             oldpts[:,maxaxis] += meddimdist
             pts = numpy.r_[pts,oldpts]
         if boxdims[maxaxis]/meddimdist > numfullcopies:
             oldpts[:,maxaxis] += meddimdist
             pts = numpy.r_[pts,oldpts[oldpts[:,maxaxis]<=boxdims[maxaxis],:]]
     return pts

def corrNonlinGetPar(linearDet,nonLinearDet,order=2,data_0=0,
    correct_0=0,plot=False,returnCorrectedDet=False):
  """ Find parameters for non linear correction
    *linearDet* should be an 1D array of the detector that is linear
    *nonLinearDet* is the detector that is sussposed to be none linear
    *data_0" is an offset to use for the data (used only if plotting)"
    *correct_0* offset of the "linear detector"""
  p =  np.polyfit(nonLinearDet,linearDet,order)
  p[-1] = p[-1]-correct_0
  if plot:
    d = corrNonlin(nonLinearDet,p,data_0=data_0,correct_0=correct_0)
    plt.plot(linearDet,nonLinearDet,".",label="before correction")
    plt.plot(linearDet,d,".",label="after correction")
    poly_lin = np.polyfit(linearDet,d,1)
    xmin = min(linearDet.min(),0)
    xtemp = np.asarray( (xmin,linearDet.max()) )
    plt.plot(xtemp,np.polyval(poly_lin,xtemp),label="linear fit")
    plt.plot(linearDet,d-np.polyval(poly_lin,linearDet),
       ".",label="difference after-linear")
    plt.xlabel("linearDet")
    plt.ylabel("nonLinearDet")
    plt.legend()
  if order>=2 and p[-3]<0:
    log.warn("corrNonlinGetPar: consistency problem, second order coefficient should \
    be > 0, please double check result (plot=True) or try inverting the data and the\
    correct arguments")

  if returnCorrectedDet:
    return corrNonlin(nonLinearDet,p,data_0=data_0,correct_0=correct_0)
  else:
    return p

def plot_fit(deg, err):
    try:
        polcoefs = np.polyfit(x_rand, y_rand, deg)
    except np.RankWarning:
        pass
    if err:
        fig = plt.figure(figsize=(15, 10)) 
        gs = gridspec.GridSpec(2, 1, height_ratios=[3, 1])
        ax = plt.subplot(gs[0])
    else:
        fig, ax = plt.subplots(1, 1, figsize=(15, 10))
    # data
    ax.plot(x_cont, np.polyval(polcoefs, x_cont))
    ax.legend(['degree {}'.format(deg)])
    plot_dat(ax, err)
    # error
    if err:
        train_error.append(sum((np.polyval(polcoefs, x_rand)-y_rand)**2)/len(x_rand))
        test_error.append(sum((np.polyval(polcoefs, x_test)-y_test)**2)/len(x_test))
        ax = plt.subplot(gs[1])
        ax.plot(deg_list[0:len(train_error)], train_error, marker='o', c=(0, 0, 1))
        ax.plot(deg_list[0:len(test_error)], test_error, marker='o', c=(0.9, 0, 0))
        ax.set_xlabel('degree')
        ax.set_ylabel('error')
        ax.legend(['train', 'test'])
        ax.set_xlim(0.9, deg_list[-1]+0.1)
        ax.set_ylim(0, None)
        ax.set_yticks([])
        # layout
        plt.tight_layout()

    def update(self, rho):
        """
        Calculate the probability function for the given state of an harmonic
        oscillator (as density matrix)
        """

        if isket(rho):
            rho = ket2dm(rho)

        self.data = np.zeros(len(self.xvecs[0]), dtype=complex)
        M, N = rho.shape

        for m in range(M):
            k_m = pow(self.omega / pi, 0.25) / \
                sqrt(2 ** m * factorial(m)) * \
                exp(-self.xvecs[0] ** 2 / 2.0) * \
                np.polyval(hermite(m), self.xvecs[0])

            for n in range(N):
                k_n = pow(self.omega / pi, 0.25) / \
                    sqrt(2 ** n * factorial(n)) * \
                    exp(-self.xvecs[0] ** 2 / 2.0) * \
                    np.polyval(hermite(n), self.xvecs[0])

                self.data += np.conjugate(k_n) * k_m * rho.data[m, n]

def rssmodelwave(grating,grang,artic,cbin,cols):
#   compute wavelengths from model (this can probably be done using pyraf spectrograph model)
    spec=np.loadtxt(datadir+"spec.txt",usecols=(1,))
    Grat0,Home0,ArtErr,T2Con,T3Con=spec[0:5]
    FCampoly=spec[5:11]
    grname=np.loadtxt(datadir+"gratings.txt",dtype=str,usecols=(0,))
    grlmm,grgam0=np.loadtxt(datadir+"gratings.txt",usecols=(1,2),unpack=True)

    grnum = np.where(grname==grating)[0][0]
    lmm = grlmm[grnum]
    alpha_r = np.radians(grang+Grat0)
    beta0_r = np.radians(artic*(1+ArtErr)+Home0)-alpha_r
    gam0_r = np.radians(grgam0[grnum])
    lam0 = 1e7*np.cos(gam0_r)*(np.sin(alpha_r) + np.sin(beta0_r))/lmm
    ww = lam0/1000. - 4.
    fcam = np.polyval(FCampoly,ww)
    disp = (1e7*np.cos(gam0_r)*np.cos(beta0_r)/lmm)/(fcam/.015)
    dfcam = 3.162*disp*np.polyval([FCampoly[x]*(5-x) for x in range(5)],ww)
    T2 = -0.25*(1e7*np.cos(gam0_r)*np.sin(beta0_r)/lmm)/(fcam/47.43)**2 + T2Con*disp*dfcam
    T3 = (-1./24.)*3162.*disp/(fcam/47.43)**2 + T3Con*disp
    T0 = lam0 + T2 
    T1 = 3162.*disp + 3*T3
    X = (np.array(range(cols))+1-cols/2)*cbin/3162.
    lam_X = T0+T1*X+T2*(2*X**2-1)+T3*(4*X**3-3*X)
    return lam_X

def get_cont(x,y,n=1,sl=1.,sh=5.):
	orilen = len(x)
	coef = np.polyfit(x,y,n)
	res = y - np.polyval(coef,x)
	IH = np.where(res>0)[0]
	IL = np.where(res<0)[0]
	dev = np.mean(res[IH])
	I = np.where((res>-sl*dev) & (res<sh*dev))[0]
	J1 = np.where(res<=-sl*dev)[0]
	J2 = np.where(res>=sh*dev)[0]
	J = np.unique(np.hstack((J1,J2)))
	cond = True
	if len(J)==0 or len(x)< .3*orilen:
		cond=False
	while cond:
		x = np.delete(x,J)
		y = np.delete(y,J)
		coef = np.polyfit(x,y,n)
		res = y - np.polyval(coef,x)
		IH = np.where(res>0)[0]
		IL = np.where(res<0)[0]
		dev = np.mean(res[IH])
		I = np.where((res>-sl*dev) & (res<sh*dev))[0]
		J1 = np.where(res<=-sl*dev)[0]
		J2 = np.where(res>=sh*dev)[0]
		J = np.unique(np.hstack((J1,J2)))
		cond = True
		if len(J)==0 or len(x)< .1*orilen:
			cond=False
	return coef

def get_ratio(sciw,rat,n=3):
	rat = scipy.signal.medfilt(rat,11)
	lori = len(sciw)
	coef = np.polyfit(sciw,rat,n)
	res = rat - np.polyval(coef,sciw)
	rms = np.sqrt(np.mean(res**2))
	I = np.where(res> 3*rms)[0]
	I2 = np.where(res< -3*rms)[0]
	I = np.sort(np.hstack((I,I2)))
	cond = True
	if len(I) == 0 or len(sciw) < .3 * lori:
		cond = False

	while cond:
		#imax = np.argmax(res**2)
		#sciw = np.delete(sciw,imax)
		#rat  = np.delete(rat,imax)
		sciw = np.delete(sciw,I)
		rat  = np.delete(rat,I)
		coef = np.polyfit(sciw,rat,n)
		res = rat - np.polyval(coef,sciw)
		rms = np.sqrt(np.mean(res**2))
		I = np.where(res> 3*rms)[0]
		I2 = np.where(res< -3*rms)[0]
		I = np.sort(np.hstack((I,I2)))
		if len(I) == 0 or len(sciw) < .3 * lori:
			cond = False

	return coef

def get_rats(ZO,ZI,ZF,pars):
	ords = []
	for i in range(sc.shape[1]):
		J1 = np.where(mw > sc[0,i,-1])[0]
		J2 = np.where(mw < sc[0,i,0])[0]
		if len(J1)>0 and len(J2)>0:
			ords.append(i)
	ords = np.array(ords)
	mf = get_full_model(pars[0],pars[1],pars[2],pars[3],RES_POW)
	tmodf = np.zeros((sc.shape[1],sc.shape[2]))
	tscif = np.zeros((sc.shape[1],sc.shape[2]))
	test_plot = np.zeros((4,sc.shape[1],sc.shape[2]))
	for i in ords:
		I = np.where((mw>sc[0,i,0]) & (mw<sc[0,i,-1]))[0]
		modw = mw[I]
		modf = mf[I]
		sciw = sc[0,i]
		scif = sc[3,i]/np.median(sc[3,i])
		modf = pixelization(modw,modf,sciw)
		#IMB = np.where(mask_bin[i]!=0)[0]
		#modf /= modf[IMB].mean()
		mscif = scipy.signal.medfilt(scif,11)
		rat = modf/mscif
		INF = np.where(mscif!=0)[0]
		coef = get_ratio(sciw[INF],rat[INF])
		scif = scif * np.polyval(coef,sciw)
		mscif = mscif * np.polyval(coef,sciw)
		coef = get_cont(sciw,mscif)
		scif = scif / np.polyval(coef,sciw)
		#plot(sciw,scif)
		coef = get_cont(sciw,modf)
		modf = modf / np.polyval(coef,sciw)
		#plot(sciw,modf)	
		tmodf[i] = modf
		tscif[i] = scif
		test_plot[0,i] = sc[0,i]
		test_plot[1,i] = scif
		test_plot[2,i] = modf
		test_plot[3,i] = mask_bin[i]
	#show()
	#print vcdx
	hdu = pyfits.PrimaryHDU(test_plot)
	os.system('rm example.fits')
	hdu.writeto('example.fits')
	rat = tscif/tmodf

	nejx = np.arange(100)/100.
	ratsout = []

	for i in range(len(ZI)):
		ejy = rat[ZO[i],ZI[i]:ZF[i]]
		ejx = np.arange(len(ejy))/float(len(ejy))
		tck = interpolate.splrep(ejx,ejy,k=3)
		if len(ratsout)==0:
			ratsout = interpolate.splev(nejx,tck)
		else:
			ratsout = np.vstack((ratsout,interpolate.splev(nejx,tck)))
			#plot(interpolate.splev(nejx,tck))
	#show()
	return ratsout

def detrend(var, ax=None, lcopy=True, ldetrend=True, ltrend=False, degree=1, rcond=None, w=None,  
            lsmooth=False, lresidual=False, window_len=11, window='hanning'): 
  ''' subtract a linear trend from a time-series array (operation is in-place) '''
  # check input
  if not isinstance(var,np.ndarray): raise NotImplementedError # too many checks
  if lcopy: var = var.copy() # make copy - not in-place!
  # fit over entire array (usually not what we want...)
  if ax is None and ldetrend: ax = np.arange(var.size) # make dummy axis, if necessary
  if var.ndim != 1:
    shape = var.shape 
    var = var.ravel() # flatten array, if necessary
  else: shape = None
  # apply optional detrending
  if ldetrend or ltrend:
    # fit linear trend
    trend = np.polyfit(ax, var, deg=degree, rcond=rcond, w=w, full=False, cov=False)
    # evaluate and subtract linear trend
    if ldetrend and ltrend: raise ArgumentError("Can either return trend/polyfit or residuals, not both.")
    elif ldetrend and not ltrend: var -= np.polyval(trend, ax) # residuals
    elif ltrend and not ldetrend: var = np.polyval(trend, ax) # residuals
  # apply optional smoothing
  if lsmooth and lresidual: raise ArgumentError("Can either return smoothed array or residuals, not both.")
  elif lsmooth: var = smooth(var, window_len=window_len, window=window)  
  elif lresidual: var -= smooth(var, window_len=window_len, window=window)
  # return detrended and/or smoothed time-series
  if shape is not None: var = var.reshape(shape)
  return var

def dualPlot(age, meanWithin, meanBetween, title):

    fig, (within, between) = plt.subplots(1, 2, sharex=True, sharey=False)

    # fitshit
    wP = np.polyfit(age, meanWithin, 1)
    bP = np.polyfit(age, meanBetween, 1)
    xnew = np.arange(age.min() - 1, age.max() + 1, 0.1)
    wFit = np.polyval(wP, xnew)
    bFit = np.polyval(bP, xnew)

    within.set_title("within network")
    between.set_title("between network")

    withinCorr, withinP = st.pearsonr(age, meanWithin)
    within.plot(age, meanWithin, "k.")
    within.plot(xnew, wFit, "r", label=(str(np.round(withinCorr, 2)) + " " + str(np.round(withinP, 4))))
    within.set_xlabel("mean connectivity")
    within.set_ylabel("age")
    within.legend()

    betweenCorr, betweenP = st.pearsonr(age, meanBetween)
    between.plot(age, meanBetween, "k.")
    between.plot(xnew, bFit, "b", label=(str(np.round(betweenCorr, 2)) + " " + str(np.round(betweenP, 4))))
    between.set_xlabel("mean connectivity")
    between.set_ylabel("age")
    between.legend()

    fig.suptitle(title)
    plt.show()
    raw_input("Press Enter to continue...")
    plt.close()

def tfit(line):
    """
    Correct for temperature systematics.  Fit a polynomial to (teff,abund)
    and require that the corrected solar value be 0.  We cut on vsini, 

    returns:
    (fitabund,fitpar,t,abund)
    fitabund - the temperature corrected abundance
    fitpar   - the parameters to the polynomial fit
    t        - the temperature array
    abund    - the non-temp-corrected abundances

    """
    deg  = 3 # fit with a 3rd degree polynomial
    #define abundance for the particular line we're looking at
    p = getelnum.Getelnum(line)
    elstr = p.elstr

    conn = sqlite3.connect(os.environ['STARSDB'])
    cur = conn.cursor()

    #pull in the abundances and the non-corrected abundances
    cmd = 'SELECT '+elstr+'_abund_nt,teff FROM mystars WHERE '+globcut(elstr)
    cur.execute(cmd)
    arr = np.array(cur.fetchall() ) 
    abund,t = arr[:,0],arr[:,1]
    abund = abund - p.abnd_sol

    #fit the points
    fitpar = np.polyfit(t,abund,deg)
    #subtract out the fit, while requiring that the solar value be 0.
    fitpar[deg] = fitpar[deg] - np.polyval(fitpar,p.teff_sol)
    fitabund = abund - np.polyval(fitpar,t)
    return (fitabund,fitpar,t,abund)

    def unravel_box(self, box):
        """Convert a rectangular represenation of the spectra back to a single
           array

        Parameters
        ----------
        box: ~numpy.ndarray
            Rectangular represnation of flux

        Returns
        -------
        data: ~numpy.ndarray
            Array of values to convert into a rectangular representation
        """
        xmax = self.region[1].max()
        xmin = 0
        ymax = self.region[0].max()
        ymin = self.region[0].min()
        ys = ymax-ymin
        xs = xmax-xmin
        data = np.zeros((ys+1,xs+1))
        coef = np.polyfit(self.region[1], self.region[0], 3)
        xarr = np.arange(xs+1)
        yarr = np.polyval(coef, xarr)-ymin
        x = self.region[1]-xmin
        y = self.region[0]-ymin - (np.polyval(coef, x) - ymin - yarr.min()).astype(int)
        data = np.zeros(self.npixels)
        data = box[y,x]
        return data

def measure_dA_dphi_fir(lock, li, tp, dA_dphi_before, dA_dphi_after):
    """Correct for impulsive phase shift at end of pulse time."""

    i_tp = np.arange(lock.t.size)[lock.t < tp][-1]
    # Use 20 data points for interpolating; this is slightly over one
    # cycle of our oscillation
    m = np.arange(-10, 11) + i_tp
    # This interpolator worked reasonably for similar, low-frequency sine waves
    interp = interpolate.KroghInterpolator(lock.t[m], lock.x[m])
    x0 = interp(tp)[()]
    # We only need t0 approximately; the precise value of f0 doesn't matter very much.
    t0 = li.t[(li.t < tp)][-1]
    f0 = li.df[(li.t < tp)][-1] + li.f0(t0)
    v0 = interp.derivative(tp)[()]
    x2 = v0 / (2*np.pi*f0)

    phi0 = np.arctan2(-x2, x0)

    ml = masklh(li.t, tp-t_fit, tp)
    mr = masklh(li.t, tp, tp + t_fit)

    A = abs(li.z_out)
    phi = np.unwrap(np.angle(li.z_out))/(2*np.pi)

    mbAl = np.polyfit(li.t[ml], A[ml], 1)
    mbAr = np.polyfit(li.t[mr], A[mr], 1)

    mb_phi_l = np.polyfit(li.t[ml], phi[ml], 1)
    mb_phi_r = np.polyfit(li.t[mr], phi[mr], 1)

    dA = np.polyval(mbAr, tp) - np.polyval(mbAl, tp)
    dphi = np.polyval(mb_phi_r, tp) - np.polyval(mb_phi_l, tp)

    return phi0, dA, dphi

def eeval(expression, w):
	""" evaluate a sympy expression at omega. return magnitude, phase."""
	num, den = e2nd(expression)
	y = numpy.polyval(num, 1j*w) / numpy.polyval(den, 1j*w)
	phase = numpy.arctan2(y.imag, y.real) * 180.0 / numpy.pi
	mag = abs(y)
	return mag, phase

    def set_limits_by_zoom(self, zoom, cx, cy, canvas=None):
        '''
        '''
        def _set_limits(axis_key, px_per_cm, cur_pos, canvas):

            if axis_key == 'x':
                d = self.width
            else:
                d = self.height

            # scale to mm
            if canvas is None:
                canvas = self.parent

            if canvas:
                d /= 2.0 * px_per_cm / 10.0
                lim = (-d + cur_pos, d + cur_pos)
                canvas.set_mapper_limits(axis_key, lim)

        cdata = self.calibration_data
        xpx_per_cm = np.polyval(cdata.get_xcoeffs(), [zoom])[0]
        ypx_per_cm = np.polyval(cdata.get_ycoeffs(), [zoom])[0]

        _set_limits('x', xpx_per_cm, cx, canvas)
        _set_limits('y', ypx_per_cm, cy, canvas)

def _precess_from_J2000_Capitaine(epoch):
    """
    Computes the precession matrix from J2000 to the given Julian Epoch.
    Expression from from Capitaine et al. 2003 as expressed in the USNO
    Circular 179.  This should match the IAU 2006 standard from SOFA.

    Parameters
    ----------
    epoch : scalar
        The epoch as a julian year number (e.g. J2000 is 2000.0)

    """
    from .angles import rotation_matrix

    T = (epoch - 2000.0) / 100.0
    # from USNO circular
    pzeta = (-0.0000003173, -0.000005971, 0.01801828, 0.2988499, 2306.083227, 2.650545)
    pz = (-0.0000002904, -0.000028596, 0.01826837, 1.0927348, 2306.077181, -2.650545)
    ptheta = (-0.0000001274, -0.000007089, -0.04182264, -0.4294934, 2004.191903, 0)
    zeta = np.polyval(pzeta, T) / 3600.0
    z = np.polyval(pz, T) / 3600.0
    theta = np.polyval(ptheta, T) / 3600.0

    return rotation_matrix(-z, 'z') *\
           rotation_matrix(theta, 'y') *\
           rotation_matrix(-zeta, 'z')