def main():
    re = []
    rm = []
    t = []
    ve = []
    vm = []
    earth = Planet()
    earth.set_e(0.0167)
    earth.set_rmin(1.471E11)
    earth.set_period(365*86400)
    earth.set_semi_major_axis(149600000)
    mars = Planet()
    mars.set_e(0.09341233)
    mars.set_rmin(2.0662E11)
    mars.set_period(687*86400)
    mars.set_semi_major_axis(227920000)
    
    for x in range(0, int(sys.argv[1])):
        re.append(earth.step())
        rm.append(mars.step())
        ve.append(earth.velocity())
        vm.append(mars.velocity())
        t.append(earth.theta)
        sys.stdout.write("\rStep " + str(earth.count) + " / " + sys.argv[1] + " " + str((float(earth.count)/float(sys.argv[1]))*100) + "% " + str(earth.theta) + " " + str(mars.theta))
        sys.stdout.flush()
 	
    print "\n"
    plt.figure(1)
    plt.polar(t, re)
    plt.polar(t, rm)

    plt.figure(2)
    plt.plot(t, ve)
    plt.plot(t, vm)
    plt.show()

def plot_phaseplot_lpf(dictdata, keys, autok, title, withn='yes'):
    colors = ['r', 'b', 'g', 'y', 'k']
    
    plt.suptitle(title, fontsize='large' )
    if autok == 'yes':
        k = dictdata.keys()
    
    for i, condition in enumerate(keys):
        datac = colors[i]
        data = dictdata[condition]
        
        try:
            n = len(data)
            theta, r = zip(*data)
        except TypeError:
            theta, r = data
            n = 1
        if withn == 'yes':
            plt.polar(theta, r, 'o', color=datac, label=condition + '\n n=' + str(n))
        if withn == 'no':
            plt.polar(theta, r, 'o', color=datac, label=condition)

        #lines, labels = plt.rgrids( (1.0, 1.4), ('', ''), angle=0 )
        tlines, tlabels = plt.thetagrids( (0, 90, 180, 270), ('0', 'pi/2', 'pi', '3pi/2') )
        leg = plt.legend(loc=(0.95,0.75))
  
        for t in leg.get_texts():
            t.set_fontsize('small')
            
        plt.subplots_adjust(top=0.85)
        plt.draw()

 def plot_SParameters(self, S='S21' , style='MA'):
     '''plot  S parameter from data
         
         Input:
             - S (String) ["S11", "S12", "S21", "S22"]: Set which parameter you want to plot
             - asked_format ("String") ["MA", "DB", "RI"]: Set in which format we would like to have plot
         
         Output:
             - matplotlib 2d figure
     '''
     
     x, y, z = self.get_SParameters(S, style)
     
     factor = self.get_frequency_unit(style='Float')
 
     if factor == 1. :
         x_label = 'Frequency [Hz]'
     elif factor == 1e3 :
         x_label = 'Frequency [kHz]'
     elif factor == 1e6 :
         x_label = 'Frequency [MHz]'
     elif factor == 1e9 :
         x_label = 'Frequency [GHz]'
     
     fig = plt.figure()
     
     
     if re.match('^[mM][aA]$', style) :
         
         ax1 = fig.add_subplot(211)
         ax1.plot(x, y, label=S)
         plt.ylabel('Amplitude [V]')
         plt.grid()
     elif re.match('^[dD][bB]$', style) :
         
         ax1 = fig.add_subplot(211)
         ax1.plot(x, y, label=S)
         plt.ylabel('Attenuation [dB]')
         plt.grid()
     
     
     if re.match('^[dD][bB]$|^[mM][aA]$', style):
     
         ax2 = fig.add_subplot(212, sharex = ax1)
         ax2.plot(x, z, label=S)
         plt.ylabel('Phase [deg]')
         plt.xlabel(x_label)
     elif re.match('^[rR][iI]$', style):
     
         ax1 = fig.add_subplot(111)
         y = y*1./y.max()
         plt.polar(y, z, label=S)
         plt.title('Normalised amplitude')
         plt.ylabel('Real part')
         plt.ylabel('Imaginary part')
     
     
     plt.legend(loc='best')
     plt.grid()
     plt.show()

def test_polar_units():
    import matplotlib.testing.jpl_units as units
    from nose.tools import assert_true
    units.register()

    pi = np.pi
    deg = units.UnitDbl( 1.0, "deg" )
    km = units.UnitDbl( 1.0, "km" )

    x1 = [ pi/6.0, pi/4.0, pi/3.0, pi/2.0 ]
    x2 = [ 30.0*deg, 45.0*deg, 60.0*deg, 90.0*deg ]

    y1 = [ 1.0, 2.0, 3.0, 4.0]
    y2 = [ 4.0, 3.0, 2.0, 1.0 ]

    fig = plt.figure()

    plt.polar( x2, y1, color = "blue" )

    # polar( x2, y1, color = "red", xunits="rad" )
    # polar( x2, y2, color = "green" )

    fig = plt.figure()

    # make sure runits and theta units work
    y1 = [ y*km for y in y1 ]
    plt.polar( x2, y1, color = "blue", thetaunits="rad", runits="km" )
    assert_true( isinstance(plt.gca().get_xaxis().get_major_formatter(), units.UnitDblFormatter) )

 def drawPair(self, pair, label):
     start_angle = self.hourToAngle(pair[1])
     end_angle = self.hourToNextAngle(pair[0], start_angle)
     new_angles = np.linspace(start_angle, end_angle, 20)
     new_points = np.ones(len(new_angles))
     plt.polar(new_angles, new_points)
     plt.fill_between(new_angles, new_points, facecolor='yellow', alpha=0.5)
     self.drawLabel(new_angles, new_points, label)

def floret_revolutions(n, **kwargs):
  r = np.arange(n)
  F = fibonacci_lim(r[-1])
  plt.figure(figsize=(6,6))
  plt.polar(np.mod(r/phi,1)*2*np.pi, r, '.', **kwargs)
  plt.polar(np.mod(F/phi,1)*2*np.pi, F, 'r.', **kwargs)
  plt.gca().set_rticks([])
  plt.show()

def display_picture(points):
    """Graphs the points given on a polar coordinate grid"""
    r_vals = []
    theta_vals = []
    for r, theta in points:
        r_vals.append(r)
        theta_vals.append(theta)

    pyplot.polar(r_vals, theta_vals, linestyle='solid', marker='o')
    pyplot.show()

	def plot_orbit(self):
		ta = self.trueAnom()
		sma = self.smAxis()[0]
		e = self.ecc()
		theta = np.linspace(0,2*math.pi,360)
		r=(sma*(1-e**2))/(1+e*np.cos(theta))
		plt.polar(theta, r)
		print(np.c_[r,theta])
		plt.savefig("Orbit.png")
		plt.show()

def my(heart):
    r, t = heart
    plt.polar(r, t, 'r', lw=5)
    tick_params = {'axis':'y', 'which':'both', 
                   'bottom': False, 'top':False, 
                   'left': False, 'right': False,
                   'labelbottom': False, 'labelleft': False}
                   
    plt.tick_params(**tick_params)
    plt.show()

 def plot_diff_cross_sect(self, **kwargs):
     """Displays a plot of the differential cross section.
     
     Arguments:
        **kwargs: Any additional arguments to plot
     Returns:
        Displays the plot
     """
     
     tgrid = np.linspace(-np.pi,np.pi,361)
     plt.polar(tgrid,self.diff_cross_sect(tgrid), marker='.', **kwargs)
     plt.show()

def show_floret_growth(theta, r, prefix='frame', grow=True, ms=10):
  n = len(r)
  savefmt = prefix + '%0' + str(len(str(n-1))) + 'i.png'
  plt.figure(figsize=(6,6))
  for i in range(n):
    plt.clf()
    if grow: s = float(n - i - 1) / n
    else:    s = 0.
    plt.polar(theta[:i+1], r[:i+1]-s, '.', mec='k', mfc='b', ms=ms)
    plt.polar(theta[i], r[i]-s, '.', mec='k', mfc='r', ms=ms)
    plt.gca().set_rlim([0, 1])
    plt.gca().set_rticks([])
    plt.savefig(savefmt %i)

def main():
    earth = Body([0, 147090000000], 2929000, 5.9726E24, Orbit(0.016710219, 149600000000, 1.495583757E8, 365))
#    print earth.orbit.radius_from_angle(math.pi)
#    print earth.orbit.find_e_anomaly(271433.6, 0.016710219)
#    print earth.orbit.find_true_anomaly(earth.orbit.find_e_anomaly(271433.6, 0.016710219))
    theta = []
    radius = []
    for time in [3]:
        earth.orbit.calc_position(time)
        theta.append(earth.orbit.true_anomaly)
        radius.append(earth.orbit.radius)
    plt.polar(theta, radius, 'o')
    plt.show()

def polar_demo():
    """Make a polar plot of some random angles.

    """

    # Sample a bunch of random angles
    A = rand.randint(0, 360, 100)
    # Plot each angle as a point with radius 1
    plt.figure()
    plt.polar(A, np.ones_like(A), 'ko')
    # Disable y-ticks, because they are distracting
    plt.yticks([], [])
    plt.title("Polar Demo", fontsize=titlesize)

def plot_phase(control=None, state0=None):
	"""generates a phase portrait and phase trajectory from initial conditions"""
	#time = np.linspace(0, 20, 100) not needed for simulate()
	
	results, time = simulate(100, 10, control, state0)
	theta, h = results[:,0], results[:,1]
	state0 = (theta[0], h[0])
	#statew = [theta[-1], h[-1]]
	#print "Final: ", statew
	
	#system trajectory
	plot.figure(1) 
	plot.plot(theta, h, color='green') #use polar(theta, h)?
	plot.figure(2)
	plot.polar(theta, h, color='green')
		
	#phase portrait (vector field)
	thetamax, hmax = max(abs(theta)), max(abs(h))
	theta, h = numpy.meshgrid(numpy.linspace(-thetamax, thetamax, 10), 
				  numpy.linspace(-hmax, hmax, 10))
	Dtheta, Dh = numpy.array(theta), numpy.array(h)
	for i in range(theta.shape[0]):
		for j in range(theta.shape[1]):
			Dtheta[i,j] = dtheta(float(theta[i,j]), float(h[i,j]))
			Dh[i,j] = dh(float(theta[i,j]), float(h[i,j]))
	
	plot.figure(1)
	plot.quiver(theta, h, Dtheta, Dh) #no polar equivalent...
	
	#optimal path mode
	h = numpy.linspace(-hmax, hmax, 100)
	plot.plot(path(h), h, color='blue') #use polar(theta, h)?
	if control is None:
		plot.savefig("optimal-satellite-cart.png", dpi=200)
	else:
		plot.savefig("satellite-nonoptimal-cart.png", dpi=200)
	plot.xlabel("angle (rad)")
	plot.ylabel("angular momentum (N-m-s)")
	
	#optimal mode in polar
	plot.figure(2)
	plot.polar(path(h), h, color='blue')
	if control is None:
		plot.savefig("optimal-satellite-polar.png", dpi=200)
	else:
		plot.savefig("satellite-nonoptimal-polar.png", dpi=200)
	plot.xlabel("angle (rad)")
	plot.ylabel("angular momentum (N-m-s)")
	
	plot.plot(state0[0], state0[1], 'o', color='green')
	plot.show()

def relative_direction_distribution(xy,verbose=False):
    """computes instantaneous directions and make an histogram centered on previous direction
    """
    dxy = xy[1:,:]-xy[0:-1,:]
    theta = npy.arctan2(dxy[:,0],dxy[:,1])
    rho = npy.sqrt(npy.sum(dxy**2,axis=1))
    dtheta = theta[1:]-theta[0:-1]
    dtheta = npy.hstack(([0],dtheta))
    #verify that theta is in [-pi,+pi]
    clip_dtheta = dtheta.copy()
    clip_dtheta[dtheta>npy.pi] = dtheta[dtheta>npy.pi] - 2.*npy.pi
    clip_dtheta[dtheta<-npy.pi] = dtheta[dtheta<-npy.pi] + 2.* npy.pi

    #resulting direction and dispersion
    (R,V,Theta,Rtot) = rayleigh(rho,clip_dtheta)

    #distribution
    N = 8
    width = 2*npy.pi/N
    offset_th = .5*width
    bins = npy.linspace(-npy.pi-offset_th,npy.pi+offset_th,N+2,endpoint=True)
    h_theta,bin_theta = npy.histogram(clip_dtheta,bins=bins,weights=rho) # ! weighted histogram
    #grouping first and last bin corresponding to the same direction
    h_theta[0]+=h_theta[-1]

    if verbose:
        import matplotlib.pyplot as plt

        fig=plt.figure()
        ax = fig.add_axes([0.1, 0.1, 0.8, 0.8], polar=True,)
        #    #plot polar histogram bins
        #    ax.bar(bin_theta[0:-1], npy.ones_like(bin_theta[0:-1]), width=width, bottom=0.0,alpha=.5)

        #plot polar distribution
        ax.bar(bin_theta[0:-2], h_theta[:-1], width=width, bottom=0.0,label='rel.direction hist.')
        #plot relative displacement
        plt.polar(clip_dtheta,rho, 'yo',label='rel.direction')

        #plot main direction and dispersion
        plt.polar(Theta,R*Rtot,'ro',label='avg.')
        ax.bar(Theta-V/2, R*Rtot*.01, color='k',width=V, bottom=R*Rtot,label='dispersion')
        ax.legend(loc='upper left')
        #plot xy trajectory
        ax = fig.add_axes([0.1, 0.1, 0.2, 0.2])
        plt.plot(xy[:,0],xy[:,1])
        plt.plot(xy[0,0],xy[0,1],'k+')
        plt.show()
    return (R,V,Theta,Rtot,clip_dtheta,rho)

def outlinebrick(brick, polar=False):
    import matplotlib.pyplot as plt
    if polar:
        n = 20
        x = np.linspace(brick.ra_min, brick.ra_max, n) * np.pi / 180
        ra = np.concatenate([x, x[-1::-1]])
        dec = np.concatenate([np.ones(n)*brick.dec_min, np.ones(n)*brick.dec_max])
        if brick.dec_min < 0:
            dec = dec + 90
        else:
            dec = 90 - dec
        plt.polar(ra, dec, 'k-', lw=2)
    else:
        x = [brick.ra_min, brick.ra_max, brick.ra_max, brick.ra_min, brick.ra_min]
        y = [brick.dec_min, brick.dec_min, brick.dec_max, brick.dec_max, brick.dec_min]
        plt.plot(x, y, 'k-', lw=2)

def process_ellipse(normPower, theta1RadFinal, figWidth, figHeigth, dir, number):
    """
    :param normPower:
    :param theta1RadFinal:
    :param figWidth: width of the figure
    :param figHeigth: height of the figure
    :param dir: full path to the directory where one wants to store the intermediate images
    :param number:
    :return:
    """
    # Combine data into [XY] to fit to an ellipse
    Mirtheta1RadFinal1 = np.concatenate([theta1RadFinal.T, (theta1RadFinal + np.pi).T])
    MirnormPower = np.concatenate([normPower.T, normPower.T])

    # Convert mirrored polar coords to cartesian coords
    xdata, ydata = pol2cart(Mirtheta1RadFinal1, MirnormPower)
    ell_data = np.vstack([xdata, ydata])
    ell_data = ell_data.T

    # Python fitting function, see EllipseDirectFit
    A, centroid = EllipseDirectFit(ell_data)
    t = orientation(A)

    # Plot Lower Left - Polar plot of angular distribution
    angDist = plt.figure(figsize=(figWidth, figHeigth))  # Creates a figure containing angular distribution.
    r_line = np.arange(0, max(MirnormPower) + .5, .5)
    th = np.zeros(len(r_line))
    for i in range(0, len(r_line)):
        th[i] = t
    th = np.concatenate([th, (th + 180)])
    r_line = np.concatenate([r_line, r_line])
    plt.polar(Mirtheta1RadFinal1, MirnormPower, color ='k', linewidth=2)
    plt.polar(th * pi / 180, r_line, color='r', linewidth=3)

    if (max(MirnormPower)<2):
        inc = 0.5
    elif (max(MirnormPower)<5):
        inc = 1
    elif max(MirnormPower)<20:
        inc = 5
    else:
        inc = 10
    plt.yticks(np.arange(inc, max(MirnormPower), inc), **ticksfont)
    plt.xticks(**ticksfont)
    angDist.savefig(dir+'angDist_' + number.__str__(), bbox_inches='tight')
    plt.close()
    return t, angDist

def trace_process(data, s11=True, plot=True):
    """Calculate the bandwidth from VNA trace data."""
    # Columns of data are: freq (Hz), mag (dB), phase (deg)
    if s11:
        peak_freq = np.argmin(data[:, 1])
    else:
        peak_freq = np.argmax(data[:, 1])

    phase_at_peak = data[peak_freq, 2]
    # subtract phase from all measurements and convert to radians
    data[:, 2] = (data[:, 2] - phase_at_peak)*np.pi*0.00555555555
    # convert dB S11 to lin S11
    data[:, 1] = 10 ** (data[:, 1] / 20)
    # convert to complex and then to imaginary
    imaginary_part = data[:, 1]*np.sin(data[:, 2])
    # find max and min imaginary
    freq_imag_min = data[np.argmin(imaginary_part), 0]
    freq_imag_max = data[np.argmax(imaginary_part), 0]

    # printing extras
    format_string = "Frequency of imaginary {} {:d}"
    flag = "minimum" if s11 else "maximum"
    # Data output
    print("Frequency at {}: {:d} Hz".format(flag, int(data[peak_freq, 0])))
    #print(format_string.format("maximum", int(freq_imag_max)))
    #print(format_string.format("minimum", int(freq_imag_min)))
    bandwidth = abs(freq_imag_max - freq_imag_min)
    print("Bandwidth {:.6f} kHz".format(bandwidth * 1e-3))
    q = data[peak_freq, 0] / bandwidth
    print("Q {:.0f}".format(q))
    if plot is False:
        return
    # Plots
    plt.figure(0, figsize=(8, 5))
    plt.plot(data[:, 0], imaginary_part)
    plt.title("Imaginary vs. frequency")
    plt.savefig("imag.png")
    plt.clf()
    plt.polar(data[:, 2], data[:, 1])  # first phase, then r
    plt.scatter(data[peak_freq, 2], data[peak_freq, 1], c='r', marker='o')
    plt.title("S11 in polar", va="bottom")
    plt.savefig("polar.png")
    plt.clf()
    plt.plot(data[:, 0], data[:, 1])
    plt.title("S11 curve in linear scale")
    plt.savefig("curve.png")
    plt.close("all")

 def PlotPolarDot(self,distance=False):
     """
     Makes a polar plot of the data, TODO: with an optional filter condition
     
     If distance is true, plots against comoving distance rather than redshift.
     """
     
     plt.clf()
     plt.title("Polar Plot of sample "+self.sample_name+" from "+self.survey_name)
     
     if not distance:
         plt.polar(self.survey_data['ra'],self.survey_data['redshift'],'.',markersize=1)
         plt.savefig(self.plots_dir+'/PolarPlot_Redshift.eps')
     else:
         plt.polar(self.survey_data['ra'],self.survey_data['c_dist'],'.',markersize=1)
         plt.savefig(self.plots_dir+'/PolarPlot_Distance.eps')
     plt.show()    

 def _map(self, param_x, param_y, data=None, method=np.mean, noplot=False, fig=None, ax=None, cax=None, bin_x=50, bin_y=50, cmap="jet", cm_min=None, cm_max=None, axescolor='w', polar=False, showmax=True, **kwargs):
     """
     Return a 2D histogram of the MC chain, showing the walker density per bin
     """
     if param_x not in self.paramstr or param_y not in self.paramstr:
         print("You must choose param_x and param_y among %s" % self.paramstr)
         return
     x = self.chain[param_x]
     if hasattr(x, 'compressed'): x = x.compressed()
     y = self.chain[param_y]
     if hasattr(y, 'compressed'): y = y.compressed()
     if data is not None:
         H, bin_y, bin_x = binned_statistic_2d(y, x, data, method, bins=(bin_y, bin_x))[:3]
     else:
         H, bin_y, bin_x = np.histogram2d(y, x, bins=(bin_y, bin_x))
     H[np.isnan(H)] = np.nanmin(H)
     maxd = np.unravel_index(np.argmax(H), H.shape)
     if cm_min is None: cm_min = np.min(H)
     if cm_max is None: cm_max = np.max(H)
     cmap, norm, mappable = _core.colorbar(cmap=cmap, cm_min=cm_min, cm_max=cm_max)
     if noplot: return H, bin_y, bin_x, cmap, norm, mappable
     if not polar:
         if fig is None: fig = plt.figure()
         if ax is None: ax = fig.add_subplot(111)
         im = mplimageNonUniformImage(ax, cmap=cmap, norm=norm, interpolation='bilinear')
         arrbin_x = _core.bins_to_array(bin_x)
         arrbin_y = _core.bins_to_array(bin_y)
         im.set_data(arrbin_x, arrbin_y, H)
         ax.images.append(im)
         if showmax is True: ax.plot(arrbin_x[maxd[1]], arrbin_y[maxd[0]], '^w', ms=7)
         ax.set_xlim(bin_x[0], bin_x[-1])
         ax.set_xlabel(param_x)
         ax.set_ylim(bin_y[0], bin_y[-1])
         ax.set_ylabel(param_y)
     else:
         if ax is None: fig, ax = plt.subplots(subplot_kw={'projection':'polar'})
         T,R = np.meshgrid(bin_x,bin_y)
         pax = ax.pcolormesh(T, R, H, cmap=cmap, norm=norm)
         ax.grid(True)
         if showmax is True: plt.polar(T[maxd], R[maxd], '^w', ms=7)
         ax.set_ylim(0, bin_y[-1])
         ax.set_title(param_x+' vs '+param_y)
     ax.grid(True, color=axescolor)
     ax.tick_params(axis='both', colors=axescolor)
     fig.colorbar(mappable, cax=cax)