def draw_force_field(x_list, y_list, width, height):
    from pprint import pprint
    '''docstring for draw_force_field''' 
    log.info("draw_force_field is working..." )
    X = []
    Y = []
    for row in range(0,height):
        one_row_x = []
        one_row_y = []
        for col in range(0, width):
            i = row * width + col
            one_row_x.append(x_list[i])
            one_row_y.append(y_list[i])
        X.append(one_row_x)
        Y.append(one_row_y)
        #X.insert(0, one_row_x)
        #Y.insert(0, one_row_y)

    Q = plt.quiver( X, Y)
    plt.quiverkey(Q, 0.5, 0.92, 2, r'')
    #l,r,b,t = plt.axis()
    #dx, dy = r-l, t-b
    ##axis([l-0.05*dx, r+0.05*dx, b-0.05*dy, t+0.05*dy])
    #plt.xticks(range(0,width +3))
    #plt.yticks(range(0,height + 3))
    plt.title('Attraction Force field')
    plt.savefig("pixel.jpg")
    plt.show()
    return 0

def plot_meancontquiv():
    data = calcwake(t1=0.4)
    y_R = data["y/R"]
    z_R = data["z/R"]
    u = data["meanu"]
    v = data["meanv"]
    w = data["meanw"]
    plt.figure(figsize=(7, 9))
    # Add contours of mean velocity
    cs = plt.contourf(y_R, z_R, u/U, 20, cmap=plt.cm.coolwarm)
    cb = plt.colorbar(cs, shrink=1, extend="both",
                      orientation="horizontal", pad=0.1)
                      #ticks=np.round(np.linspace(0.44, 1.12, 10), decimals=2))
    cb.set_label(r"$U/U_{\infty}$")
    # Make quiver plot of v and w velocities
    Q = plt.quiver(y_R, z_R, v/U, w/U, angles="xy", width=0.0022,
                   edgecolor="none", scale=3.0)
    plt.xlabel(r"$y/R$")
    plt.ylabel(r"$z/R$")
    plt.quiverkey(Q, 0.8, 0.21, 0.1, r"$0.1 U_\infty$",
               labelpos="E",
               coordinates="figure",
               fontproperties={"size": "small"})
    ax = plt.axes()
    ax.set_aspect(1)
    # Plot circle to represent turbine frontal area
    circ = plt.Circle((0, 0), radius=1, facecolor="none", edgecolor="gray",
                      linewidth=3.0)
    ax.add_patch(circ)
    plt.tight_layout()

def PLOT_E(v1,v2,v3,v4,scale,title="",flag_out=0,name_out="output.png",flag_show=0):
	if not len(v1) == len (v2):
		print("v1 and v2 have different length")
		return
	if not len(v1) == len (v3):
		print("v1 and v3 have different length")
		return
	if not len(v1) == len (v4):
		print("v1 and v4 have different length")
		return
	E1 = []
	E2 = []
	for i in range(len(v1)):
		E=math.hypot(v3[i],v4[i])
		T=math.atan2(v4[i],v3[i])
		E1.append( E*math.cos(0.5*T))
		E2.append( E*math.sin(0.5*T))

	plt.title(title)
	A=plt.quiver(v1,v2,E1,E2,scale=(float)(scale),width=0.005,headwidth=0,pivot='middle')
	plt.quiverkey(A,0.95,0.95,0.1,"0.1")
	plt.axis('scaled')
	if flag_out == 1:
		plt.savefig(name_out)
	if flag_show == 1:
		plt.show()
	plt.close()

def main():

    plot_utils.apply_plot_params(width_cm=20, height_cm=20, font_size=10)

    high_hles_years = [1993, 1995, 1998]
    low_hles_years = [1997, 2001, 2006]
    data_path = "/BIG1/skynet1_rech1/diro/sample_obsdata/eraint/eraint_uvslp_years_198111_201102_NDJmean_ts.nc"





    with xr.open_dataset(data_path) as ds:
        print(ds)


        u = get_composit_for_name(ds, "u10", high_years_list=high_hles_years, low_years_list=low_hles_years)
        v = get_composit_for_name(ds, "v10", high_years_list=high_hles_years, low_years_list=low_hles_years)
        msl = get_composit_for_name(ds, "msl", high_years_list=high_hles_years, low_years_list=low_hles_years)

        lons = ds["longitude"].values
        lats = ds["latitude"].values

        print(lats)
        print(msl.shape)
        print(lons.shape, lats.shape)


        lons2d, lats2d = np.meshgrid(lons, lats)


    fig = plt.figure()

    map = Basemap(llcrnrlon=-130, llcrnrlat=22, urcrnrlon=-28,
                  urcrnrlat=65, projection='lcc', lat_1=33, lat_2=45,
                  lon_0=-95, resolution='i', area_thresh=10000)


    clevs = np.arange(-11.5, 12, 1)
    cmap = cm.get_cmap("bwr", len(clevs) - 1)
    bn = BoundaryNorm(clevs, len(clevs) - 1)

    x, y = map(lons2d, lats2d)
    im = map.contourf(x, y, msl / 100, levels=clevs, norm=bn, cmap=cmap) # convert to mb (i.e hpa)
    map.colorbar(im)


    stride = 2
    ux, vy = map.rotate_vector(u, v, lons2d, lats2d)
    qk = map.quiver(x[::stride, ::stride], y[::stride, ::stride], ux[::stride, ::stride], vy[::stride, ::stride],
               scale=10, width=0.01, units="inches")
    plt.quiverkey(qk, 0.5, -0.1, 2, "2 m/s", coordinates="axes")


    map.drawcoastlines(linewidth=0.5)
    map.drawcountries()
    map.drawstates()
    #plt.show()

    fig.savefig("hles_wind_compoosits.png", bbox_inches="tight", dpi=300)

    def whiskerplot(cat,col,fig):

        if col == 'psf':
            key = r'e_{PSF}'
        if col == 'e':
            key = r'e'
        if col == 'dpsf':
            key = r'\Delta e_{PSF}'

        scale=0.02

        y,x,mw,e1,e2,e=field.field.whisker_calc(cat,col=col)
        pos0=0.5*np.arctan2(e2/mw,e1/mw)
        e/=mw
        for i in range(len(x)):
            y[i,:,:],x[i,:,:]=field.field_methods.ccd_to_field(i,y[i,:,:]-2048,x[i,:,:]-1024)

        print 'y,x',y[i,:,:],x[i,:,:]

        plt.figure(fig)
        print np.shape(x),np.shape(y),np.shape(np.sin(pos0)*e),np.shape(np.cos(pos0)*e)
        Q = plt.quiver(np.ravel(y),np.ravel(x),np.ravel(np.sin(pos0)*e),np.ravel(np.cos(pos0)*e),units='width',pivot='middle',headwidth=0,width=.0005)
        plt.quiverkey(Q,0.2,0.2,scale,str(scale)+' '+key,labelpos='E',coordinates='figure',fontproperties={'weight': 'bold'})
        plt.savefig('plots/y1/whisker_'+col+'.pdf', dpi=500, bbox_inches='tight')
        plt.close(fig)

        return

def add_std_vector_to_12months_cycle_figure(plt, fig, CF, parameter):
	std_length	= parameter[0]
	scale		= parameter[1]
	qkeyx		= parameter[2]
	qkeyy		= parameter[3]
	plt.quiverkey(CF, qkeyx, qkeyy, std_length, str(std_length))
	return plt

def plot_meancontquiv(turbine="turbine2", save=False):
    mean_u = load_vel_map(turbine=turbine, component="u")
    mean_v = load_vel_map(turbine=turbine, component="v")
    mean_w = load_vel_map(turbine=turbine, component="w")
    y_R = np.round(np.asarray(mean_u.columns.values, dtype=float), decimals=4)
    z_R = np.asarray(mean_u.index.values, dtype=float)
    plt.figure(figsize=(7.5, 4.1))
    # Add contours of mean velocity
    cs = plt.contourf(y_R, z_R, mean_u / U, 20, cmap=plt.cm.coolwarm)
    cb = plt.colorbar(cs, orientation="vertical")
    cb.set_label(r"$U/U_{\infty}$")
    # Make quiver plot of v and w velocities
    Q = plt.quiver(y_R, z_R, mean_v / U, mean_w / U, angles="xy", width=0.0022, edgecolor="none", scale=3.0)
    plt.xlabel(r"$y/R$")
    plt.ylabel(r"$z/R$")
    plt.quiverkey(
        Q, 0.65, 0.045, 0.1, r"$0.1 U_\infty$", labelpos="E", coordinates="figure", fontproperties={"size": "small"}
    )
    ax = plt.axes()
    ax.set_aspect(1)
    # Plot circle to represent turbine frontal area
    circ = plt.Circle((0, 0), radius=1, facecolor="none", edgecolor="gray", linewidth=3.0)
    ax.add_patch(circ)
    plt.tight_layout()
    if save:
        plt.savefig("figures/" + turbine + "-meancontquiv.pdf")

def plot_grid(EIC_grid,sat_track,sat_krig,title):
    '''
    This function plots a scatter map of the EICS grid and its horizontal components, and the krigged value for the
    satellite.

    :param EIC_grid: The EICS grid
    :param satPos: The position of the satellite
    :param ptz_u: The krigged u component of the satellite
    :param ptz_v: The krigged v component of the satllite
    :param title: Timestamp of the satellite
    :return: The figure
    '''

    # Define the size of the figure in inches
    plt.figure(figsize=(18,18))

    # The horizontal components of the Ionospheric current from the EICS grid
    u = EIC_grid[:,2]
    v = EIC_grid[:,3]

    '''
    The m variable defines the basemap of area that is going to be displayed.
    1) width and height is the area in pixels of the area to be displayed.
    2) resolution is the resolution of the boundary dataset being used 'c' for crude and 'l' for low
    3) projection is type of projection of the basemape, in this case it is a Lambert Azimuthal Equal Area projection
    4) lat_ts is the latitude of true scale,
    5) lat_0 and lon_0 is the latitude and longitude of the central point of the basemap
    '''
    m = Basemap(width=8000000, height=8000000, resolution='l', projection='lcc',\
             lat_0=60,lon_0=-100.)

    m.drawcoastlines() #draw the coastlines on the basemap

    # draw parallels and meridians and label them
    m.drawparallels(np.arange(-80.,81.,20.),labels=[1,0,0,0],fontsize=10)
    m.drawmeridians(np.arange(-180.,181.,20.),labels=[0,0,0,1],fontsize=10)

    # Project the inputted grid into x,y values defined of the projected basemap parameters
    x,y =m(EIC_grid[:,1],EIC_grid[:,0])
    satx,saty = m(sat_track[::10,7],sat_track[::10,6])
    satkrigx,satkrigy = m(sat_krig[1],sat_krig[0])


    '''
    Plot the inputted grid as a quiver plot on the basemap,
    1) x,y are the projected latitude and longitude coordinates of the grid
    2) u and v are the horizontal components of the current
    3) the EICS grid values are plotted in blue color where as the satellite krigged values are in red
    '''
    eic = m.quiver(x,y,u,v,width = 0.004, scale=10000,color='#0000FF')
    satkrig = m.quiver(satkrigx,satkrigy,sat_krig[2],sat_krig[3],color='#FF0000',width=0.004, scale = 10000)
    satpos = m.scatter(satx,saty,s=100,marker='.',c='#009933',edgecolors='none',label='Satellite Track')
    sat_halo = m.scatter(satkrigx,satkrigy,s=400,facecolors='none',edgecolors='#66FF66',linewidth='5')

    plt.title(title)
    plt.legend([satpos],['Satellite Track'],loc='upper right',scatterpoints=1)
    plt.quiverkey(satkrig,0.891,0.948,520,u'\u00B1' +'520 mA/m',labelpos='E')

    plt.savefig('EICS_20110311_002400.png',bbox_inches='tight',pad_inches=0.2)

def water_quiver1(current_data):
    u = water_u1(current_data)
    v = water_v1(current_data)

    plt.hold(True)
    Q = plt.quiver(current_data.x[::2, ::2], current_data.y[::2, ::2], u[::2, ::2], v[::2, ::2])
    max_speed = np.max(np.sqrt(u ** 2 + v ** 2))
    label = r"%s m/s" % str(np.ceil(0.5 * max_speed))
    plt.quiverkey(Q, 0.15, 0.95, 0.5 * max_speed, label, labelpos="W")
    plt.hold(False)

 def animate(n): # the function of the animation
     ax.cla()
     plt.title('Drifter: {0} {1}'.format(drifter_ID,point['time'][n].strftime("%F %H:%M")))
     draw_basemap(ax, points)
     ax.plot(drifter_points['lon'],drifter_points['lat'],'bo-',markersize=6,label='Drifter')
     ax.annotate(an2,xy=(dr_points['lon'][-1],dr_points['lat'][-1]),xytext=(dr_points['lon'][-1]+0.01*track_days,
                 dr_points['lat'][-1]+0.01*track_days),fontsize=6,arrowprops=dict(arrowstyle="fancy"))
     ax.plot(model_points['lon'][:n+1],model_points['lat'][:n+1],'ro-',markersize=6,label=MODEL)
     #M = np.hypot(U[n], V[n])
     Q = ax.quiver(X,Y,U[n],V[n],color='black',pivot='tail',units='xy')
     plt.quiverkey(Q, 0.5, 0.92, 1, r'$1 \frac{m}{s}$', labelpos='E',fontproperties={'weight': 'bold','size':18})

def plotResPosArrow2D(ccdSet, iexp, matchVec, sourceVec, outputDir):
    import matplotlib.pyplot as plt
    _xm = []
    _ym = []
    _dxm = []
    _dym = []
    for m in matchVec:
        if (m.good == True and m.iexp == iexp):
            _xm.append(m.u)
            _ym.append(m.v)
            _dxm.append((m.xi_fit - m.xi)*3600)
            _dym.append((m.eta_fit - m.eta)*3600)
    _xs = []
    _ys = []
    _dxs = []
    _dys = []
    if (len(sourceVec) != 0):
        for s in sourceVec:
            if (s.good == True and s.iexp == iexp):
                _xs.append(s.u)
                _ys.append(s.v)
                _dxs.append((s.xi_fit - s.xi)*3600)
                _dys.append((s.eta_fit - s.eta)*3600)

    xm = numpy.array(_xm)
    ym = numpy.array(_ym)
    dxm = numpy.array(_dxm)
    dym = numpy.array(_dym)
    xs = numpy.array(_xs)
    ys = numpy.array(_ys)
    dxs = numpy.array(_dxs)
    dys = numpy.array(_dys)

    plt.clf()
    plt.rc("text", usetex=USETEX)
    plt.rc('xtick', labelsize=10)
    plt.rc('ytick', labelsize=10)

    plotCcd(ccdSet)
    q = plt.quiver(xm, ym, dxm, dym, units="inches", angles="xy", scale=1, color="green", label="external")
    if len(xm) != 0 and len(ym) != 0:
        xPos = round(xm.min() + (xm.max() - xm.min())*0.002, -2)
        yPos = round(ym.max() + (ym.max() - ym.min())*0.025, -2)
        plt.quiverkey(q, xPos, yPos, 0.1, "0.1 arcsec", coordinates="data", color="blue", labelcolor="blue",
                      labelpos='E', fontproperties={'size': 10})
    plt.quiver(xs, ys, dxs, dys, units="inches", angles="xy", scale=1, color="red", label="internal")

    plt.axes().set_aspect("equal")
    plt.legend(fontsize=8)
    plt.title("LSST: %d" % (iexp))
    plt.savefig(os.path.join(outputDir, "ResPosArrow2D_%d.png" % (iexp)), format="png")

  def plot_whisker(x,y,e1,e2,name='',label='',scale=.01,key='',chip=False):

    plt.figure()
    Q = plt.quiver(x,y,e1,e2,units='width',pivot='middle',headwidth=0,width=.0005)
    if chip:
      plt.quiverkey(Q,0.2,0.125,scale,str(scale)+' '+key,labelpos='E',coordinates='figure',fontproperties={'weight': 'bold'})
      plt.xlim((-250,4250))
      plt.ylim((-200,2100))
    else:
      plt.quiverkey(Q,0.2,0.2,scale,str(scale)+' '+key,labelpos='E',coordinates='figure',fontproperties={'weight': 'bold'})
    plt.savefig('plots/footprint/whisker_'+name+'_'+label+'.png', dpi=500,bbox_inches='tight')
    plt.close()

    return

def quick_static(gflist,datapath,run_name,run_num,c):
    '''
    Make quick quiver plot of static fields
    
    IN:
        gflist: Tod ecide which stations to plot
        datapath: Where are the data files
    '''
    import matplotlib.pyplot as plt
    from numpy import genfromtxt,where,zeros,meshgrid,linspace

    
    GF=genfromtxt(gflist,usecols=3)
    sta=genfromtxt(gflist,usecols=0,dtype='S')
    lon=genfromtxt(gflist,usecols=1,dtype='f')
    lat=genfromtxt(gflist,usecols=2,dtype='f')
    #Read coseismcis
    i=where(GF!=0)[0]
    lon=lon[i]
    lat=lat[i]
    n=zeros(len(i))
    e=zeros(len(i))
    u=zeros(len(i))
    #Get data
    if run_name!='' or run_num!='':
        run_name=run_name+'.'
        run_num=run_num+'.'
    for k in range(len(i)):
        neu=genfromtxt(datapath+run_name+run_num+sta[i[k]]+'.static.neu')
        n[k]=neu[0]#/((neu[0]**2+neu[1]**2)**0.5)
        e[k]=neu[1]#/((neu[0]**2+neu[1]**2)**0.5)
        u[k]=neu[2]#/(2*abs(neu[2]))

            
    #Plot
    plt.figure()
    xi = linspace(min(lon), max(lon), 500)
    yi = linspace(min(lat), max(lat), 500)
    X, Y = meshgrid(xi, yi)
    #c=Colormap('bwr')
    #plt.contourf(X,Y,Z,100)
    #plt.colorbar()
    Q=plt.quiver(lon,lat,e,n,width=0.001,color=c)
    plt.scatter(lon,lat,color='b')
    plt.grid()
    plt.title(datapath+run_name+run_num)
    plt.show()
    qscale_en=1
    plt.quiverkey(Q,X=0.1,Y=0.9,U=qscale_en,label=str(qscale_en)+'m')

def overdraw_vec_with_axis(plt, datax, datay, xgrid, ygrid, std_length=0.1, scale=1, intvl=3, qkeyx=0.8, qkeyy=0.8, lw = 2):
	import matplotlib.pyplot as plt
	if datax.shape!=datay.shape:
		raise Exception('datax and datay don\'t have same shape!')

	if datax.shape[0] != ygrid.size:
		raise Exception('data size x does not match ygrid size!')

	if datay.shape[1] != xgrid.size:
		raise Exception('data size y does not match xgrid size!')

	X, Y=np.meshgrid(xgrid, ygrid)
	Q=plt.quiver(X[::intvl, ::intvl], Y[::intvl, ::intvl], datax[::intvl, ::intvl], datay[::intvl, ::intvl], angles='xy', scale=scale, lw = lw)
	plt.quiverkey(Q,  qkeyx,  qkeyy,  std_length,  str(std_length))
	return plt

    def plotmap(self, domain = [0., 360., -90., 90.], res='c', stepp=2, scale=20):

        latitudes = self.windspeed.latitudes.data
        longitudes = self.windspeed.longitudes.data

        m = bm(projection='cyl',llcrnrlat=latitudes.min(),urcrnrlat=latitudes.max(),\
        llcrnrlon=longitudes.min(),urcrnrlon=longitudes.max(),\
        lat_ts=0, resolution=res)

        lons, lats = np.meshgrid(longitudes, latitudes)

        cmap = palettable.colorbrewer.sequential.Oranges_9.mpl_colormap

        f, ax = plt.subplots(figsize=(10,6))

        m.ax = ax

        x, y = m(lons, lats)

        im = m.pcolormesh(lons, lats, self.windspeed.data, cmap=cmap)

        cb = m.colorbar(im)
        cb.set_label('wind speed (m/s)', fontsize=14)

        Q = m.quiver(x[::stepp,::stepp], y[::stepp,::stepp], \
                     self.uanoms.data[::stepp,::stepp], self.vanoms.data[::stepp,::stepp], \
                     pivot='middle', scale=scale)

        l,b,w,h = ax.get_position().bounds

        qk = plt.quiverkey(Q, l+w-0.1, b-0.03, 5, "5 m/s", labelpos='E', fontproperties={'size':14}, coordinates='figure')

        m.drawcoastlines()

        return f

def my_plot(u,v,t,daystr,levels):
    #boston light swim
    ax= [-71.10, -70.10, 41.70, 42.70] # region to plot
    vel_arrow = 0.2 # velocity arrow scale
    subsample = 8  # subsampling of velocity vectors

    # find velocity points in bounding box
    ind = np.argwhere((lonc >= ax[0]) & (lonc <= ax[1]) & (latc >= ax[2]) & (latc <= ax[3]))

    np.random.shuffle(ind)
    Nvec = int(len(ind) / subsample)
    idv = ind[:Nvec]
    # tricontourf plot of water depth with vectors on top
    plt.figure(figsize=(20,10))
    plt.subplot(111,aspect=(1.0/np.cos(lat[:].mean()*np.pi/180.0)))
    #tricontourf(tri, t,levels=levels,shading='faceted',cmap=plt.cm.gist_earth)
    plt.tricontourf(tri, t,levels=levels,shading='faceted')
    plt.axis(ax)
    plt.gca().patch.set_facecolor('0.5')
    cbar=plt.colorbar()
    cbar.set_label('Forecast Surface Temperature (C)', rotation=-90)
    plt.tricontour(tri, t,levels=[0])
    Q = plt.quiver(lonc[idv],latc[idv],u[idv],v[idv],scale=10)
    maxstr='%3.1f m/s' % vel_arrow
    qk = plt.quiverkey(Q,0.92,0.08,vel_arrow,maxstr,labelpos='W')
    plt.title('NECOFS Surface Velocity, Layer %d, %s UTC' % (ilayer, daystr))
    plt.plot(lon_track,lat_track,'m-o')
    plt.plot(lon_buoy,lat_buoy,'y-o')

    def velocidadlimpio(self):
        #sino no llega al 2, asi funciona        
        v = np.arange(-3, 3,1)
        [x,y] = np.meshgrid(v,v)

        z=np.multiply(x,np.exp( -np.power(x,2) - np.power(y,2) ))

        #Matplotlib t invierte el orden de las matrices a diferenciade matlab
        [py,px] = np.gradient(z,1,1)

        print 'x '+str(x)
        print 'y '+str(y)
        print 'z '+str(z)
        print 'px '+str(px)
        print 'py '+str(py)
        
        #q = plt.quiver(X, Y, u, v, angles='xy', scale=40, color=['r'])
        #p = plt.quiverkey(q,1,16.5,50,"50 m/s",coordinates='data',color='r')        

        #primero los rangos, despues los valores q contiene
        q = plt.quiver(x,y, px, py)
        p = plt.quiverkey(q,1,16.5,50,"50 m/s",coordinates='data',color='r')
        plt.title('Velocidad')
        plt.show()
        print 'Fourth plot loaded...'  

def plot_meancontquiv(t1_fraction=0.5, cb_orientation="vertical", save=False):
    """Plot contours of normalized mean streamwise velocity and vectors of
    cross-stream and vertical mean velocity.
    """
    data = load_probe_data(t1_fraction=t1_fraction)
    y_R = data["y_R"]
    z_H = data["z_H"]
    mean_u = data["mean_u"]
    mean_v = data["mean_v"]
    mean_w = data["mean_w"]
    scale = 7.5/10.0
    plt.figure(figsize=(10*scale, 3*scale))
    # Add contours of mean velocity
    cs = plt.contourf(y_R, z_H, mean_u, 20, cmap=plt.cm.coolwarm)
    if cb_orientation == "horizontal":
        cb = plt.colorbar(cs, shrink=1, extend="both",
                          orientation="horizontal", pad=0.14)
    elif cb_orientation == "vertical":
        cb = plt.colorbar(cs, shrink=0.88, extend="both",
                          orientation="vertical", pad=0.02)
    cb.set_label(r"$U/U_{\infty}$")
    # Make quiver plot of v and w velocities
    Q = plt.quiver(y_R, z_H, mean_v, mean_w, width=0.0022,
                   edgecolor="none", scale=3)
    plt.xlabel(r"$y/R$")
    plt.ylabel(r"$z/H$")
    plt.ylim(-0.15, 0.85)
    plt.xlim(-1.68/R, 1.68/R)
    if cb_orientation == "horizontal":
        plt.quiverkey(Q, 0.65, 0.26, 0.1, r"$0.1 U_\infty$",
                      labelpos="E",
                      coordinates="figure",
                      fontproperties={"size": "small"})
    elif cb_orientation == "vertical":
        plt.quiverkey(Q, 0.65, 0.09, 0.1, r"$0.1 U_\infty$",
                      labelpos="E",
                      coordinates="figure",
                      fontproperties={"size": "small"})
    plot_turb_lines()
    ax = plt.axes()
    ax.set_aspect(H/R)
    plt.yticks([0, 0.13, 0.25, 0.38, 0.5, 0.63, 0.75])
    plt.grid(True)
    plt.tight_layout()
    if save:
        plt.savefig("figures/meancontquiv.pdf")
        plt.savefig("figures/meancontquiv.png", dpi=300)

def plot_meancontquiv(save=False, show=False, savetype=".pdf",
                      cb_orientation="vertical"):
    """Plot mean contours/quivers of velocity."""
    mean_u = load_vel_map("u")
    mean_v = load_vel_map("v")
    mean_w = load_vel_map("w")
    y_R = np.round(np.asarray(mean_u.columns.values, dtype=float), decimals=4)
    z_H = np.asarray(mean_u.index.values, dtype=float)
    plt.figure(figsize=(7.5, 4.8))
    # Add contours of mean velocity
    cs = plt.contourf(y_R, z_H, mean_u,
                      np.arange(0.15, 1.25, 0.05), cmap=plt.cm.coolwarm)
    if cb_orientation == "horizontal":
        cb = plt.colorbar(cs, shrink=1, extend="both",
                          orientation="horizontal", pad=0.14)
    elif cb_orientation == "vertical":
        cb = plt.colorbar(cs, shrink=1, extend="both",
                          orientation="vertical", pad=0.02)
    cb.set_label(r"$U/U_{\infty}$")
    plt.hold(True)
    # Make quiver plot of v and w velocities
    Q = plt.quiver(y_R, z_H, mean_v, mean_w, width=0.0022,
                   edgecolor="none", scale=3.0)
    plt.xlabel(r"$y/R$")
    plt.ylabel(r"$z/H$")
    # plt.ylim(-0.2, 0.78)
    # plt.xlim(-3.2, 3.2)
    if cb_orientation == "horizontal":
        plt.quiverkey(Q, 0.65, 0.26, 0.1, r"$0.1 U_\infty$",
                      labelpos="E",
                      coordinates="figure",
                      fontproperties={"size": "small"})
    elif cb_orientation == "vertical":
        plt.quiverkey(Q, 0.65, 0.055, 0.1, r"$0.1 U_\infty$",
                      labelpos="E",
                      coordinates="figure",
                      fontproperties={"size": "small"})
    plot_turb_lines()
    plot_exp_lines()
    ax = plt.axes()
    ax.set_aspect(2.0)
    plt.yticks(np.around(np.arange(-1.125, 1.126, 0.125), decimals=2))
    plt.tight_layout()
    if show:
        plt.show()
    if save:
        plt.savefig("figures/meancontquiv"+savetype)

def Vector(plt, X, Y, data, parameter, label_flg = 0):
	# std_length,scaleの目安として、
	# 流速の場合、1, 10.0
	# 風応力の場合、0.1, 1.0
	std_length	= parameter[0]
	scale		= parameter[1]
	intvl		= parameter[2]
	qkeyx		= parameter[3]
	qkeyy		= parameter[4]
	data[np.where(abs(data)==np.inf)]=np.nan
	datax = data[0, :, :]
	datay = data[1, :, :]
	Q=plt.quiver(X[::intvl,::intvl],Y[::intvl,::intvl],datax[::intvl,::intvl],datay[::intvl,::intvl],angles='xy',scale=scale)
	if label_flg == 0:
		plt.quiverkey(Q, qkeyx, qkeyy, std_length, str(std_length))

	return plt, Q