code: """
  Generates a cluster using a plummer model with a salpeter Initial Mass Function.
  Compares the generated IMF against the expected line.
  """
  
  import numpy 
  from matplotlib import pyplot
  from amuse.units import units
  from amuse.units import nbody_system
  from amuse.ic.plummer import new_plummer_model
  from amuse.ic.salpeter import new_salpeter_mass_distribution
  
  def new_cluster(number_of_stars = 1000):
      masses = new_salpeter_mass_distribution(
          number_of_stars, 
          mass_min = 0.1 | units.MSun,
          mass_max = 125.0 | units.MSun, 
          alpha = -2.35
      )
      nbody_converter = nbody_system.nbody_to_si(masses.sum(), 1 | units.parsec)
      particles = new_plummer_model(number_of_stars, nbody_converter)
      particles.mass = masses
      particles.move_to_center()
      return particles
  
  def plot_particles_and_mass_distribution(particles):
      figure = pyplot.figure(figsize= (12,6))
      
      subplot = figure.add_subplot(1, 2, 1)
          
      subplot.scatter(
          particles.x.value_in(units.parsec),
          particles.y.value_in(units.parsec),
          s = particles.mass.value_in(units.MSun),# * len(particles),
          edgecolors = 'red',
          facecolors = 'red'
      )
      
      subplot.set_xlim(-4,4)
      subplot.set_ylim(-4,4)
      subplot.set_xlabel('x (parsec)')
      subplot.set_ylabel('y (parsec)')
      
      subplot = figure.add_subplot(1, 2, 2)
      
      masses = particles.mass.value_in(units.MSun)
      
      bins = 10**numpy.linspace(-1, 2, 100)
      number_of_particles, bin_edges= numpy.histogram(masses, bins = bins)
      
      bin_sizes = bin_edges[1:] - bin_edges[:-1]
      
      y = number_of_particles / bin_sizes
      x = (bin_edges[1:] + bin_edges[:-1]) / 2.0
      
      y = y[number_of_particles > 10.0] 
      x = x[number_of_particles > 10.0]
      subplot.scatter(x, y)
      
      c = ((0.1**-1.35) - (125.0**-1.35)) / 1.35
      subplot.plot(x, len(particles)/ c * (x**-2.35))
      
      subplot.set_xscale('log')
      subplot.set_yscale('log')
      
      subplot.set_xlabel(u'M [M\u2299]')
      subplot.set_ylabel('N')
      
      pyplot.show()
      
  if __name__ == "__main__":
      particles = new_cluster(20000)
      plot_particles_and_mass_distribution(particles)
   """
  Evolves two stars dynamically (hermit, nbody code) each star will 
  lose mass during the evolution (evtwin, stellar evolution code)
  
  We start with two stars, one 10.0 and one 1.0 solar mass star. These
  stars start orbiting with a stable kepler orbit. 
  After 2 orbital periods the stars will begin to lose mass and the binary
  will become unstable.
  """
  
  from amuse.plot import scatter, xlabel, ylabel, plot
  from matplotlib import pyplot
  from math import pi
  from amuse.units.optparse import OptionParser
  
  from amuse.units import units
  from amuse.units import constants
  from amuse.units.nbody_system import nbody_to_si
  from amuse.community.evtwin.interface import EVtwin
  from amuse.community.sse.interface import SSE
  from amuse.community.hermite0.interface import Hermite
  
  from amuse.datamodel import Particles
  
  
  def set_up_initial_conditions(orbital_period, kinetic_to_potential_ratio):
      print("Setting up initial conditions")
      stars =  Particles(2)
      stars.mass = [10.0, 1.0] | units.MSun
      stars.radius = 0 | units.RSun
      stars.position = [0.0, 0.0, 0.0] | units.AU
      stars.velocity = [0.0, 0.0, 0.0] | units.km / units.s
      
      print("Binary with masses: "+str(stars.mass)+", and orbital period: ", orbital_period)
      semimajor_axis = ((constants.G * stars.total_mass() * (orbital_period / (2 * pi))**2.0)**(1.0/3.0))
      separation = 2 * semimajor_axis * (1 - kinetic_to_potential_ratio)
      print("Initial separation:", separation.as_quantity_in(units.AU))
      relative_velocity = ( (kinetic_to_potential_ratio / (1.0 - kinetic_to_potential_ratio)) * 
          constants.G * stars.total_mass() / semimajor_axis).sqrt()
      print("Initial relative velocity:", relative_velocity.as_quantity_in(units.km / units.s))
      
      stars[0].x = separation
      stars[0].vy = relative_velocity
      stars.move_to_center()
      return stars
  
  def set_up_stellar_evolution_code(stars):
      stellar_evolution = EVtwin()
      stellar_evolution.initialize_code()
      # if you run with mesa, you can play with the wind efficiency
      # stellar_evolution.parameters.RGB_wind_scheme = 1
      # stellar_evolution.parameters.reimers_wind_efficiency = 1.0e6 # ridiculous, but instructive
      stellar_evolution.particles.add_particles(stars)
      return stellar_evolution
      
  def set_up_gravitational_dynamics_code(stars):
      convert_nbody = nbody_to_si(11.0 | units.MSun, 10.0 | units.AU)
      gravitational_dynamics = Hermite(convert_nbody)
      gravitational_dynamics.parameters.epsilon_squared = 0.0 | units.AU ** 2
      view_on_the_primary = gravitational_dynamics.particles.add_particle(stars[0])
      gravitational_dynamics.particles.add_particle(stars[1])
      return gravitational_dynamics, view_on_the_primary
      
  
  def simulate_binary_evolution(binary, orbital_period, t_offset_stars, t_end):
      distance = [] | units.AU
      mass = [] | units.MSun
      time = [] | units.yr
      
      stellar_evolution = set_up_stellar_evolution_code(binary)
      gravitational_dynamics, primary = set_up_gravitational_dynamics_code(binary)
      from_se_to_gd = stellar_evolution.particles.new_channel_to(gravitational_dynamics.particles)
      
      current_time = 0.0 * t_end
      
          
      print("Evolving with stellar wind")
      while current_time < t_end:
          current_time += orbital_period / 10
          gravitational_dynamics.evolve_model(current_time)
          stellar_evolution.evolve_model(current_time + t_offset_stars)
          from_se_to_gd.copy_attributes(['mass'])
          separation = (gravitational_dynamics.particles[0].position - gravitational_dynamics.particles[1].position).length()
          distance.append(separation)
          mass.append(primary.mass)
          time.append(current_time)
          print("System evolved to time: ", current_time, ", primary mass:", primary.mass.as_quantity_in(units.MSun), ", separation:", separation.as_quantity_in(units.AU))
      
      print("Evolution done")
      return distance, mass, time
  
  def orbit_plot(distance, mass, time):
      figure = pyplot.figure(figsize = (6, 10), dpi = 100)
      subplot = figure.add_subplot(2, 1, 1)
      plot(time, distance)
      xlabel('t')
      ylabel('separation')
      pyplot.margins(0.05)
      subplot = figure.add_subplot(2, 1, 2)
      plot(time, ((mass - mass[0]) / mass[0]) * 100.0)
      xlabel('t')
      ylabel('mass')
      pyplot.margins(0.05)
      pyplot.show()
  
  def main(
          orbital_period = 1000.0 | units.yr, 
          kinetic_to_potential_ratio = 0.8, 
          periods = 10,
          age = 10 | units.Myr
      ):
      
      t_offset_stars = age
      t_end = periods * orbital_period
      
      binary = set_up_initial_conditions(orbital_period, kinetic_to_potential_ratio)
      distance, mass, time = simulate_binary_evolution(binary, orbital_period, t_offset_stars, t_end)
      orbit_plot(distance, mass, time)
      
  def new_option_parser():
      result = OptionParser()
      result.add_option(
          "-o", "--orbitalperiod", 
          default = 1000 | units.yr,
          dest="orbital_period",
          help="initial orbital period of the binary (in years)",
          type="float",
          unit=units.yr
      )
      
      result.add_option(
          "-k", "--kpratio", 
          default = 0.8,
          dest="kinetic_to_potential_ratio",
          help="kinetec to potential energy ratio, values less than 1.0 correspond to bound systems",
          type="float"
      )
      result.add_option(
          "--periods", 
          default = 10,
          dest="periods",
          help="number of orbital periods to evolve the binary",
          type="int"
      )
      result.add_option(
          "--age", 
          default = 10 | units.Myr,
          dest="age",
          help="initial age of the stars to start the simulation with",
          type="float",
          unit=units.Myr
      )
      
      
      return result
      
      
  if __name__ == "__plot__":
      main(1000 | units.yr, 0.8, 10, 10 | units.Myr)
      
  if __name__ == "__main__":
      options, args = new_option_parser().parse_args()
      main(**options.__dict__)
  
apis: numpy.histogram matplotlib.pyplot.figure matplotlib.pyplot.show numpy.linspace
apis: matplotlib.pyplot.figure matplotlib.pyplot.show matplotlib.pyplot.margins