-
Notifications
You must be signed in to change notification settings - Fork 115
Expand file tree
/
Copy pathclump_finding.py
More file actions
102 lines (76 loc) · 2.67 KB
/
Copy pathclump_finding.py
File metadata and controls
102 lines (76 loc) · 2.67 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
# -*- coding: ascii -*-
"""
Generates a random distribution of particles and uses Hop to determine clumps
in the particle set.
"""
from __future__ import print_function
# import numpy
from numpy import random
from matplotlib import pyplot
from amuse.units import units
from amuse.units import nbody_system
from amuse.community.hop.interface import Hop
from amuse.datamodel.particles import Particles
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
)
particles = Particles(number_of_stars)
particles.mass = masses
particles.x = units.parsec(random.gamma(2.0, 1.0, number_of_stars))
particles.y = units.parsec(random.gamma(1.0, 1.0, number_of_stars))
particles.z = units.parsec(random.random(number_of_stars))
return particles
def find_clumps(particles, unit_converter):
hop = Hop(unit_converter)
hop.particles.add_particles(particles)
hop.calculate_densities()
hop.do_hop()
result = [x.get_intersecting_subset_in(particles) for x in hop.groups()]
hop.stop()
return result
def plot_clumps(groups, total_mass):
number_of_particles_in_group = []
fraction_of_mass_in_group = []
for group in groups:
number_of_particles_in_group.append(len(group))
fraction = (group.mass.sum() / total_mass)
fraction_of_mass_in_group.append(fraction)
figure = pyplot.figure(figsize=(12, 6))
subplot = figure.add_subplot(1, 2, 1)
colormap = pyplot.cm.Paired
for index, group in enumerate(groups):
color = colormap(1.0 * index / len(groups))
subplot.scatter(
group.x.value_in(units.parsec),
group.y.value_in(units.parsec),
s=group.mass.value_in(units.MSun),
edgecolors=color,
facecolors=color
)
subplot.set_xlim(0, 1)
subplot.set_ylim(0, 1)
subplot.set_xlabel('x (parsec)')
subplot.set_ylabel('y (parsec)')
subplot = figure.add_subplot(1, 2, 2)
subplot.plot(
number_of_particles_in_group,
fraction_of_mass_in_group,
)
subplot.set_xscale('log')
subplot.set_yscale('log')
subplot.set_xlabel('N')
subplot.set_ylabel('df/d(Log_10 N)')
figure.savefig('x.png')
pyplot.show()
if __name__ == "__main__":
number_of_stars = 10000
stars = new_cluster(number_of_stars)
total_mass = stars.mass.sum()
unit_converter = nbody_system.nbody_to_si(total_mass, 1 | units.parsec)
groups = find_clumps(stars, unit_converter)
plot_clumps(groups, total_mass)