def plot_vector_diff_F160W():
    """
    Using the xym2mat analysis on the F160W filter in the 2010 data set,
    plot up the positional offset vectors for each reference star in
    each star list relative to the average list. This should show
    us if there are systematic problems with the distortion solution
    or PSF variations.
    """
    x, y, m, xe, ye, me, cnt = load_catalog()

    dx = x - x[0]
    dy = y - y[0]
    dr = np.hypot(dx, dy)

    py.clf()
    py.subplots_adjust(left=0.05, bottom=0.05, right=0.95, top=0.95)

    for ii in range(1, x.shape[0]):
        idx = np.where((x[ii,:] != 0) & (y[ii,:] != 0) & (dr[ii,:] < 0.1) &
                       (xe < 0.05) & (ye < 0.05))[0]

        py.clf()
        q = py.quiver(x[0,idx], y[0,idx], dx[ii, idx], dy[ii, idx], scale=1.0)
        py.quiverkey(q, 0.9, 0.9, 0.03333, label='2 mas', color='red')
        py.title('{0} stars in list {1}'.format(len(idx), ii))
        py.xlim(0, 4500)
        py.ylim(0, 4500)

        foo = input('Continue?')
        if foo == 'q' or foo == 'Q':
            break

def check_vpd_ks2_astrometry():
    """
    Check the VPD and quiver plots for our KS2-extracted, re-transformed astrometry.
    """
    catFile = workDir + '20.KS2_PMA/wd1_catalog.fits'
    tab = atpy.Table(catFile)

    good = (tab.xe_160 < 0.05) & (tab.ye_160 < 0.05) & \
        (tab.xe_814 < 0.05) & (tab.ye_814 < 0.05) & \
        (tab.me_814 < 0.05) & (tab.me_160 < 0.05)

    tab2 = tab.where(good)

    dx = (tab2.x_160 - tab2.x_814) * ast.scale['WFC'] * 1e3
    dy = (tab2.y_160 - tab2.y_814) * ast.scale['WFC'] * 1e3

    py.clf()
    q = py.quiver(tab2.x_814, tab2.y_814, dx, dy, scale=5e2)
    py.quiverkey(q, 0.95, 0.85, 5, '5 mas', color='red', labelcolor='red')
    py.savefig(workDir + '20.KS2_PMA/vec_diffs_ks2_all.png')

    py.clf()
    py.plot(dy, dx, 'k.', ms=2)
    lim = 30
    py.axis([-lim, lim, -lim, lim])
    py.xlabel('Y Proper Motion (mas)')
    py.ylabel('X Proper Motion (mas)')
    py.savefig(workDir + '20.KS2_PMA/vpd_ks2_all.png')

    idx = np.where((np.abs(dx) < 10) & (np.abs(dy) < 10))[0]
    print('Cluster Members (within dx < 10 mas and dy < 10 mas)')
    print(('   dx = {dx:6.2f} +/- {dxe:6.2f} mas'.format(dx=dx[idx].mean(),
                                                        dxe=dx[idx].std())))
    print(('   dy = {dy:6.2f} +/- {dye:6.2f} mas'.format(dy=dy[idx].mean(),
                                                        dye=dy[idx].std())))

def plot_residuals_year_2pos(year, filt, pos1, pos2, refresh=False):
    year = str(year) 
    
    dir_xym = year + '_' + filt + '/01.XYM/'

    # Try to read from a FITS table of compiled data if possible.
    if refresh == True:
        make_residuals_table_year_2pos(year, filt, pos1, pos2)

    d = read_residuals_table_year_2pos(year, filt, pos1, pos2)

    print d.xstd_p.min(), d.xstd_p.max()

    # Remember, this was all aligned to F814W in ACS. So there is a
    # relative scale change between the transformed and the raw pixels.
    # We need to put evertyhing back on the same scale.
    # The plate scale was 0.411095 (pulled from TRANS.xym2mat.ref4)
    scale = 0.411095

    dx1 = d.xmean_p[0, :] - d.stars['x']
    dx2 = d.xmean_p[1, :] - d.stars['x']
    dy1 = d.ymean_p[0, :] - d.stars['y']
    dy2 = d.ymean_p[1, :] - d.stars['y']
    lim = 0.06
    idx = np.where((np.abs(dx1) < lim) & (np.abs(dy1) < lim) & 
                   (np.abs(dx2) < lim) & (np.abs(dy2) < lim))[0]

    # Plot the offsets for each position on a common coordinate system.
    py.clf()
    q1 = py.quiver(d.xmean[idx], d.ymean[idx], dx1[idx], dy1[idx], 
                   color='red', scale=1.8)
    q1 = py.quiver(d.xmean[idx], d.ymean[idx], dx2[idx], dy2[idx], 
                   color='blue', scale=1.8)
    py.quiverkey(q1, -0.1, -0.12, 0.02/scale, '0.02 WFC3IR pixels')
    py.axis('equal')
    py.xlabel("X (pixels)")
    py.ylabel("Y (pixels)")
    py.title(year + ',' + filt + ',' + pos1 + ',' + pos2)
    py.savefig('plots/resid_final_{0}_{1}_{2}_{3}.png'.format(year, filt, pos1, pos2))

    # Plot the offsets for each position on the raw coordinate system.
    py.clf()
    xraw1 = d.xraw[0,:,:].mean(axis=0)
    yraw1 = d.yraw[0,:,:].mean(axis=0)
    xraw2 = d.xraw[1,:,:].mean(axis=0)
    yraw2 = d.yraw[1,:,:].mean(axis=0)
    
    q1 = py.quiver(xraw1[idx], yraw1[idx], dx1[idx], dy1[idx], 
                   color='red', scale=1.8)
    q1 = py.quiver(xraw2[idx], yraw2[idx], dx2[idx], dy2[idx], 
                   color='blue', scale=1.8)
    py.quiverkey(q1, -0.1, -0.12, 0.02/scale, '0.02 WFC3IR pixels')
    py.axis('equal')
    py.xlabel("X (pixels)")
    py.ylabel("Y (pixels)")
    py.title(year + ',' + filt + ',' + pos1 + ',' + pos2)
    py.savefig('plots/resid_raw_{0}_{1}_{2}_{3}.png'.format(year, filt, pos1, pos2))

    def PlotBaseFly(self):
        self.BaseFly()
        # Creates lists of x and y coordinates
        x_coords = [0, (self.displ * numpy.cos(self.head * (2 * numpy.pi) / 360))]
        y_coords = [0, (self.displ * numpy.sin(self.head * (2 * numpy.pi) / 360))]

        fig = pylab.figure()
        # Adjusts margins
        fig.subplots_adjust(left=None, bottom=0.2, right=None, top=None, wspace=None, hspace=None)

        # Plots start point
        pylab.plot(x_coords[0], y_coords[0], self.marker1, markersize=10, label='Start')
        # Plots end point
        pylab.plot(x_coords[1], y_coords[1], self.marker2, markersize=10, label='End')

        # Generates 2D wind vector field
        if self.displ < 10:
            X, Y = numpy.meshgrid(numpy.arange(-30, 30, 6), numpy.arange(-30, 30, 6))
            Q = pylab.quiver(X, Y, self.wind_U, self.wind_V, label='Wind Direction')
            QK = pylab.quiverkey(Q, 0.1, -0.08, 10, 'Wind Direction', labelpos='W',
                                 fontproperties={'weight': 'bold'})
        else:
            X, Y = numpy.meshgrid(numpy.arange(-self.displ - 20, self.displ + 20, self.displ / 5), numpy.arange(-self.displ - 20, self.displ + 20, self.displ / 5))
            Q = pylab.quiver(X, Y, self.wind_U, self.wind_V, label='Wind Direction')
            QK = pylab.quiverkey(Q, 0.1, -0.08, self.displ / 5, 'Wind Direction', labelpos='W',
                                 fontproperties={'weight': 'bold'})

        # Shows displ, head, time
        # Refer to http://docs.python.org/2/library/string.html for formatting
        # and http://stackoverflow.com/questions/19261089/how-to-use-new-style-string-formatting-in-matplotlib-figure-text
        '''
        Text = '{0} {1:10}\n'.format('Displacement: ', int(self.displ)) + \
               '{0} {1:10}\n'.format('Heading:         ', int(self.head)) + \
               '{0} {1:10}\n'.format('Flight Time:  ', int(self.time))

        fig.text(0.25, -0.01, Text, family='monospace')
        '''

        TextL = 'Displacement (m):\nHeading (Degrees):\nDescent Time (s):'
        TextR = str(int(self.displ)) + '\n' + str(int(self.head)) + '\n' + str(int(self.time))
        fig.text(0.4, 0.01, TextL, multialignment='left')
        fig.text(0.6, 0.01, TextR, multialignment='right')

        pylab.title('UAV Displacement in Y against UAV Displacement in X')
        pylab.xlabel('UAV Displacement in X')
        pylab.ylabel('UAV Displacement in Y')
        pylab.xlim(-self.displ - 20, self.displ + 20)
        pylab.ylim(-self.displ - 20, self.displ + 20)

        pylab.legend()

        pylab.show()

def show_shwfs_vecs(
    shiftvec, subaps, refpos=None, img=None, extent=None, title=None, scale=10, pause=False, fignum=None, keep=True
):
    """
	Show a SHWFS measurements with subapertures and potentially a background 
	image. If **refpos** is given, SHWFS shift vectors will be plotted at 
	those locations, otherwise at the center of the sub apertures.

	Additionally, a background image (wavefront, sensor image) can be 
	plotted along with the measurements in the background. If the coordinate 
	system is different than that of the sub apertures, this can be fixed by 
	using the **extent** keyword.

	**Scale** tweaks the arrow length, unity being the approximate real size.
	The default is 10 to enhance the visibility.

	@param [in] shiftvec (N,2) vector with shifts
	@param [in] subaps (N,4) vector with sub aperture positions as (low0, high0, low1, high1)
	@param [in] refpos Positions to plot shift vectors at (if None use subaps center)
	@param [in] img Background image, e.g. a wavefront
	@param [in] extent Extent for the image for imshow()
	@param [in] title Plot title
	@param [in] scale Scale used for arrow width, unity is approximately
	@param [in] pause Pause before continuing
	@param [in] fignum Figure number to use for figure()
	@param [in] keep Keep plot window open after returning
	"""

    import pylab as plt
    from matplotlib.collections import PatchCollection

    sasize = np.r_[subaps[:, :2].ptp(1).mean(), subaps[:, 2:].ptp(1).mean()]
    subaps_im = [
        plt.Rectangle((subap[1], subap[0]), sasize[0], sasize[1], fc="none", ec="k") for subap in subaps[:, ::2]
    ]

    plt.figure(fignum)
    plt.clf()
    plt.title(title or "SHWFS vectors on MLA grid")

    if img != None:
        plt.imshow(img, extent=extent)

    if refpos == None:
        refpos = subaps[:, ::2] + sasize / 2

    q = plt.quiver(
        refpos[:, 1], refpos[:, 0], shiftvec[:, 1], shiftvec[:, 0], angles="xy", scale=refpos.ptp() / 10.0, color="r"
    )
    p = plt.quiverkey(q, refpos.min(0)[1], refpos.min(0)[0], 5, "5 pix.", coordinates="data", color="r")

    thisgca = plt.gca()
    thisgca.add_collection(PatchCollection(subaps_im, match_original=True))

    if pause:
        raw_input("Press any key to continue...")
    if not keep:
        plt.close()

def test_darPlusDistortion():
    data_dir = module_dir + '/distortion/'
    file_geox_darunfix = data_dir + 'nirc2dist_xgeoim.fits'
    file_geoy_darunfix = data_dir + 'nirc2dist_ygeoim.fits'

    data_dir = '/u/ghezgroup/data/m92_test/08jul_new_on/'
    file_geox_darfix = data_dir + 'reduce/kp/gc_f1/ce0249geo_x.fits'
    file_geoy_darfix = data_dir + 'reduce/kp/gc_f1/ce0249geo_y.fits'

    xon = pyfits.getdata(file_geox_darfix)
    yon = pyfits.getdata(file_geoy_darfix)
    xoff = pyfits.getdata(file_geox_darunfix)
    yoff = pyfits.getdata(file_geoy_darunfix)

    # Make arrays with the coordinates for each 
    imgsize = 1024
    axisX = np.arange(imgsize, dtype=float)
    axisY = np.arange(imgsize, dtype=float)
    xcoo2d, ycoo2d = np.meshgrid(axisX, axisY)

    # Lets trim so that we only keep every 20th pixel
    idx = np.arange(25, imgsize, 50)
    xon = xon.take(idx, axis=0).take(idx, axis=1)
    yon = yon.take(idx, axis=0).take(idx, axis=1)
    xoff = xoff.take(idx, axis=0).take(idx, axis=1)
    yoff = yoff.take(idx, axis=0).take(idx, axis=1)
    xcoo2d = xcoo2d.take(idx, axis=0).take(idx, axis=1)
    ycoo2d = ycoo2d.take(idx, axis=0).take(idx, axis=1)

    # Calculate differences
    xdiff = xon - xoff
    ydiff = yon - yoff

    # Make vector plots
    py.clf()
    qvr = py.quiver2([xcoo2d], [ycoo2d], [xdiff], [ydiff],
                     units='width', scale=5, 
                     width=0.005, headwidth=3, headlength=3, 
                     headaxislength=3)
    py.quiverkey(qvr, 100, 1120, 1.0, '1 pixel', coordinates='data', color='r')
    py.xlabel('NIRC2 X (pixel)')
    py.ylabel('NIRC2 Y (pixel)')
    py.title('Arrows point to DAR Fix')
    #py.savefig('plots/vector_daroffon.png')
    py.show()

def plot_arches_pm():
    """
    Plot the proper motion vectors for the Arches (Clarkson+ 2012).
    """
    dataFile = "/u/jlu/data/arches/clarkson2012_table5.txt"

    d = atpy.Table(dataFile, type="mrt")

    # Put into a "field" reference frame.
    idx = np.where(d.FMem > 0.9)[0]
    f_pmX = d.pmX[idx].mean()
    f_pmY = d.pmY[idx].mean()

    d.pmX -= f_pmX
    d.pmY -= f_pmY

    py.clf()
    Q = py.quiver(d.DelX, d.DelY, d.pmX, d.pmY)
    py.quiverkey(Q, 95, 95, 1, "40 km/s", coordinates="figure")

def plot_vpd_align(alignRoot, epoch1, epoch2, outRoot='plot_align'):
    """
    Read in an aligned data set. Plot a vector-point-diagram for the 
    positional differences between two different epochs within the data set.
    """
    s = starset.StarSet(alignRoot)

    x1 = s.getArrayFromEpoch(epoch1, 'xorig')
    y1 = s.getArrayFromEpoch(epoch1, 'yorig')
    x2 = s.getArrayFromEpoch(epoch2, 'xorig')
    y2 = s.getArrayFromEpoch(epoch2, 'yorig')

    x1e = s.getArrayFromEpoch(epoch1, 'xpixerr_p')
    y1e = s.getArrayFromEpoch(epoch1, 'ypixerr_p')
    x2e = s.getArrayFromEpoch(epoch2, 'xpixerr_p')
    y2e = s.getArrayFromEpoch(epoch2, 'ypixerr_p')

    dx = x2 - x1
    dy = y2 - y1
    dxe = np.hypot(x1e, x2e)
    dye = np.hypot(y1e, y2e)
    dr = np.hypot(dx, dy)

    idx = np.where((dxe < 0.02) & (dye < 0.02) & (dr < 3.0))[0]
    print len(x1), len(idx)

    drmax = 1.5

    py.figure(1)
    py.clf()
    py.plot(dx[idx], dy[idx], 'k.', ms=2)
    py.xlabel('Delta X (pixels)')
    py.ylabel('Delta Y (pixels)')
    py.axis([-drmax, drmax, -drmax, drmax])
    py.savefig(outRoot + '_vpd.png')

    py.figure(2)
    py.clf()
    q = py.quiver(x1[idx], y1[idx], dx[idx], dy[idx], scale=10)
    py.quiverkey(q, 0.5, 0.95, 0.05, '0.05 pix', color='red')
    py.savefig(outRoot + '_dr_vec.png')

def whisker4QReduce(X2WIN_IMAGE=None,Y2WIN_IMAGE=None,XYWIN_IMAGE=None):
    """
    This function make the whisker plot based on the sextractor export from QuickReduce
    J.Hao, 12/4/2012
    """
    xpos = np.genfromtxt('xpos_ypos_fp.txt').T[0]
    ypos = np.genfromtxt('xpos_ypos_fp.txt').T[1]
    temp = np.zeros(62).astype('complex')
    temp.real = X2WIN_IMAGE - Y2WIN_IMAGE
    temp.imag = 2*XYWIN_IMAGE
    data=np.array([xpos, ypos,X2WIN_IMAGE + Y2WIN_IMAGE,temp]).T
    data = averageN30new(data)
    data = subMeanAll(data)
    pl.figure(figsize=(11,5.5))
    pl.subplot(1,2,1)
    phi22 = 0.5*np.arctan2(data[:,3].imag,data[:,3].real)
    x = data[:,0].real
    y = data[:,1].real
    u = np.abs(data[:,3])*np.cos(phi22)
    v = np.abs(data[:,3])*np.sin(phi22)
    qvr = pl.quiver(x,y,u,v,width = 0.004, color='r',pivot='middle',headwidth=0.,headlength=0.,headaxislength=0.,scale_units='width')
    qk = pl.quiverkey(qvr, -150,-240,0.3,str(0.3)+' pix^2',coordinates='data',color='blue')
    pl.plot(x,y,'b,')
    pl.xlim(-250,250)
    pl.ylim(-250,250)
    pl.grid(color='g')
    pl.xlabel('Camera WEST [mm]')
    pl.ylabel('Camera NORTH [mm]')
    pl.title('M22')
    pl.subplot(1,2,2)
    m20sqr = np.sqrt(data[:,2].real)
    x = data[:,0].real
    y = data[:,1].real
    m20sqr_med = np.median(m20sqr) 
    m20sqr_diff = m20sqr - m20sqr_med
    m20sqr_diff_absmed = np.median(np.abs(m20sqr_diff))
    plotScale = 1./m20sqr_diff_absmed*100
    pos = m20sqr_diff >=0
    neg = m20sqr_diff < 0
    pl.scatter(x[pos],y[pos],s=m20sqr_diff[pos]*plotScale,c='r',alpha=0.5)
    pl.scatter(x[neg],y[neg],s=-m20sqr_diff[neg]*plotScale,c='b',alpha=0.5)
    pl.scatter(-230,-210,s=0.01*plotScale,c='b',alpha=0.5)
    pl.text(-200,-215,'-'+str(0.01)+' pix')
    pl.scatter(-230,-230,s=0.01*plotScale,c='r',alpha=0.5)
    pl.text(-200,-235,str(0.01)+' pix')
    pl.plot(x,y,'y,')
    pl.grid(color='g')
    pl.xlim(-250,250)
    pl.ylim(-250,250)
    pl.xlabel('Camera WEST [mm]')
    pl.ylabel('Camera NORTH [mm]')
    pl.title('median '+r'$\sqrt{M20}$: '+str(round(0.27*m20sqr_med,3))+' [arcsec]')
    return '---done!--'

def plot_pos_diff_xym1mat(xym1mat_cfile, scale=8, errMax=0.05):
    """
    Send in a match_c.txt file from xym1mat (3rd file) that has been trimmed
    of all bogus entries. Then plot of the differences in the positions
    from the 1st and 2nd starlists that went into the list matching process.
    This allows me to see any large scale systematics due to residual distortion
    or alignment errors.

    This assumes that the input starlists were in the XYMEEE format. We will
    use the frame #2 astrometric errors to trim out bad points.
    """
    t = atpy.Table(xym1mat_cfile, type='ascii')

    if errMax != None:
        x2e = t.col17
        y2e = t.col18

        err = np.hypot(x2e, y2e)

        origLen = len(t)
        t = t.where(err < errMax)
        newLen = len(t)

        print 'Trimmed %d of %d sources with errors > %.2f pix' % \
            (origLen-newLen, origLen, errMax)
        
    x1 = t.col8
    y1 = t.col9
    x2 = t.col11
    y2 = t.col12
    dx = x1 - x2
    dy = y1 - y2


    py.clf()
    q = py.quiver(x1, y1, dx, dy, scale=scale)
    py.quiverkey(q, 0.9, 0.95, 0.1, '0.1 pix', coordinates='axes', color='red')
    py.xlabel('X (pix)')
    py.ylabel('Y (pix)')
    py.savefig('pos_diff_vecs_' + xym1mat_cfile.replace('txt', 'png'))

  def sliceuv(self,ind,time=0,plot=False,**opts):
    savename=False
    if savename in opts.keys(): savename=opts['savename']

    if ind=='bar':
      isK,uname,vname,ind = True, 'ubar','vbar', 9999
    elif ind in ('s','surf','surface'):
      isK,uname,vname,ind = True,  'u','v', -1
    elif ind>=0:
      isK,uname,vname,ind = True,  'u','v', ind
    elif ind <0:
      isK=False
      isK,uname,vname,ind = False, 'u','v', ind


    if isK:
      xu,yu,zu,u=self.slicek(uname,ind,time)
      xv,yv,zv,v=self.slicek(vname,ind,time)
    else:
      xu,yu,zu,u=self.slicez(uname,ind,time)
      xv,yv,zv,v=self.slicez(vname,ind,time)

    z=(zu[1:,:]+zu[:-1,:])/2.
    u=( u[1:,:]+ u[:-1,:])/2.
    v=( v[:,1:]+ v[:,:-1])/2.

    x=self.grid.lonp
    y=self.grid.latp

    # rotate uv:
    ang=rt.rho2uvp(self.grid.angle,'p')
    u,v=calc.rot2d(u,v,-ang)


    if plot:
      p=self.grid.plot(bathy=None)
      xm, ym = p(x,y)
      mm=0*m
      mm[::3,::3]=1
      mm=mm*m==1
      s=np.sqrt(u**2+v**2)
      q=pylab.quiver(xm[mm],ym[mm],u[mm],v[mm],s[mm])
      pylab.colorbar(shrink=.7)
      pylab.axis([p.xmin, p.xmax, p.ymin, p.ymax])
      qk = pylab.quiverkey(q, .05, .01, 0.2, '0.2 m/s',labelpos='E',
                                              coordinates='axes')
      if savename:
        pylab.savefig(savename,dpi=pylab.gcf().dpi)
        pylab.close(pylab.gcf())

    return x,y,z,u,v

def plot_residuals_year_pos(year, filt, pos, refresh=False):
    year = str(year) 
    
    dir_xym = year + '_' + filt + '_' + pos + '/01.XYM/'

    # Try to read from a FITS table of compiled data if possible.
    if refresh == True:
        make_residuals_table_year_pos(year, filt, pos)

    d = read_residuals_table_year_pos(year, filt, pos)

    print d.xstd.min(), d.xstd.max()
    
    # Look at the variance in the raw X pixels.
    py.clf()
    q = py.quiver(d.xmean, d.ymean, d.xstd, d.ystd, scale=0.7)
    py.quiverkey(q, -0.1, -0.12, 0.02, '0.02 pixels', color='red')
    py.xlim(0, 1050)
    py.xlabel("X (pixels)")
    py.ylabel("Y (pixels)")
    py.title(year + ',' + filt + ',' + pos)
    
    py.savefig('plots/resid_errors_{0}_{1}_{2}.png'.format(year, filt, pos))

    def sliceuv(self, ind, time=0, plot=False, **opts):
        savename = False
        if savename in opts.keys():
            savename = opts["savename"]

        if ind == "bar":
            isK, uname, vname, ind = True, "ubar", "vbar", 9999
        elif ind in ("s", "surf", "surface"):
            isK, uname, vname, ind = True, "u", "v", -1
        elif ind >= 0:
            isK, uname, vname, ind = True, "u", "v", ind
        elif ind < 0:
            isK = False
            isK, uname, vname, ind = False, "u", "v", ind

        if isK:
            xu, yu, zu, u = self.slicek(uname, ind, time)
            xv, yv, zv, v = self.slicek(vname, ind, time)
        else:
            xu, yu, zu, u = self.slicez(uname, ind, time)
            xv, yv, zv, v = self.slicez(vname, ind, time)

        z = (zu[1:, :] + zu[:-1, :]) / 2.0
        u = (u[1:, :] + u[:-1, :]) / 2.0
        v = (v[:, 1:] + v[:, :-1]) / 2.0

        x = self.grid.lonp
        y = self.grid.latp

        # rotate uv:
        ang = rt.rho2uvp(self.grid.angle, "p")
        u, v = calc.rot2d(u, v, -ang)

        if plot:
            p = self.grid.plot(bathy=None)
            xm, ym = p(x, y)
            mm = 0 * m
            mm[::3, ::3] = 1
            mm = mm * m == 1
            s = np.sqrt(u ** 2 + v ** 2)
            q = pylab.quiver(xm[mm], ym[mm], u[mm], v[mm], s[mm])
            pylab.colorbar(shrink=0.7)
            pylab.axis([p.xmin, p.xmax, p.ymin, p.ymax])
            qk = pylab.quiverkey(q, 0.05, 0.01, 0.2, "0.2 m/s", labelpos="E", coordinates="axes")
            if savename:
                pylab.savefig(savename, dpi=pylab.gcf().dpi)
                pylab.close(pylab.gcf())

        return x, y, z, u, v

def display_2nd_moments(data=None):
    # remove the mean for all moments
    datamean = data.mean(axis = 0)
    data = subMeanAll(data)
    pl.figure(figsize=(11,5.5))
    pl.subplot(1,2,1)
    phi22 = 0.5*np.arctan2(data[:,3].imag,data[:,3].real)
    x = data[:,0].real
    y = data[:,1].real
    #phi22[x<0] = phi22+np.deg2rad(180)
    u = np.abs(data[:,3])*np.cos(phi22)
    v = np.abs(data[:,3])*np.sin(phi22)
    qvr = pl.quiver(x,y,u,v,width = 0.004, color='r',pivot='middle',headwidth=0.,headlength=0.,headaxislength=0.,scale_units='width')
    #qk = pl.quiverkey(qvr, -150,-240,np.max(np.sqrt(u**2+v**2)),str(round(np.max(np.sqrt(u**2+v**2)),3))+' pix^2',coordinates='data',color='blue')
    qk = pl.quiverkey(qvr, -150,-240,0.3,str(0.3)+' pix^2',coordinates='data',color='blue')
    pl.plot(x,y,'b,')
    pl.xlim(-250,250)
    pl.ylim(-250,250)
    pl.grid(color='g')
    pl.xlabel('Camera WEST [mm]')
    pl.ylabel('Camera NORTH [mm]')
    pl.title('mean |M22|= '+str(round(abs(datamean[:,3]),5))+' pix^2')
    pl.subplot(1,2,2)
    m20sqr = np.sqrt(data[:,2].real)
    x = data[:,0].real
    y = data[:,1].real
    m20sqr_med = np.median(m20sqr) 
    m20sqr_diff = m20sqr - m20sqr_med
    m20sqr_diff_absmed = np.median(np.abs(m20sqr_diff))
    plotScale = 1./m20sqr_diff_absmed*100
    pos = m20sqr_diff >=0
    neg = m20sqr_diff < 0
    pl.scatter(x[pos],y[pos],s=m20sqr_diff[pos]*plotScale,c='r',alpha=0.5)
    pl.scatter(x[neg],y[neg],s=-m20sqr_diff[neg]*plotScale,c='b',alpha=0.5)
    pl.scatter(-230,-210,s=0.01*plotScale,c='b',alpha=0.5)
    pl.text(-200,-215,'-'+str(0.01)+' pix')
    pl.scatter(-230,-230,s=0.01*plotScale,c='r',alpha=0.5)
    pl.text(-200,-235,str(0.01)+' pix')
    pl.plot(x,y,'y,')
    pl.grid(color='g')
    pl.xlim(-250,250)
    pl.ylim(-250,250)
    pl.xlabel('Camera WEST [mm]')
    pl.ylabel('Camera NORTH [mm]')
    pl.title('median '+r'$\sqrt{M20}$: '+str(round(scale*m20sqr_med,3))+' [arcsec]')
    return '---done!--'

def add_arrows(fieldx,fieldy):
	norm=abs(np.concatenate((fieldx,fieldy)).max())/1. #1d norm
	freq=4
	M = py.sqrt(pow(fieldx, 2) + pow(fieldy, 2))
	X,Y = py.meshgrid( py.arange(0,nx,1),py.arange(0,ny,1) )

        if(1==1):
            Q = py.quiver( 
		X[::freq, ::freq], Y[::freq, ::freq], 
		fieldx[::freq, ::freq], fieldy[::freq, ::freq], 
                #		M[::freq, ::freq],cmap=py.cm.jet,
                #		,color='yellow'
		#M[::freq, ::freq],cmap=py.cm.gray, alpha=0.2,
#		facecolors=mycmap(M[::freq, ::freq]),
                #		width=0.4,linewidths=(.5), edgecolors=('k'), 
#                color=colors,
		width=(0.2*py.sqrt(freq)),linewidths=(.2),
#		pivot='mid', 
#		pivot='tail', 
		units='x',
		headaxislength=5 ,
		scale=norm/freq)
            colors=mycmap_rgba(M[::freq, ::freq])
            Q.set_facecolor(colors)
            Q.set_edgecolor('k')
            Q.set_edgecolor(colors)
        else:
            Q = py.quiver( 
		X[::freq, ::freq], Y[::freq, ::freq], 
		fieldx[::freq, ::freq], fieldy[::freq, ::freq], 
                #		M[::freq, ::freq],cmap=py.cm.gray, alpha=1,
		width=0.2,
		units='x',
		headaxislength=5 ,
		scale=norm/2)

        legend='{:.2g}km/s'.format(norm)
        qk = py.quiverkey(Q, 0.9, 1.05, norm, legend,
                          labelpos='E',
                          fontproperties={'weight': 'bold'})

    def PlotFlightContour(self):

        # Creates list for UAV directions from 0 to 359
        uav_direction_0_359 = numpy.arange(0, 360, 0.1)

        # Creates empty list to store distance travelled by UAV
        displacement = []

        # Creates empty lists to store x and y coordinates
        x_coords = []
        y_coords = []

        # For-loop to run through the possible UAV directions
        for x in uav_direction_0_359:
            self.uav_U = self.uav_speed * numpy.cos(x * (2 * numpy.pi) / 360)
            self.uav_V = self.uav_speed * numpy.sin(x * (2 * numpy.pi) / 360)
            y = self.BaseFly()
            # Appends results of BaseFly into relevant lists
            displacement.append(y[0])
            x_coords.append(y[0] * numpy.cos(y[1] * (2 * numpy.pi) / 360))
            y_coords.append(y[0] * numpy.sin(y[1] * (2 * numpy.pi) / 360))

        # Generates 2D wind vector field
        X, Y = numpy.meshgrid(numpy.arange(min(x_coords) - 80, max(x_coords) + 80, (max(x_coords) - min(x_coords)) / 5),
                              numpy.arange(min(y_coords) - 80, max(y_coords) + 80, (max(y_coords) - min(y_coords)) / 5))
        Q = pylab.quiver(X, Y, self.wind_U, self.wind_V, label='Wind Direction')
        QK = pylab.quiverkey(Q, 0.1, -0.08, 10, 'Wind Direction', labelpos='W',
                             fontproperties={'weight': 'bold'})

        pylab.plot(x_coords, y_coords, self.marker3, label='Displacement Contour')
        pylab.plot(0, 0, self.marker1, label="Start")
        pylab.title('UAV Displacement in Y against UAV Displacement in X\n(Varying UAV Direction from 0 to 359)')
        pylab.xlabel('UAV Displacement in X')
        pylab.ylabel('UAV Displacement in Y')

        pylab.legend()

        pylab.show()

def makeVectorField(f, xmin, xmax, ymin, ymax, step):
    X, Y = meshgrid(arange(xmin, xmax, step), arange(ymin, ymax, step))
    U, V = zip(*map(lambda xx: f(*xx), zip(X, Y)))
    Q = quiver(X, Y, U, V, units='width')
    quiverkey(Q, 0, 0, 4, '', coordinates='figure', labelpos='W')

def plot_quiver_one_pass(reread=False):
    """
    Plot a quiver vector plot between F814W in 2005 and F160W in 2010 and 2013.
    This is just to check that we don't hae any gross flows between the two cameras.
    """
    if reread:
        t2005 = starlists.read_matchup('MATCHUP.XYMEEE.F814W.2005.ref5')
        t2010 = starlists.read_matchup('MATCHUP.XYMEEE.F160W.2010.ref5')
        t2013 = starlists.read_matchup('MATCHUP.XYMEEE.F160W.2013.ref5')

        t2005.write('MATCHUP.XYMEEE.F814W.2005.ref5.fits')
        t2010.write('MATCHUP.XYMEEE.F160W.2010.ref5.fits')
        t2013.write('MATCHUP.XYMEEE.F160W.2013.ref5.fits')
    else:
        t2005 = Table.read('MATCHUP.XYMEEE.F814W.2005.ref5.fits')
        t2010 = Table.read('MATCHUP.XYMEEE.F160W.2010.ref5.fits')
        t2013 = Table.read('MATCHUP.XYMEEE.F160W.2013.ref5.fits')

    good = np.where((t2005['m'] < -8) & (t2010['m'] < -8) & (t2013['m'] < -8) &
                    (t2005['xe'] < 0.05) & (t2010['xe'] < 0.05) & (t2013['xe'] < 0.05) & 
                    (t2005['ye'] < 0.05) & (t2010['ye'] < 0.05) & (t2013['ye'] < 0.05) & 
                    (t2005['me'] < 0.05) & (t2010['me'] < 0.05) & (t2013['me'] < 0.05))[0]

    g2005 = t2005[good]
    g2010 = t2010[good]
    g2013 = t2013[good]

    dx_05_10 = (g2010['x'] - g2005['x']) * ast.scale['WFC'] * 1e3
    dy_05_10 = (g2010['y'] - g2005['y']) * ast.scale['WFC'] * 1e3

    dx_05_13 = (g2013['x'] - g2005['x']) * ast.scale['WFC'] * 1e3
    dy_05_13 = (g2013['y'] - g2005['y']) * ast.scale['WFC'] * 1e3

    dx_10_13 = (g2013['x'] - g2010['x']) * ast.scale['WFC'] * 1e3
    dy_10_13 = (g2013['y'] - g2010['y']) * ast.scale['WFC'] * 1e3

    small = np.where((np.abs(dx_05_10) < 20) & (np.abs(dy_05_10) < 20) & 
                     (np.abs(dx_05_13) < 20) & (np.abs(dy_05_13) < 20) & 
                     (np.abs(dx_10_13) < 20) & (np.abs(dy_10_13) < 20))[0]

    print(len(g2005), len(small), len(dx_05_10))
    g2005 = g2005[small]
    g2010 = g2010[small]
    g2013 = g2013[small]
    dx_05_10 = dx_05_10[small]
    dx_05_13 = dx_05_13[small]
    dx_10_13 = dx_10_13[small]
    dy_05_10 = dy_05_10[small]
    dy_05_13 = dy_05_13[small]
    dy_10_13 = dy_10_13[small]
    print(len(g2005), len(small), len(dx_05_10))
    
    py.clf()
    q = py.quiver(g2005['x'], g2005['y'], dx_05_10, dy_05_10, scale=2e2)
    py.quiverkey(q, 0.95, 0.95, 5, '5 mas', color='red', labelcolor='red')
    py.title('2010 - 2005')
    py.savefig('vec_diff_ref5_05_10.png')

    py.clf()
    q = py.quiver(g2005['x'], g2005['y'], dx_05_13, dy_05_13, scale=2e2)
    py.quiverkey(q, 0.95, 0.95, 5, '5 mas', color='red', labelcolor='red')
    py.title('2013 - 2005')
    py.savefig('vec_diff_ref5_05_13.png')

    py.clf()
    q = py.quiver(g2005['x'], g2005['y'], dx_10_13, dy_10_13, scale=1e2)
    py.quiverkey(q, 0.95, 0.95, 5, '5 mas', color='red', labelcolor='red')
    py.title('2013 - 2010')
    py.savefig('vec_diff_ref5_10_13.png')

    py.clf()
    py.plot(dx_05_10, dy_05_10, 'k.', ms=2)
    lim = 10
    py.axis([-lim, lim, -lim, lim])
    py.xlabel('X Proper Motion (mas)')
    py.ylabel('Y Proper Motion (mas)')
    py.title('2010 - 2005')
    py.savefig('pm_diff_ref5_05_10.png')

    py.clf()
    py.plot(dx_05_13, dy_05_13, 'k.', ms=2)
    py.axis([-lim, lim, -lim, lim])
    py.xlabel('X Proper Motion (mas)')
    py.ylabel('Y Proper Motion (mas)')
    py.title('2013 - 2005')
    py.savefig('pm_diff_ref5_05_13.png')

    py.clf()
    py.plot(dx_10_13, dy_10_13, 'k.', ms=2)
    py.axis([-lim, lim, -lim, lim])
    py.xlabel('X Proper Motion (mas)')
    py.ylabel('Y Proper Motion (mas)')
    py.title('2013 - 2010')
    py.savefig('pm_diff_ref5_10_13.png')

    print('2010 - 2005')
    print('   dx = {dx:6.2f} +/- {dxe:6.2f} mas'.format(dx=dx_05_10.mean(), dxe=dx_05_10.std()))
    print('   dy = {dy:6.2f} +/- {dye:6.2f} mas'.format(dy=dy_05_10.mean(), dye=dy_05_10.std()))

    print('2013 - 2005')
#.........这里部分代码省略.........

def plot_vpd_across_field(nside=4, interact=False):
    """
    Plot the VPD at different field positions so we can see if there are
    systematic discrepancies due to residual distortions.
    """
    # Read in matched and aligned star lists from the *.ref5 analysis.
    # Recall these data sets are in the F814W reference frame with a 50 mas plate scale.
    t2005 = starlists.read_matchup(workDir + '02.MAT/MATCHUP.XYMEEE.F814W.ref5')
    t2010 = starlists.read_matchup(workDir + '02.MAT/MATCHUP.XYMEEE.F125W.ref5')

    scale = 50.0 # mas per pixel
    

    # Trim down to only those stars that are detected in both epochs.
    # Also make cuts on astrometric/photometric errors, etc.
    # We only want the well measured stars in this analysis.
    perrLim814 = 1.0 / scale
    perrLim125 = 4.0 / scale
    merrLim814 = 0.05
    merrLim125 = 0.1
    
    cond = ((t2005.m != 0) & (t2010.m != 0) &
            (t2005.xe < perrLim814) & (t2005.ye < perrLim814) &
            (t2010.xe < perrLim125) & (t2010.ye < perrLim125) &
            (t2005.me < merrLim814) & (t2010.me < merrLim125))

    t2005 = t2005.where(cond)
    t2010 = t2010.where(cond)

    # Calculate proper motions
    dt = years['2010_F125W'] - years['2005_F814W']
    dx = t2010.x - t2005.x
    dy = t2010.y - t2005.y
    pmx = dx * scale / dt
    pmy = dy * scale / dt
    pm = np.hypot(pmx, pmy)

    t2005.add_column('pmx', pmx)
    t2005.add_column('pmy', pmy)
    t2005.add_column('pm', pm)


    # Divide up the region into N x N boxes and plot up the VPD for each.
    xlo = math.floor(t2005.x.min())
    xhi = math.ceil(t2005.x.max())
    ylo = math.floor(t2005.y.min())
    yhi = math.ceil(t2005.y.max())
    xboxsize = round((xhi - xlo) / nside)
    yboxsize = round((yhi - ylo) / nside)

    # Setup colors
    jet = py.get_cmap('jet')
    cNorm = colors.Normalize(vmin=0, vmax=nside**2)
    colorMap = py.cm.ScalarMappable(norm=cNorm, cmap=jet)

    # Save the average proper motions in each box
    pmx = np.zeros((nside, nside), dtype=float)
    pmy = np.zeros((nside, nside), dtype=float)
    pmxe = np.zeros((nside, nside), dtype=float)
    pmye = np.zeros((nside, nside), dtype=float)
    xcen = np.zeros((nside, nside), dtype=float)
    ycen = np.zeros((nside, nside), dtype=float)

    pmCut = 1.0

    # Calculate the global mean proper motion
    # Start by trimming down to a 1 mas/yr radius
    idx2 = np.where(pm < pmCut)[0]
    pmx_all = np.median( t2005.pmx[idx2] )
    pmy_all = np.median( t2005.pmy[idx2] )
    
    out = 'All X:{0:5.0f}-{1:5.0f}  Y:{2:5.0f}-{3:5.0f}  '
    out += 'PMX:{4:5.2f} +/- {5:5.2f} PMY:{6:5.2f} +/- {7:5.2f}  '
    out += 'N:{8:5d}'
    print((out.format(xlo, xhi, ylo, yhi, pmx_all, 0.0, pmy_all, 0.0, len(idx2))))

    # Make a global proper motion diagram of star with a proper motion within
    # 1 mas/yr. This is mainly to see systematic flows due to residual distortion.
    pmTot = np.hypot(t2005.pmx, t2005.pmy)
    clust = np.where(pmTot < pmCut)[0]
    py.clf()
    q = py.quiver(t2005.x[clust], t2005.y[clust], t2005.pmx[clust], t2005.pmy[clust],
                  scale=18)
    py.quiverkey(q, 0.5, 0.98, 1, '1 mas/yr', color='red', labelcolor='red')
    py.xlabel('X Position (pixels)')
    py.ylabel('Y Position (pixels)')
    py.xlim(xlo, xhi)
    py.ylim(ylo, yhi)
    out = '{0}/plots/vec_proper_motion_all.png'
    py.savefig(out.format(workDir))
    
    py.clf()
    for xx in range(nside):
        for yy in range(nside):
            xlo_box = xlo + xx * xboxsize
            ylo_box = ylo + yy * yboxsize
            xhi_box = xlo + (1+xx) * xboxsize
            yhi_box = ylo + (1+yy) * yboxsize

            idx = np.where((t2005.x > xlo_box) & (t2005.x <= xhi_box) &
#.........这里部分代码省略.........

def velocity_image(sim, width="10 kpc", vector_color='black', edgecolor='black', quiverkey_bg_color=None,
                   vector_resolution=40, scale=None, mode='quiver', key_x=0.3, key_y=0.9,
                   key_color='white', key_length="100 km s**-1", quiverkey=True, density=1.0, **kwargs):
    """

    Make an SPH image of the given simulation with velocity vectors overlaid on top.

    For a description of additional keyword arguments see :func:`~pynbody.plot.sph.image`,
    or see the `tutorial <http://pynbody.github.io/pynbody/tutorials/pictures.html#velocity-vectors>`_.

    **Keyword arguments:**

    *vector_color* (black): The color for the velocity vectors

    *edgecolor* (black): edge color used for the lines - using a color
     other than black for the *vector_color* and a black *edgecolor*
     can result in poor readability in pdfs

    *vector_resolution* (40): How many vectors in each dimension (default is 40x40)

    *quiverkey_bg_color* (none): The color for the legend (scale) background

    *scale* (None): The length of a vector that would result in a displayed length of the
    figure width/height.

    *mode* ('quiver'): make a 'quiver' or 'stream' plot

    *key_x* (0.3): Display x (width) position for the vector key (quiver mode only)

    *key_y* (0.9): Display y (height) position for the vector key (quiver mode only)

    *key_color* (white): Color for the vector key (quiver mode only)

    *key_length* (100 km/s): Velocity to use for the vector key (quiver mode only)

    *density* (1.0): Density of stream lines (stream mode only)

    *quiverkey* (True): Whether or not to inset the key

    """

    subplot = kwargs.get('subplot', False)
    av_z = kwargs.get('av_z',None)
    if subplot:
        p = subplot
    else:
        import matplotlib