def get_data(frame_i, frame, modes = default_modes):
    """ frame_i is ith frame, frame is frame number """

    # Load Data Files
    normalized_density = (fromfile("gasdens%d.dat" % frame).reshape(num_rad, num_theta)) / surface_density
    averagedDensity = np.average(normalized_density, axis = 1)

    vrad = (fromfile("gasvrad%d.dat" % frame).reshape(num_rad, num_theta))
    vtheta = (fromfile("gasvtheta%d.dat" % frame).reshape(num_rad, num_theta))

    vorticity = util.velocity_curl(vrad, vtheta, rad, theta, frame = ref_frame)
    vortensity = vorticity / normalized_density[1:, 1:]

    # Find Peak in Radial Profile (in Outer Disk)
    peak_rad, peak_density = find_peak(averagedDensity)
    min_rad, min_density = find_min(averagedDensity, peak_rad)

    # Gather Azimuthal Profiles
    num_profiles = 5
    spread = 1.0 * scale_height # half-width

    azimuthal_radii = np.linspace(peak_rad - spread, peak_rad + spread, num_profiles)
    azimuthal_indices = [np.searchsorted(rad, this_radius) for this_radius in azimuthal_radii]
    azimuthal_profiles = np.array([np.fft.fft(vortensity[azimuthal_index, :]) for azimuthal_index in azimuthal_indices])

    # Normalize by m = 0 mode (integral of density), Take Absolute Value
    azimuthal_profiles = np.array([np.abs(azimuthal_profile / azimuthal_profile[0]) for azimuthal_profile in azimuthal_profiles])

    for m, mode in enumerate(modes):
        modes_over_time[m, frame_i] = np.max(azimuthal_profiles[:, mode])

    print "%d: %.4f, %.4f, %.4f, %.4f, %.4f" % (frame, np.max(azimuthal_profiles[:, 1]), np.max(azimuthal_profiles[:, 2]), np.max(azimuthal_profiles[:, 3]), np.max(azimuthal_profiles[:, 4]), np.max(azimuthal_profiles[:, 5]))

def find_extrema(args):
    # Unwrap Args
    i, frame = args

    # Data
    density = (fromfile("gasdens%d.dat" % frame).reshape(num_rad, num_theta)) / surface_density
    averagedDensity = np.average(density, axis = 1)

    vrad = (fromfile("gasvrad%d.dat" % frame).reshape(num_rad, num_theta))
    vtheta = (fromfile("gasvtheta%d.dat" % frame).reshape(num_rad, num_theta))

    vorticity = util.velocity_curl(vrad, vtheta, rad, theta, frame = ref_frame)
    vortensity = vorticity / density[1:, 1:]

    # Find Peak in Radial Profile (in Outer Disk)
    peak_rad, peak_rad_index, peak_density = find_radial_peak(averagedDensity)

    # Find Peaks in Azimuthal Profiles (and center around that peak)
    density_peak_theta, density_theta_index, max_density = find_azimuthal_extrema(density[peak_rad_index], maximum = True) # Max
    vortensity_peak_theta, vortensity_peak_theta_index, min_vortensity = find_azimuthal_extrema(vortensity[peak_rad_index], maximum = False) # Min

    print "%d, %.2f, %.3f" % (frame, max_density, min_vortensity)

    #return max_density, min_vortensity
    maximum_densities[i] = max_density
    minimum_vortensities[i] = min_vortensity

def get_data(frame):
    # Load Data Files
    normalized_density = (fromfile("gasdens%d.dat" % frame).reshape(num_rad, num_theta)) / surface_density
    averagedDensity = np.average(normalized_density, axis = 1)

    vrad = (fromfile("gasvrad%d.dat" % frame).reshape(num_rad, num_theta))
    vtheta = (fromfile("gasvtheta%d.dat" % frame).reshape(num_rad, num_theta))

    vorticity = util.velocity_curl(vrad, vtheta, rad, theta, frame = ref_frame)
    vortensity = vorticity / normalized_density[1:, 1:]

    # Find Peak in Radial Profile (in Outer Disk)

    peak_rad, peak_density = find_peak(averagedDensity)
    min_rad, min_density = find_min(averagedDensity, peak_rad)

    # Gather Azimuthal Profiles
    num_profiles = 5
    spread = 2.0 * scale_height # half-width

    azimuthal_radii = np.linspace(peak_rad - spread, peak_rad + spread, num_profiles)
    azimuthal_indices = [np.searchsorted(rad, this_radius) for this_radius in azimuthal_radii]
    azimuthal_profiles = [vortensity[azimuthal_index, :] for azimuthal_index in azimuthal_indices]

    return azimuthal_radii, azimuthal_profiles

def sum_vorticity(args):
    i, frame = args

    # Get Data
    density = (fromfile("gasdens%d.dat" % frame).reshape(num_rad, num_theta)) / surface_density_zero
    vrad = np.array(fromfile("gasvrad%d.dat" % frame).reshape(num_rad, num_theta))
    vtheta = np.array(fromfile("gasvtheta%d.dat" % frame).reshape(num_rad, num_theta))

    # Get Background Data
    vtheta_keplerian = np.array(fromfile("gasvtheta0.dat").reshape(num_rad, num_theta))

    # Subtract off Keplerian velocity (and rotate back into non-rotating frame???)
    vtheta -= (vtheta_keplerian)
    vorticity = util.velocity_curl(vrad, vtheta, rad, theta, frame = 0)

    #vorticity = np.abs(vorticity) # should be all negative

    # Add up vorticity
    dr = rad[1] - rad[0] # assumes arithmetic grid
    d_phi = theta[1] - theta[0]

    radial_vorticity = np.average(vorticity, axis = 1)
    total_vorticity = np.sum((dr * d_phi) * rad[:-1, None] * radial_vorticity)
    
    # Print Update
    print "%d: %.4f" % (frame, total_vorticity)

    # Store Data
    vorticity_over_time[i] = total_vorticity

def find_vortensity_min(frame):
    """ returns radius with maximum azimuthally averaged density """
    i = frame

    # Find density max (Search for vortensity minimum only near the density max)
    density_max = find_density_max(frame)
    start = density_max - (0.2 * (scale_height / 0.06))
    stop = density_max + (0.1 * (scale_height / 0.06))

    # Data
    truncated_rad = truncate(rad, start = start, stop = stop)

    density = (fromfile("gasdens%d.dat" % i).reshape(num_rad, num_theta))
    normalized_density = density / surface_density_zero

    vrad = (fromfile("gasvrad%d.dat" % i).reshape(num_rad, num_theta))
    vtheta = (fromfile("gasvtheta%d.dat" % i).reshape(num_rad, num_theta))

    vorticity = curl(vrad, vtheta, rad, theta)
    vortensity = vorticity / normalized_density[1:, 1:]
    avg_vortensity = truncate(np.average(vortensity, axis = 1), start = start, stop = stop)

    kernel_size = 1
    smoothed_avg_vortensity = smooth(avg_vortensity, kernel_size)

    # Retrieve radius with min value
    arg_min = np.argmin(smoothed_avg_vortensity)
    radius_min = truncated_rad[arg_min]

    return radius_min

def get_vortensity_variation(frame, vortex_location, figure = False):
    """
    return variation in vortensity (defined to be min / avg at a particular radius)
    """
    i = frame
    # Data
    density = (fromfile("gasdens%d.dat" % i).reshape(num_rad, num_theta))
    normalized_density = density / surface_density_zero

    vrad = (fromfile("gasvrad%d.dat" % i).reshape(num_rad, num_theta))
    vtheta = (fromfile("gasvtheta%d.dat" % i).reshape(num_rad, num_theta))

    vorticity = curl(vrad, vtheta, rad, theta)
    vortensity = vorticity / normalized_density[1:, 1:]

    # Azimuthal Profile (Average over width of vortex)
    arg_vortex = np.searchsorted(rad, vortex_location)
    half_width = int(16 * (num_rad / 512.0))
    kernel_size = num_theta / 200
    vortensity_at_vortex = smooth(np.average(vortensity[(arg_vortex - half_width):(arg_vortex + half_width), :], 
                                             weights = weights(half_width), axis = 0), kernel_size)
    # Replace values above maximum with maximum
    max_value = 0.3
    vortensity_at_vortex[vortensity_at_vortex > max_value] = max_value

    # Plot
    if figure:
        fig = plot.figure()

        # Axis
        angles = np.linspace(0, 2 * np.pi, 7)
        degree_angles = ["%d" % d_a for d_a in np.linspace(0, 360, 7)]

        plot.xlim(0, 2 * np.pi)
        plot.xticks(angles, degree_angles)

        # Plot
        plot.plot(theta[1:], vortensity_at_vortex, linewidth = linewidth)

        # Annotate
        plot.xlabel("Theta", fontsize = fontsize)
        plot.ylabel("Vortensity", fontsize = fontsize)

        # Save and Close
        directory = "azimuthalVortensity"

        plot.savefig("%s/azimuthalVortensity_%04d.png" % (directory, i), bbox_inches = 'tight', dpi = my_dpi)
        #plot.show()
        plot.close(fig) # Close Figure (to avoid too many figures)

    # Variation
    min_vortensity = np.min(vortensity_at_vortex)
    avg_vortensity = np.average(vortensity_at_vortex)

    #print min_vortensity

    variation = (1 + avg_vortensity - min_vortensity) / (2 * avg_vortensity)
    return variation

    def choose_axis(i, axis):
        # Orbit Number
        time = float(fargo_par["Ninterm"]) * float(fargo_par["DT"])
        orbit = int(round(time / (2 * np.pi), 0)) * i

        # Set up figure
        fig = plot.figure(figsize = (700 / my_dpi, 600 / my_dpi), dpi = my_dpi)
        ax = fig.add_subplot(111)

        # Axis
        angles = np.linspace(0, 2 * np.pi, 7)
        degree_angles = ["%d" % d_a for d_a in np.linspace(0, 360, 7)]

        plot.ylim(0, 2 * np.pi)
        plot.yticks(angles, degree_angles)

        if axis == "zoom":
            x = (rad - 1) / scale_height
            prefix = "zoom_"
            plot.xlim(0, 40) # to match the ApJL paper
            xlabel = r"($r - r_p$) $/$ $h$"
        else:
            x = rad
            prefix = ""
            plot.xlim(float(fargo_par["Rmin"]), float(fargo_par["Rmax"]))
            xlabel = "Radius"

        # Data
        vrad = np.array(fromfile("gasvrad%d.dat" % i).reshape(num_rad, num_theta))
        vtheta = np.array(fromfile("gasvtheta%d.dat" % i).reshape(num_rad, num_theta))

        vtheta_keplerian = np.array(fromfile("gasvtheta0.dat").reshape(num_rad, num_theta))

        # Subtract off Keplerian velocity (and rotate back into non-rotating frame???)
        vtheta -= (vtheta_keplerian)

        vorticity = util.velocity_curl(vrad, vtheta, rad, theta, frame = 0)

        # Divide out Angular Frequency (for Rossby Number)
        vorticity = vorticity / (np.array([r**(-1.5) for r in rad[:-1]]))[:, None]

        ### Plot ###
        result = ax.pcolormesh(x, theta, np.transpose(vorticity), cmap = cmap)
    
        fig.colorbar(result)
        result.set_clim(clim[0], clim[1])

        # Annotate
        plot.xlabel(xlabel, fontsize = fontsize)
        plot.ylabel(r"$\phi$", fontsize = fontsize)
        plot.title("Vorticity (v - vK) Map at Orbit %d" % orbit, fontsize = fontsize + 1)

        # Save and Close
        plot.savefig("%s/%sexcessVorticityMap_%03d.png" % (save_directory, prefix, i), bbox_inches = 'tight', dpi = my_dpi)
        if show:
            plot.show()
        plot.close(fig) # Close Figure (to avoid too many figures)

    def choose_axis(i, axis):
        # Orbit Number
        time = float(fargo_par["Ninterm"]) * float(fargo_par["DT"])
        orbit = int(round(time / (2 * np.pi), 0)) * i

        # Set up figure
        fig = plot.figure(figsize = (700 / my_dpi, 600 / my_dpi), dpi = my_dpi)
        ax = fig.add_subplot(111)

        # Axis
        angles = np.linspace(0, 2 * np.pi, 7)
        degree_angles = ["%d" % d_a for d_a in np.linspace(0, 360, 7)]

        plot.ylim(0, 2 * np.pi)
        plot.yticks(angles, degree_angles)

        if axis == "zoom":
            x = (rad - 1) / scale_height
            prefix = "zoom_"
            plot.xlim(0, 40) # to match the ApJL paper
            xlabel = r"($r - r_p$) $/$ $h$"
        else:
            x = rad
            prefix = ""
            plot.xlim(float(fargo_par["Rmin"]), float(fargo_par["Rmax"]))
            xlabel = "Radius"

        # Data
        density = (fromfile("gasdens%d.dat" % i).reshape(num_rad, num_theta))
        normalized_density = density / surface_density_zero

        vrad = (fromfile("gasvrad%d.dat" % i).reshape(num_rad, num_theta))
        vtheta = (fromfile("gasvtheta%d.dat" % i).reshape(num_rad, num_theta))

        vortensity = util.velocity_curl(vrad, vtheta, rad, theta, frame = ref_frame) / normalized_density[1:, 1:]
        avg_vortensity = np.average(vortensity, axis = 1)

        excess_vortensity = avg_vortensity[:, None] - vortensity #### Vortex is Positive ####

        ### Plot ###
        result = ax.pcolormesh(x, theta, np.transpose(excess_vortensity), cmap = cmap)
    
        fig.colorbar(result)
        result.set_clim(clim[0], clim[1])

        # Annotate
        plot.xlabel(xlabel, fontsize = fontsize)
        plot.ylabel(r"$\phi$", fontsize = fontsize)
        plot.title("Excess Vortensity Map at Orbit %d" % orbit, fontsize = fontsize + 1)

        # Save and Close
        plot.savefig("%s/%sexcessVortensityMap_%03d.png" % (save_directory, prefix, i), bbox_inches = 'tight', dpi = my_dpi)
        if show:
            plot.show()
        plot.close(fig) # Close Figure (to avoid too many figures)

def sum_vorticity(args):
    i, frame = args

    # Get Data
    density = (fromfile("gasdens%d.dat" % frame).reshape(num_rad, num_theta)) / surface_density_zero
    vrad = np.array(fromfile("gasvrad%d.dat" % frame).reshape(num_rad, num_theta))
    vtheta = np.array(fromfile("gasvtheta%d.dat" % frame).reshape(num_rad, num_theta))

    # Get Background Data
    vtheta_keplerian = np.array(fromfile("gasvtheta0.dat").reshape(num_rad, num_theta))

    # Subtract off Keplerian velocity (and rotate back into non-rotating frame???)
    vtheta -= (vtheta_keplerian)
    vorticity = util.velocity_curl(vrad, vtheta, rad, theta, frame = 0)

    # Mask Non-Vortex Regions
    min_vorticity = -0.65

    vorticity[density[:-1, :-1] < 0.45] = 0 # if density is too low, it's not in the vortex
    vorticity[vorticity > 0] = 0 # if vorticity is positive, it's not in the vortex
    vorticity[vorticity < min_vorticity] = min_vorticity # if vorticity is too low, it's not in the vortex

    vorticity = np.abs(vorticity) # everything that remains should be negative. switch to positive.

    # Extract Near Vortex
    averagedDensity = np.average(density, axis = 1)
    peak_rad, peak_density = find_peak(averagedDensity)

    vortex_start = np.max([1.0, peak_rad - 5.0 * scale_height])
    vortex_end = peak_rad + 5.0 * scale_height

    vortex_start_i = np.searchsorted(rad, vortex_start)
    vortex_end_i = np.searchsorted(rad, vortex_end)

    vortex_rad = rad[vortex_start_i : vortex_end_i]
    vortex_vorticity_grid = vorticity[vortex_start_i : vortex_end_i]

    vortex_vorticity = np.average(vortex_vorticity_grid, axis = 1)

    # Add up vorticity
    dr = rad[1] - rad[0] # assumes arithmetic grid
    d_phi = theta[1] - theta[0]

    total_vorticity = np.sum((dr * d_phi) * vortex_rad[:, None] * vortex_vorticity)

    ##### Instead: take max vorticity (90th percentile???) #####
    
    # Print Update
    print "%d: %.4f" % (frame, total_vorticity)

    # Store Data
    vorticity_over_time[i] = total_vorticity

def gather_vortex_over_time():
    """ add up vortex mass over time """    
    for frame in times:
        density = (fromfile("gasdens%d.dat" % frame).reshape(num_rad, num_theta))
        normalized_density = density / surface_density_zero

        vrad = (fromfile("gasvrad%d.dat" % frame).reshape(num_rad, num_theta))
        vtheta = (fromfile("gasvtheta%d.dat" % frame).reshape(num_rad, num_theta))

        vorticity = curl(vrad, vtheta, rad, theta)
        vortensity = vorticity / normalized_density[1:, 1:]

        vortex_mass_i, _ = vortex_mass(rad, theta, density, vortensity)
        vortex_masses.append(vortex_mass_i)

def map_one_vortex(frame):
    """ get 2-D grid showing vortex """
    density = (fromfile("gasdens%d.dat" % frame).reshape(num_rad, num_theta))
    normalized_density = density / surface_density_zero

    vrad = (fromfile("gasvrad%d.dat" % frame).reshape(num_rad, num_theta))
    vtheta = (fromfile("gasvtheta%d.dat" % frame).reshape(num_rad, num_theta))

    vorticity = curl(vrad, vtheta, rad, theta)
    vortensity = vorticity / normalized_density[1:, 1:]

    this_vortex_mass, vortex_map = vortex_mass(rad, theta, density, vortensity)

    print "Vortex Mass at Frame %d: %.8f" % (frame, this_vortex_mass)
    return vortex_map

    def choose_axis(i, axis):
        # Orbit Number
        time = float(fargo_par["Ninterm"]) * float(fargo_par["DT"])
        orbit = int(round(time / (2 * np.pi), 0)) * i

        # Set up figure
        fig = plot.figure(figsize = (700 / my_dpi, 600 / my_dpi), dpi = my_dpi)
        ax = fig.add_subplot(111)

        # Axis
        angles = np.linspace(0, 2 * np.pi, 7)
        degree_angles = ["%d" % d_a for d_a in np.linspace(0, 360, 7)]

        plot.ylim(0, 2 * np.pi)
        plot.yticks(angles, degree_angles)
        if axis == "zoom":
            x = (rad - 1) / scale_height
            prefix = "zoom_"
            plot.xlim(0, 40) # to match the ApJL paper
            xlabel = r"($r - r_p$) $/$ $h$"
        else:
            x = rad
            prefix = ""
            plot.xlim(float(fargo_par["Rmin"]), float(fargo_par["Rmax"]))
            xlabel = "Radius"

        # Data
        blank_density = (fromfile("gasdens%d.dat" % blank_frame).reshape(num_rad, num_theta))
        density = (fromfile("gasdens%d.dat" % i).reshape(num_rad, num_theta))

        normalized_density = (density - blank_density) / surface_density_zero # Diff

        ### Plot ###
        result = ax.pcolormesh(x, theta, np.transpose(normalized_density), cmap = cmap)
        fig.colorbar(result)
        result.set_clim(clim[0], clim[1])

        # Annotate
        this_title = readTitle()
        plot.xlabel(xlabel, fontsize = fontsize)
        plot.ylabel(r"$\phi$", fontsize = fontsize)
        plot.title("Gas Difference Density Map at Orbit %d (- Orbit %d)\n%s" % (orbit, blank_frame, this_title), fontsize = fontsize + 1)

        # Save and Close
        plot.savefig("%s/%sdiff_densityMap_%04d-%04d.png" % (save_directory, prefix, blank_frame, i), bbox_inches = 'tight', dpi = my_dpi)
        if show:
            plot.show()
        plot.close(fig) # Close Figure (to avoid too many figures)

def make_plot(frame, show = False):
    # Orbit Number
    time = float(fargo_par["Ninterm"]) * float(fargo_par["DT"])
    orbit = int(round(time / (2 * np.pi), 0)) * frame

    # Set up figure
    fig = plot.figure()
    ax = fig.add_subplot(111, polar = True) 

    # Axis
    rmax = 1.0 # to match ApJL paper

    plot.xticks([], []) # Angles
    plot.yticks([rmax], ["%.1f" % rmax]) # Max Radius
    plot.ylim(0, rmax)

    # Data
    density = (fromfile("gasdens%d.dat" % frame).reshape(num_rad, num_theta))
    normalized_density = density / surface_density_zero

    ### Plot ###
    result = ax.pcolormesh(theta, rad, normalized_density, cmap = cmap)
    fig.colorbar(result)
    result.set_clim(clim[0], clim[1])

    # Annotate
    this_title = readTitle()
    plot.title("Gas Density Map at Orbit %d: %s" % (orbit, this_title), fontsize = fontsize + 1)

    # Save and Close
    plot.savefig("%s/innerDensityMap_%04d.png" % (save_directory, frame), bbox_inches = 'tight', dpi = my_dpi)
    if show:
        plot.show()
    plot.close(fig) # Close Figure (to avoid too many figures)

def get_data(frame):
    # Find Peak in Radial Profile (in Outer Disk)
    density = (fromfile("gasdens%d.dat" % frame).reshape(num_rad, num_theta)) / surface_density
    averagedDensity = np.average(density, axis = 1)

    peak_rad, peak_index, peak_density = find_peak(averagedDensity)
    min_rad, min_density = find_min(averagedDensity, peak_rad)

    peak_theta, peak_theta_index, peak_azimuthal_density = find_azimuthal_peak(density[peak_index])

    if len(sys.argv) > 2:
        # Supply central theta as an argument
        peak_theta = float(sys.argv[2])

    # Gather Azimuthal Profiles
    pressure = density * (scale_height)**2 * np.power(rad, -1)[:, None]

    num_profiles = 5
    spread = 30.0 # half-width

    radial_theta = np.linspace(peak_theta - spread, peak_theta + spread, num_profiles)
    radial_theta[radial_theta < 0] += 360.0
    radial_theta[radial_theta > 360] -= 360.0

    radial_indices = [np.searchsorted(theta, this_theta * (np.pi / 180.0)) for this_theta in radial_theta]
    radial_profiles = [pressure[:, radial_index] for radial_index in radial_indices]

    return radial_theta, radial_profiles

def read_nse(fin):
  """Read single electrode spike record.
  Inputs:
    fin - file handle
    only_timestamps - if true, only load the timestamps, ignoring the waveform and feature data

  Output: Dictionary with fields
    'header'  - header info
    'packets' - pylab data structure with the spike data

  Notes:
    0. What is spike acquizition entity number? Ask neuralynx
    1. Removing LOAD_ATTR overhead by defining time_stamp_append = [time_stamp[n].append for n in xrange(100)] and using
       time_stamp_append[dwCellNumber](qwTimeStamp) in the loop does not seem to improve performance.
       It reduces readability so I did not use it.
    2. Disabling garbage collection did not help
    3. In general, dictionary lookups slow things down
    4. using numpy arrays, with preallocation is slower than the dumb, straightforward python list appending
  """

  hdr = read_header(fin)
  nse_packet = pylab.dtype([
    ('timestamp', 'Q'),
    ('saen', 'I'),
    ('cellno', 'I'),
    ('Features', '8I'),
    ('waveform', '32h')
  ])
  data = pylab.fromfile(fin, dtype=nse_packet, count=-1)
  return {'header': hdr, 'packets': data}

def get_data(frame):
    # Find Peak in Radial Profile (in Outer Disk)
    density = (fromfile("gasdens%d.dat" % frame).reshape(num_rad, num_theta)) / surface_density
    averagedDensity = np.average(density, axis = 1)

    peak_rad, peak_density = find_peak(averagedDensity)
    min_rad, min_density = find_min(averagedDensity, peak_rad)

    ### Gather Azimuthal Profiles ###
    num_profiles = 5
    spread = 1.0 * scale_height # half-width

    azimuthal_radii = np.linspace(peak_rad - spread, peak_rad + spread, num_profiles)
    azimuthal_indices = [np.searchsorted(rad, this_radius) for this_radius in azimuthal_radii]
    azimuthal_profiles = np.array([np.fft.fft(density[azimuthal_index, :]) for azimuthal_index in azimuthal_indices])

    # Normalize by m = 0 mode (integral of density), Take Absolute Value
    azimuthal_profiles = [np.abs(azimuthal_profile / azimuthal_profile[0]) for azimuthal_profile in azimuthal_profiles]

    ### Gather Averaged Profiles ###
    start_half = azimuthal_indices[1]
    end_half = azimuthal_indices[-2]
    avg_half_profile = np.average(density[start_half : end_half, :], axis = 0)
    avg_half = np.fft.fft(avg_half_profile)

    start_full = azimuthal_indices[0]
    end_full = azimuthal_indices[-1]
    avg_full_profile = np.average(density[start_full : end_full, :], axis = 0)
    avg_full = np.fft.fft(avg_full_profile)

    # Normalize
    avg_half = np.abs(avg_half / avg_half[0])
    avg_full = np.abs(avg_full / avg_full[0])

    return azimuthal_radii, azimuthal_profiles, avg_half, avg_full

    def choose_axis(i, axis):
        # Orbit Number
        time = float(fargo_par["Ninterm"]) * float(fargo_par["DT"])
        orbit = int(round(time / (2 * np.pi), 0)) * i

        # Axis
        if axis == "zoom":
            x = (rad - 1) / scale_height
            prefix = "zoom_"
            plot.xlim(0, 40) # to match the ApJL paper
            #plot.ylim(0, 1.2)
            xlabel = r"($r - r_p$) $/$ $h$"
        else:
            x = rad
            prefix = ""
            plot.xlim(float(fargo_par["Rmin"]) - 0.05, float(fargo_par["Rmax"]) + 0.05)
            #plot.ylim(0, 5.0)
            xlabel = "Radius"
            
        # Data
        vrad = (fromfile("gasvrad%d.dat" % i).reshape(num_rad, num_theta))
        averaged_vrad = np.average(vrad, axis = 1)

        ### Plot ###
        plot.plot(x, averaged_vrad, linewidth = linewidth, label = "%d" % frame)

        # Annotate
        this_title = readTitle()
        plot.xlabel(xlabel, fontsize = fontsize)
        plot.ylabel(r"Azimuthally Averaged $V_{rad}$", fontsize = fontsize)
        plot.title("%s" % this_title, fontsize = fontsize + 1)

def get_data(frame_i, frame, modes = default_modes):
    """ frame_i is ith frame, frame is frame number """

    # Find Peak in Radial Profile (in Outer Disk)
    density = (fromfile("gasdens%d.dat" % frame).reshape(num_rad, num_theta)) / surface_density
    averagedDensity = np.average(density, axis = 1)

    peak_rad, peak_density = find_peak(averagedDensity)
    min_rad, min_density = find_min(averagedDensity, peak_rad)

    # Gather Azimuthal Profiles
    num_profiles = 7
    spread = 1.0 * scale_height # half-width

    azimuthal_radii = np.linspace(peak_rad - spread, peak_rad + spread, num_profiles)
    azimuthal_indices = [np.searchsorted(rad, this_radius) for this_radius in azimuthal_radii]
    azimuthal_profiles = np.array([np.fft.fft(density[azimuthal_index, :]) for azimuthal_index in azimuthal_indices])

    # Normalize by m = 0 mode (integral of density), Take Absolute Value
    azimuthal_profiles = np.array([np.abs(azimuthal_profile / azimuthal_profile[0]) for azimuthal_profile in azimuthal_profiles])

    for m, mode in enumerate(modes):
        modes_over_time[m, frame_i] = np.max(azimuthal_profiles[:, mode]) #np.sqrt(np.sum(np.power(azimuthal_profiles[:, mode], 2))) 

    single_mode_strength[frame_i] = modes_over_time[0, frame_i] / np.max(modes_over_time[1:, frame_i]) # m = 1 / Max of Higher Number Modes
    single_mode_concentration[frame_i] = np.std(azimuthal_profiles[:, 1]) / modes_over_time[0, frame_i]

    print "%d: %.4f, %.4f, %.4f, %.4f, %.4f - [%.4f, %.4f]" % (frame, np.max(azimuthal_profiles[:, 1]), np.max(azimuthal_profiles[:, 2]), np.max(azimuthal_profiles[:, 3]), np.max(azimuthal_profiles[:, 4]), np.max(azimuthal_profiles[:, 5]), single_mode_strength[frame_i], single_mode_concentration[frame_i])

    def choose_axis(i, axis):
        # Orbit Number
        time = float(fargo_par["Ninterm"]) * float(fargo_par["DT"])
        orbit = int(round(time / (2 * np.pi), 0)) * i

        # Set up figure
        fig = plot.figure(figsize = (700 / my_dpi, 600 / my_dpi), dpi = my_dpi)

        # Axis
        plot.ylim(0, 1.3)
        if axis == "zoom":
            x = (rad - 1) / scale_height
            prefix = "zoom_"
            plot.xlim(0, 40) # to match the ApJL paper
            #plot.ylim(0, 1.3)
            xlabel = r"($r - r_p$) $/$ $h$"
        else:
            x = rad
            prefix = ""
            plot.xlim(float(fargo_par["Rmin"]), float(fargo_par["Rmax"]))
            xlabel = "Radius"

        # Data
        density = (fromfile("gasdens%d.dat" % i).reshape(num_rad, num_theta))
        normalized_density = density / surface_density_zero

        vrad = (fromfile("gasvrad%d.dat" % i).reshape(num_rad, num_theta))
        vtheta = (fromfile("gasvtheta%d.dat" % i).reshape(num_rad, num_theta))

        vorticity = util.velocity_curl(vrad, vtheta, rad, theta, average = True, frame = ref_frame)
        vortensity = vorticity / normalized_density[1:, 1:]
        averaged_w = np.average(vortensity, axis = 1)

        ### Plot ###
        plot.plot(x[1:], averaged_w, linewidth = linewidth)

        # Annotate
        this_title = readTitle()
        plot.xlabel(xlabel, fontsize = fontsize)
        plot.ylabel("Azimuthally Averaged Vortensity", fontsize = fontsize)
        plot.title("Orbit %d: %s" % (orbit, this_title), fontsize = fontsize + 1)

        # Save and Close
        plot.savefig("%s/%saveragedVortensity_%03d.png" % (save_directory, prefix, i), bbox_inches = 'tight', dpi = my_dpi)
        if show:
            plot.show()
        plot.close(fig) # Close Figure (to avoid too many figures)

def make_plot(frame, show = False):
    # Orbit Number
    time = float(fargo_par["Ninterm"]) * float(fargo_par["DT"])
    orbit = int(round(time / (2 * np.pi), 0)) * frame

    # Set up figure
    fig = plot.figure()
    ax = fig.add_subplot(111)

    # Axis
    angles = np.linspace(0, 2 * np.pi, 7)
    degree_angles = ["%d" % d_a for d_a in np.linspace(0, 360, 7)]

    plot.ylim(0, 2 * np.pi)
    plot.yticks(angles, degree_angles)

    plot.xlim(float(fargo_par["Rmin"]), 1.0)

    # Data
    x = rad
    density = (fromfile("gasdens%d.dat" % frame).reshape(num_rad, num_theta))
    normalized_density = density / surface_density_zero

    vrad = (fromfile("gasvrad%d.dat" % frame).reshape(num_rad, num_theta))
    vtheta = (fromfile("gasvtheta%d.dat" % frame).reshape(num_rad, num_theta))

    vorticity = curl(vrad, vtheta, rad, theta) / normalized_density[1:, 1:]

    ### Plot ###
    result = ax.pcolormesh(x, theta, np.transpose(vorticity), cmap = cmap)

    fig.colorbar(result)
    result.set_clim(clim[0], clim[1])

    # Annotate
    this_title = readTitle()
    plot.xlabel("Radius", fontsize = fontsize)
    plot.ylabel(r"$\phi$", fontsize = fontsize)
    plot.title("Vortensity Map at Orbit %d: %s" % (orbit, this_title), fontsize = fontsize + 1)

    # Save and Close
    plot.savefig("%s/vorticityMap_%03d.png" % (save_directory, frame), bbox_inches = 'tight', dpi = my_dpi)
    if show:
        plot.show()
    plot.close(fig) # Close Figure (to avoid too many figures)