Python BHTree.evolve_model方法代码示例

# 需要导入模块: from amuse.community.bhtree.interface import BHTree [as 别名]
# 或者: from amuse.community.bhtree.interface.BHTree import evolve_model [as 别名]
    def test19(self):
        particles = datamodel.Particles(2)
        particles.x = [0.0,10.0] | nbody_system.length
        particles.y = 0.0 | nbody_system.length
        particles.z = 0.0 | nbody_system.length
        particles.radius = 0.005 | nbody_system.length
        particles.vx =  0.0 | nbody_system.speed
        particles.vy =  0.0 | nbody_system.speed
        particles.vz =  0.0 | nbody_system.speed
        particles.mass = 1.0 | nbody_system.mass

        very_short_time_to_evolve = 1 | units.s
        very_long_time_to_evolve = 1e9 | nbody_system.time
       
        instance = BHTree()
        instance.initialize_code()
        instance.parameters.stopping_conditions_timeout = very_short_time_to_evolve 
        self.assertEquals(instance.parameters.stopping_conditions_timeout, very_short_time_to_evolve)
        instance.parameters.epsilon_squared = (0.01 | nbody_system.length)**2
        instance.particles.add_particles(particles) 
        instance.stopping_conditions.timeout_detection.enable()
        start = time.time()
        instance.evolve_model(very_long_time_to_evolve)
        end = time.time()
        self.assertTrue(instance.stopping_conditions.timeout_detection.is_set())
        self.assertTrue((end-start) < very_short_time_to_evolve.value_in(units.s) + 2)#2 = some overhead compensation
        instance.stop()

# 需要导入模块: from amuse.community.bhtree.interface import BHTree [as 别名]
# 或者: from amuse.community.bhtree.interface.BHTree import evolve_model [as 别名]
    def test17(self):
        print "Testing BHTree collision_detection"
        particles = datamodel.Particles(7)
        particles.mass = 0.001 | nbody_system.mass
        particles.radius = 0.01 | nbody_system.length
        particles.x = [-101.0, -100.0, -0.5, 0.5, 100.0, 101.0, 104.0] | nbody_system.length
        particles.y = 0 | nbody_system.length
        particles.z = 0 | nbody_system.length
        particles.velocity = [[2, 0, 0], [-2, 0, 0]]*3 + [[-4, 0, 0]] | nbody_system.speed
        
        instance = BHTree(redirection='none')
        instance.initialize_code()
        instance.parameters.set_defaults()
        
        # Uncommenting any of the following two lines will suppress collision detection
#~        instance.parameters.use_self_gravity = 0
#~        instance.parameters.epsilon_squared = 0.0 | nbody_system.length**2
        
        instance.parameters.opening_angle = 0.1
        instance.particles.add_particles(particles)
        collisions = instance.stopping_conditions.collision_detection
        collisions.enable()
        instance.evolve_model(1.0 | nbody_system.time)
        
        self.assertTrue(collisions.is_set())
        self.assertTrue(instance.model_time < 0.5 | nbody_system.time)
        self.assertEquals(len(collisions.particles(0)), 3)
        self.assertEquals(len(collisions.particles(1)), 3)
        self.assertEquals(len(particles - collisions.particles(0) - collisions.particles(1)), 1)
        self.assertEquals(abs(collisions.particles(0).x - collisions.particles(1).x) < 
                (collisions.particles(0).radius + collisions.particles(1).radius),
                [True, True, True])
        
        sticky_merged = datamodel.Particles(len(collisions.particles(0)))
        sticky_merged.mass = collisions.particles(0).mass + collisions.particles(1).mass
        sticky_merged.radius = collisions.particles(0).radius
        for p1, p2, merged in zip(collisions.particles(0), collisions.particles(1), sticky_merged):
            merged.position = (p1 + p2).center_of_mass()
            merged.velocity = (p1 + p2).center_of_mass_velocity()
        
        print instance.model_time
        print instance.particles
        instance.particles.remove_particles(collisions.particles(0) + collisions.particles(1))
        instance.particles.add_particles(sticky_merged)
        
        instance.evolve_model(1.0 | nbody_system.time)
        print
        print instance.model_time
        print instance.particles
        self.assertTrue(collisions.is_set())
        self.assertTrue(instance.model_time < 1.0 | nbody_system.time)
        self.assertEquals(len(collisions.particles(0)), 1)
        self.assertEquals(len(collisions.particles(1)), 1)
        self.assertEquals(len(instance.particles - collisions.particles(0) - collisions.particles(1)), 2)
        self.assertEquals(abs(collisions.particles(0).x - collisions.particles(1).x) < 
                (collisions.particles(0).radius + collisions.particles(1).radius),
                [True])
        instance.stop()

# 需要导入模块: from amuse.community.bhtree.interface import BHTree [as 别名]
# 或者: from amuse.community.bhtree.interface.BHTree import evolve_model [as 别名]
    def test2(self):
        #not completed 
        convert_nbody = nbody_system.nbody_to_si(1.0 | units.MSun, 149.5e6 | units.km)

        instance = BHTree(convert_nbody)
        #instance.dt_dia = 1
        instance.parameters.epsilon_squared = 0.001 | units.AU**2
        #instance.timestep = 0.0001
        #instance.use_self_gravity = 0
        instance.commit_parameters()
        
        stars = datamodel.Stars(2)
        sun = stars[0]
        sun.mass = units.MSun(1.0)
        sun.position = units.m(numpy.array((0.0,0.0,0.0)))
        sun.velocity = units.ms(numpy.array((0.0,0.0,0.0)))
        sun.radius = units.RSun(1.0)

        earth = stars[1]
        earth.mass = units.kg(5.9736e24)
        earth.radius = units.km(6371) 
        earth.position = units.km(numpy.array((149.5e6,0.0,0.0)))
        earth.velocity = units.ms(numpy.array((0.0,29800,0.0)))

        instance.particles.add_particles(stars)
        instance.commit_particles()
        self.assertAlmostRelativeEquals(sun.radius, instance.particles[0].radius)
    
        for x in range(1,2000,10):
            instance.evolve_model(x | units.day)
            instance.particles.copy_values_of_all_attributes_to(stars)
            stars.savepoint()
            
        if HAS_MATPLOTLIB:
            figure = pyplot.figure()
            plot = figure.add_subplot(1,1,1)
            
            x_points = earth.get_timeline_of_attribute("x")
            y_points = earth.get_timeline_of_attribute("y")
            
            x_points_in_AU = map(lambda (t,x) : x.value_in(units.AU), x_points)
            y_points_in_AU = map(lambda (t,x) : x.value_in(units.AU), y_points)
            
            plot.scatter(x_points_in_AU,y_points_in_AU, color = "b", marker = 'o')
            
            plot.set_xlim(-1.5, 1.5)
            plot.set_ylim(-1.5, 1.5)
               
            
            test_results_path = self.get_path_to_results()
            output_file = os.path.join(test_results_path, "bhtree-earth-sun.svg")
            figure.savefig(output_file)    
        
        instance.cleanup_code()
        instance.stop()

# 需要导入模块: from amuse.community.bhtree.interface import BHTree [as 别名]
# 或者: from amuse.community.bhtree.interface.BHTree import evolve_model [as 别名]
def nbody_integrator(Ncl, mcl, rcl, t_end, n_steps, escape_velocity_fraction, R):
    converter = nbody_system.nbody_to_si(mcl, rcl)
    bodies = new_plummer_model(Ncl, convert_nbody=converter)

    #estimate of milky way mass by "Mass models of the Milky Way", McMillan
    blackhole_mass = 1.26e12 | units.MSun
    blackhole = Particle(mass=blackhole_mass)
    blackhole.position = [0,0,0] | units.m

    cluster_velocity = [0,0,0] | units.m / units.s
    cluster_position = [0,0,0] | units.parsec
    cluster_position[0] = R
    G_si = converter.to_si(nbody_system.G)
    escape_v = (2*G_si*blackhole_mass/R).sqrt().as_quantity_in(units.m/units.s)
    V = escape_v * escape_velocity_fraction
    cluster_velocity[1] = V
    bodies.move_to_center()
    bodies.velocity += cluster_velocity
    bodies.position += cluster_position
    bodies.add_particle(blackhole)
        

    gravity = BHTree(converter)
    gravity.particles.add_particles(bodies)
    channel_from_gravity_to_framework = gravity.particles.\
        new_channel_to(bodies)

    time = zero
    dt = t_end / float(n_steps)
    x = 0
    base_path = "encounter_plots/"+str(R.value_in(units.parsec))+"_"+\
        str(escape_velocity_fraction)+"_"
    while time < t_end:
        plot_cluster(bodies, base_path+str(x),time, rcl, V)
        time += dt
        gravity.evolve_model(time)
        channel_from_gravity_to_framework.copy()
        x+=1
    plot_cluster(bodies, base_path+str(x),time, rcl, V)

    gravity.stop()
    return V, bodies

# 需要导入模块: from amuse.community.bhtree.interface import BHTree [as 别名]
# 或者: from amuse.community.bhtree.interface.BHTree import evolve_model [as 别名]
 def test18(self):
     particles = datamodel.Particles(2)
     particles.x = [0.0,10.0] | nbody_system.length
     particles.y = 0 | nbody_system.length
     particles.z = 0 | nbody_system.length
     particles.radius = 0.005 | nbody_system.length
     particles.vx =  0 | nbody_system.speed
     particles.vy =  0 | nbody_system.speed
     particles.vz =  0 | nbody_system.speed
     particles.mass = 1.0 | nbody_system.mass
    
     instance = BHTree()
     instance.initialize_code()
     instance.parameters.stopping_conditions_number_of_steps = 2
     self.assertEquals(instance.parameters.stopping_conditions_number_of_steps, 2)
     instance.parameters.epsilon_squared = (0.01 | nbody_system.length)**2
     instance.particles.add_particles(particles) 
     instance.stopping_conditions.number_of_steps_detection.enable()
     instance.evolve_model(10 | nbody_system.time)
     self.assertTrue(instance.stopping_conditions.number_of_steps_detection.is_set())
     self.assertTrue(instance.model_time < 10 | nbody_system.time)
     instance.stop()

# 需要导入模块: from amuse.community.bhtree.interface import BHTree [as 别名]
# 或者: from amuse.community.bhtree.interface.BHTree import evolve_model [as 别名]
    def test3(self):
        convert_nbody = nbody_system.nbody_to_si(1.0 | units.MSun, 149.5e6 | units.km)

        instance = BHTree(convert_nbody)
        #instance.dt_dia = 1
        instance.parameters.epsilon_squared = 0.001 | units.AU**2
        #instance.timestep = 0.0001
        #instance.use_self_gravity = 0
        instance.commit_parameters()
        
        
        stars = datamodel.Stars(2)
        star1 = stars[0]
        star2 = stars[1]

        star1.mass = units.MSun(1.0)
        star1.position = units.AU(numpy.array((-.10,0.0,0.0)))
        star1.velocity = units.AUd(numpy.array((0.0,0.0,0.0)))
        star1.radius = units.RSun(1.0)

        star2.mass = units.MSun(1.0)
        star2.position = units.AU(numpy.array((.10,0.0,0.0)))
        star2.velocity = units.AUd(numpy.array((0.0,0.0,0.0)))
        star2.radius = units.RSun(100.0)
        
        instance.particles.add_particles(stars)
        instance.commit_particles()
    
        for x in range(1,200,1):
            instance.evolve_model(x | units.day)
            instance.particles.copy_values_of_all_attributes_to(stars)
            #instance.get_indices_of_colliding_particles()
            #print stars[0].position-stars[1].position
            stars.savepoint()
            
        instance.cleanup_code()
        instance.stop()

# 需要导入模块: from amuse.community.bhtree.interface import BHTree [as 别名]
# 或者: from amuse.community.bhtree.interface.BHTree import evolve_model [as 别名]
    def test23(self):
        particles = datamodel.Particles(2)
        particles.x = [0.0,10.0] | nbody_system.length
        particles.y = 0.0 | nbody_system.length
        particles.z = 0.0 | nbody_system.length
        particles.vx =  1.0 | nbody_system.speed
        particles.vy =  0.0 | nbody_system.speed
        particles.vz =  0.0 | nbody_system.speed
        particles.mass = 0.1 | nbody_system.mass

        instance = BHTree(redirection="none")
        instance.particles.add_particles(particles) 
        instance.commit_particles()
        instance.evolve_model(0.1 | nbody_system.time)
        self.assertFalse(instance.particles[0].vy > 0| nbody_system.speed)
        self.assertAlmostRelativeEquals(instance.particles[0].x , 0.1 | nbody_system.length, 4)
        instance.particles.new_channel_to(particles).copy()
        particles.vy = 1| nbody_system.speed
        particles.new_channel_to(instance.particles).copy()
        
        instance.evolve_model(0.2 | nbody_system.time)
        self.assertTrue(instance.particles[0].vy > 0| nbody_system.speed)
        self.assertAlmostRelativeEquals(instance.particles[0].y , 0.1 | nbody_system.length, 4)
        instance.stop()

# 需要导入模块: from amuse.community.bhtree.interface import BHTree [as 别名]
# 或者: from amuse.community.bhtree.interface.BHTree import evolve_model [as 别名]
    def test15(self):
        print "Test15: Testing effect of BHTree parameter epsilon_squared"
        convert_nbody = nbody_system.nbody_to_si(1.0 | units.MSun, 1.0 | units.AU)
        
        particles = datamodel.Particles(2)
        sun = particles[0]
        sun.mass = 1.0 | units.MSun
        sun.position = [0.0, 0.0, 0.0] | units.AU
        sun.velocity = [0.0, 0.0, 0.0] | units.AU / units.yr
        sun.radius = 1.0 | units.RSun

        earth = particles[1]
        earth.mass = 5.9736e24 | units.kg
        earth.radius = 6371.0 | units.km
        earth.position = [0.0, 1.0, 0.0] | units.AU
        earth.velocity = [2.0*numpy.pi, -0.0001, 0.0] | units.AU / units.yr
        
        initial_direction = math.atan((earth.velocity[0]/earth.velocity[1]))
        final_direction = []
        for log_eps2 in range(-9,10,2):
            instance = BHTree(convert_nbody)
            instance.initialize_code()
            instance.parameters.epsilon_squared = 10.0**log_eps2 | units.AU ** 2
            instance.particles.add_particles(particles)
            instance.commit_particles()
            instance.evolve_model(0.25 | units.yr)
            final_direction.append(math.atan((instance.particles[1].velocity[0]/
                instance.particles[1].velocity[1])))
            instance.stop()
        # Small values of epsilon_squared should result in normal earth-sun dynamics: rotation of 90 degrees
        self.assertAlmostEquals(abs(final_direction[0]), abs(initial_direction+math.pi/2.0), 2)
        # Large values of epsilon_squared should result in ~ no interaction
        self.assertAlmostEquals(final_direction[-1], initial_direction, 2)
        # Outcome is most sensitive to epsilon_squared when epsilon_squared = d(earth, sun)^2
        delta = [abs(final_direction[i+1]-final_direction[i]) for i in range(len(final_direction)-1)]
        self.assertEquals(delta[len(final_direction)//2 -1], max(delta))

# 需要导入模块: from amuse.community.bhtree.interface import BHTree [as 别名]
# 或者: from amuse.community.bhtree.interface.BHTree import evolve_model [as 别名]
    def test1(self):
        convert_nbody = nbody_system.nbody_to_si(1.0 | units.MSun, 149.5e6 | units.km)

        instance = BHTree(convert_nbody)
        instance.parameters.epsilon_squared = 0.001 | units.AU**2
        
        stars = datamodel.Stars(2)
        
        sun = stars[0]
        sun.mass = units.MSun(1.0)
        sun.position = [0.0,0.0,0.0] | units.m
        sun.velocity = [0.0,0.0,0.0] | units.ms
        sun.radius = units.RSun(1.0)

        earth = stars[1]
        earth.mass = units.kg(5.9736e24)
        earth.radius = units.km(6371) 
        earth.position = [149.5e6, 0.0, 0.0] | units.km
        earth.velocity = [0.0, 29800, 0.0] | units.ms

        #instance.particles.add_particles(stars)
        instance.particles.add_particles(stars)
        
        postion_at_start = earth.position.value_in(units.AU)[0]
        
        instance.evolve_model(365.0 | units.day)
        instance.particles.copy_values_of_all_attributes_to(stars)
        
        postion_after_full_rotation = earth.position.value_in(units.AU)[0]
       
        self.assertAlmostEqual(postion_at_start, postion_after_full_rotation, 3)
        
        instance.evolve_model(365.0 + (365.0 / 2) | units.day)
        
        instance.particles.copy_values_of_all_attributes_to(stars)
        
        postion_after_half_a_rotation = earth.position.value_in(units.AU)[0]
        self.assertAlmostEqual(-postion_at_start, postion_after_half_a_rotation, 2)
        
        
        instance.evolve_model(365.0 + (365.0 / 2) + (365.0 / 4)  | units.day)
         
        instance.particles.copy_values_of_all_attributes_to(stars)
        
        postion_after_half_a_rotation = earth.position.value_in(units.AU)[1]
        
        self.assertAlmostEqual(-postion_at_start, postion_after_half_a_rotation, 1)
        instance.cleanup_code()
        instance.stop()

# 需要导入模块: from amuse.community.bhtree.interface import BHTree [as 别名]
# 或者: from amuse.community.bhtree.interface.BHTree import evolve_model [as 别名]
def assignment_2d():
    current_cluster_mass = 400 | units.MSun
    initial_mass_fraction = 0.84
    desired_initial_mass = current_cluster_mass / initial_mass_fraction

    masses = new_salpeter_mass_distribution(100000)
    mean_salpeter_mass = masses.mean()
    print "mean salpeter mass", mean_salpeter_mass
    N = int(desired_initial_mass / mean_salpeter_mass)
    print "N", N


    Rvir = 10 | units.lightyear
    z = 0.17
    masses = new_salpeter_mass_distribution(N)
    converter = nbody_system.nbody_to_si(masses.sum(), Rvir)
    G_SI = converter.to_si(nbody_system.G)
    bodies = new_plummer_sphere(N, convert_nbody=converter)
    bodies.mass = masses
    bodies.metalicity = z

    # start the gravity solver
    gravity = BHTree(converter)
    gravity.initialize_code()
    gravity.parameters.timestep = 0.1 | units.Myr

    # start the stellar evolution solver
    stellar = SSE()
    stars = stellar.particles.add_particles(bodies)
    from_stellar_evolution_to_model \
        = stellar.particles.new_channel_to(bodies)
    from_stellar_evolution_to_model.copy_attributes(["mass"])

    bodies.scale_to_standard(converter)
    gravity.particles.add_particles(bodies)

    from_model_to_gravity = bodies.new_channel_to(gravity.particles)
    from_gravity_to_model = gravity.particles.new_channel_to(bodies)
    gravity.commit_particles()

    end_time = 1000 | units.Myr
    current_time = 0 | units.Myr
    cluster = "Hyades"
    bound_stars_counts = []
    main_sequence_stars_counts = []
    giant_stars_counts = []
    remnant_stars_counts = []
    max_radii = [] | units.parsec
    virial_radii = [] | units.parsec
    times = [] | units.Myr
    while current_time < end_time:
        name_of_the_figure = "isochrone_with_grav_"+str(int(current_time.value_in(units.Myr)))+".png"

        gravity.evolve_model(current_time)
        stellar.evolve_model(current_time)


        from_gravity_to_model.copy()
        from_stellar_evolution_to_model.copy_attributes(["mass", "radius"])
        from_model_to_gravity.copy_attributes(["mass"])


        remnant_count = 0
        main_sequence_count = 0
        giant_count = 0
        for star in stars:
            if stellar_remnant_state(star):
                remnant_count += 1
            if stellar_giant_state(star):
                giant_count += 1
            if stellar_main_sequence_state(star):
                main_sequence_count += 1
        max_radius = bodies.total_radius()
        virial_radius = bodies.virial_radius()
        bound_star_count = len(bodies.bound_subset(unit_converter=converter, G=G_SI))
        print "bound stars:", bound_star_count
        print "main sequence stars:", main_sequence_count
        print "giant stars:", giant_count
        print "remnant stars:", remnant_count
        print "cluster radius(max from centre):", max_radius
        print "virial radius:", virial_radius
        print current_time

        times.append(current_time)
        remnant_stars_counts.append(remnant_count)
        giant_stars_counts.append(giant_count)
        main_sequence_stars_counts.append(main_sequence_count)
        max_radii.append(max_radius)
        virial_radii.append(virial_radius)
        bound_stars_counts.append(bound_star_count)


        temperatures = stars.temperature
        luminosities = stars.luminosity

        plot_HR_diagram(temperatures, luminosities,
                        cluster+"/",
                        name_of_the_figure, current_time)
        current_time += 10 | units.Myr
    data = {}
#.........这里部分代码省略.........

# 需要导入模块: from amuse.community.bhtree.interface import BHTree [as 别名]
# 或者: from amuse.community.bhtree.interface.BHTree import evolve_model [as 别名]
def simulate_small_cluster(number_of_stars, end_time = 40 | units.Myr,
                           name_of_the_figure = "test-2.svg"):
    #numpy.random.seed(1)
    
    salpeter_masses = new_salpeter_mass_distribution(number_of_stars)
    total_mass = salpeter_masses.sum()
    
    convert_nbody = nbody_system.nbody_to_si(total_mass, 1.0 | units.parsec)
    
    particles = new_plummer_model(number_of_stars, convert_nbody);
   
    gravity = BHTree(convert_nbody)
    gravity.initialize_code()
    #gravity.parameters.set_defaults()
    #print gravity.parameters.timestep.as_quantity_in(units.Myr)
    gravity.parameters.timestep = 0.0001 | units.Myr	# tiny!
    gravity.parameters.epsilon_squared \
        = (float(number_of_stars)**(-0.333333) | units.parsec) ** 2
        
    stellar_evolution = SSE()
    stellar_evolution.initialize_module_with_default_parameters() 
    
    print "setting masses of the stars"
    particles.radius = 0.0 | units.RSun
    particles.mass = salpeter_masses
    
    print "initializing the particles"
    stellar_evolution.particles.add_particles(particles)
    from_stellar_evolution_to_model \
        = stellar_evolution.particles.new_channel_to(particles)
    from_stellar_evolution_to_model.copy_attributes(["mass"])
    
    print "centering the particles"
    particles.move_to_center()
    print "scaling particles to viridial equilibrium"
    particles.scale_to_standard(convert_nbody)
    
    gravity.particles.add_particles(particles)
    from_model_to_gravity = particles.new_channel_to(gravity.particles)
    from_gravity_to_model = gravity.particles.new_channel_to(particles)
    
    gravity.commit_particles()
        
    time = 0.0 | units.Myr    
    particles.savepoint(time)
    
    
    total_energy_at_t0 = gravity.kinetic_energy + gravity.potential_energy   
    
    print "evolving the model until t = " + str(end_time)
    while time < end_time:
        time += 0.25 | units.Myr
        
        print "gravity evolve step starting"
        gravity.evolve_model(time)
        print "gravity evolve step done"
        
        print "stellar evolution step starting"
        stellar_evolution.evolve_model(time)
        print "stellar evolution step done"

        from_gravity_to_model.copy()
        from_stellar_evolution_to_model.copy_attributes(["mass", "radius"])

        particles.savepoint(time)  
        
        from_model_to_gravity.copy_attributes(["mass"])
        
        total_energy_at_this_time \
            = gravity.kinetic_energy + gravity.potential_energy   
        print_log(time, gravity, particles,
                  total_energy_at_t0, total_energy_at_this_time)

    
    test_results_path = get_path_to_results()
    output_file = os.path.join(test_results_path, "small.hdf5")
    if os.path.exists(output_file):
        os.remove(output_file)
    storage = store.StoreHDF(output_file)
    storage.store(particles)
   
    gravity.stop()
    stellar_evolution.stop()
    
    plot_particles(particles, name_of_the_figure)