Skip to content

Commit 4a5500b

Browse files
LourensVeenLourensVeen
committed
Fix formatting
I don't have my Fortran safety goggles with me, so this will avoid eye injury and OSHA violations.
1 parent 230ca21 commit 4a5500b

1 file changed

Lines changed: 117 additions & 105 deletions

File tree

‎src/amuse/ext/protodisk.py‎

Lines changed: 117 additions & 105 deletions
Original file line numberDiff line numberDiff line change
@@ -8,48 +8,54 @@
88

99
from amuse.datamodel import Particles
1010
from amuse.datamodel import ParticlesWithUnitsConverted
11+
12+
1113
def approximate_inverse_error_function(x):
12-
a=8*(numpy.pi-3)/3*numpy.pi*(4-numpy.pi)
13-
return numpy.sign(x)*numpy.sqrt(
14-
numpy.sqrt((2/numpy.pi/a+numpy.log(1-x**2)/2)**2-numpy.log(1-x**2)/a)-(2/numpy.pi/a+numpy.log(1-x**2)/2)
15-
)
16-
17-
class uniform_unit_cylinder(object):
18-
def __init__(self,targetN, base_grid=None):
19-
cube_cylinder_ratio=numpy.pi*0.5**2
20-
self.targetN=targetN
21-
self.estimatedN=targetN/cube_cylinder_ratio
14+
a = 8 * (numpy.pi - 3) / 3 * numpy.pi * (4 - numpy.pi)
15+
16+
return numpy.sign(x) * numpy.sqrt(
17+
numpy.sqrt(
18+
(2 / numpy.pi / a + numpy.log(1 - x**2) / 2)**2 - numpy.log(1 - x**2) / a
19+
) - (2 / numpy.pi / a + numpy.log(1 - x**2) / 2)
20+
)
21+
22+
class uniform_unit_cylinder:
23+
def __init__(self, targetN, base_grid=None):
24+
cube_cylinder_ratio = numpy.pi*0.5**2
25+
self.targetN = targetN
26+
self.estimatedN = targetN / cube_cylinder_ratio
27+
2228
if base_grid is None:
23-
self.base_grid=uniform_random_unit_cube
29+
self.base_grid = uniform_random_unit_cube
2430
else:
25-
self.base_grid=base_grid
26-
27-
def cutout_cylinder(self,x,y,z):
28-
r=x**2+y**2
29-
selection=r < numpy.ones_like(r)
30-
x=x.compress(selection)
31-
y=y.compress(selection)
32-
z=z.compress(selection)
33-
return x,y,z
31+
self.base_grid = base_grid
32+
33+
def cutout_cylinder(self, x, y, z):
34+
r = x**2 + y**2
35+
selection = r < numpy.ones_like(r)
36+
x = x.compress(selection)
37+
y = y.compress(selection)
38+
z = z.compress(selection)
39+
return x, y, z
3440

3541
def make_xyz(self):
36-
if(self.base_grid==uniform_random_unit_cube):
37-
estimatedN=self.estimatedN
38-
x=[]
42+
if (self.base_grid == uniform_random_unit_cube):
43+
estimatedN = self.estimatedN
44+
x = []
3945
while len(x) < self.targetN:
40-
estimadedN=estimatedN*1.1+1
41-
x,y,z=self.cutout_cylinder(*(self.base_grid(estimatedN)).make_xyz())
42-
return x[0:self.targetN],y[0:self.targetN],z[0:self.targetN]
46+
estimadedN = estimatedN * 1.1 + 1
47+
x, y, z = self.cutout_cylinder(*(self.base_grid(estimatedN)).make_xyz())
48+
49+
return x[0:self.targetN], y[0:self.targetN], z[0:self.targetN]
4350
else:
4451
return self.cutout_cylinder(*(self.base_grid(self.estimatedN)).make_xyz())
4552

4653

4754
class ProtoPlanetaryDisk:
48-
4955
def __init__(
5056
self, targetN, convert_nbody=None, discfraction=0.1,
5157
densitypower=1., thermalpower=0.5, radius_min=1, radius_max=100,
52-
gamma=1.,q_out=2.,base_grid=None, Rmin=None, Rmax=None,
58+
gamma=1., q_out=2., base_grid=None, Rmin=None, Rmax=None,
5359
):
5460
if Rmin is not None:
5561
warnings.warn(
@@ -62,8 +68,10 @@ def __init__(
6268
"this is only allowed if one of them is None"
6369
)
6470
radius_min = Rmin
71+
6572
if radius_min is None:
6673
raise ValueError("radius_min must be set")
74+
6775
if Rmax is not None:
6876
warnings.warn(
6977
"Rmax is deprecated, use radius_max instead",
@@ -75,90 +83,94 @@ def __init__(
7583
"this is only allowed if one of them is None"
7684
)
7785
radius_max = Rmax
86+
7887
if radius_max is None:
7988
raise ValueError("radius_max must be set")
8089

81-
self.targetN=targetN
82-
self.convert_nbody=convert_nbody
83-
self.densitypower=densitypower
84-
self.thermalpower=thermalpower
85-
self.Rmin=radius_min
86-
self.Rmax=radius_max
87-
self.gamma=gamma
88-
self.q_out=q_out
89-
self.discfraction=discfraction
90-
91-
self.a=self.thermalpower
92-
self.a2=self.thermalpower/2
93-
self.g=densitypower
94-
self.g2=2-densitypower
95-
self.k_out=((1+discfraction)/Rmax**3)**0.5
96-
self.sigma_out=self.g2*discfraction/(2*numpy.pi*Rmax**self.g*(Rmax**self.g2-Rmin**self.g2))
97-
self.cs_out=self.q_out*numpy.pi*self.sigma_out/self.k_out
98-
99-
self.base_cylinder=uniform_unit_cylinder(targetN,base_grid)
100-
101-
102-
def sigma(self,r):
103-
return self.sigma_out*(self.Rmax/r)**self.g
104-
105-
def csound(self,r):
106-
return self.cs_out*(self.Rmax/r)**self.a2
107-
108-
def cmass(self,r):
109-
return self.discfraction*(r**self.g2-self.Rmin**self.g2)/(self.Rmax**self.g2-self.Rmin**self.g2)
110-
111-
def mass_encl(self,r):
112-
return 1+self.cmass(r)
113-
114-
def kappa(self,r):
115-
return (self.mass_encl(r)/r**3)**0.5
116-
117-
def toomreQ(self,r):
118-
return self.csound(r)*self.kappa(r)/numpy.pi/self.sigma(r)
119-
120-
def getradius(self,f):
121-
return ((self.Rmax**self.g2-self.Rmin**self.g2)*f+self.Rmin**self.g2)**(1./self.g2)
122-
123-
def zscale(self,r):
124-
return self.csound(r)/self.kappa(r)
125-
126-
def u(self,r):
127-
if self.gamma ==1.:
90+
self.targetN = targetN
91+
self.convert_nbody = convert_nbody
92+
self.densitypower = densitypower
93+
self.thermalpower = thermalpower
94+
self.Rmin = radius_min
95+
self.Rmax = radius_max
96+
self.gamma = gamma
97+
self.q_out = q_out
98+
self.discfraction = discfraction
99+
100+
self.a = self.thermalpower
101+
self.a2 = self.thermalpower/2
102+
self.g = densitypower
103+
self.g2 = 2 - densitypower
104+
self.k_out = ((1 + discfraction) / Rmax**3)**0.5
105+
self.sigma_out = self.g2 * discfraction / (
106+
2 * numpy.pi * Rmax**self.g * (Rmax**self.g2 - Rmin**self.g2))
107+
self.cs_out = self.q_out * numpy.pi * self.sigma_out / self.k_out
108+
109+
self.base_cylinder = uniform_unit_cylinder(targetN, base_grid)
110+
111+
def sigma(self, r):
112+
return self.sigma_out * (self.Rmax / r)**self.g
113+
114+
def csound(self, r):
115+
return self.cs_out * (self.Rmax / r)**self.a2
116+
117+
def cmass(self, r):
118+
return self.discfraction * (r**self.g2 - self.Rmin**self.g2) / (
119+
self.Rmax**self.g2 - self.Rmin**self.g2)
120+
121+
def mass_encl(self, r):
122+
return 1 + self.cmass(r)
123+
124+
def kappa(self, r):
125+
return (self.mass_encl(r) / r**3)**0.5
126+
127+
def toomreQ(self, r):
128+
return self.csound(r) * self.kappa(r) / numpy.pi / self.sigma(r)
129+
130+
def getradius(self, f):
131+
return (
132+
(self.Rmax**self.g2 - self.Rmin**self.g2) * f + self.Rmin**self.g2)**(
133+
1./self.g2)
134+
135+
def zscale(self, r):
136+
return self.csound(r) / self.kappa(r)
137+
138+
def u(self, r):
139+
if self.gamma == 1.:
128140
return self.csound(r)**2
129141
else:
130-
return self.csound(r)**2/(self.gamma-1)
142+
return self.csound(r)**2 / (self.gamma - 1)
143+
144+
def vcirc(self, r):
145+
return (self.mass_encl(r) / r)**0.5
131146

132-
def vcirc(self,r):
133-
return (self.mass_encl(r)/r)**0.5
134-
135147
def new_model(self):
136-
x,y,z=self.base_cylinder.make_xyz()
137-
self.actualN=len(x)
138-
f=x**2+y**2
139-
r=f**0.5
140-
rtarget=self.getradius(f)
141-
142-
mass=self.discfraction*numpy.ones_like(x)/self.actualN
143-
internal_energy=self.u(rtarget)
144-
zscale=self.zscale(rtarget)
145-
r=r.clip(1.e-8,2.)
146-
x=x/r
147-
y=y/r
148-
149-
vx=-y*self.vcirc(rtarget)
150-
vy=x*self.vcirc(rtarget)
151-
vz=numpy.zeros_like(x)
152-
153-
x=rtarget*x
154-
y=rtarget*y
155-
z=approximate_inverse_error_function(z)*zscale*2.**0.5
148+
x, y, z = self.base_cylinder.make_xyz()
149+
self.actualN = len(x)
150+
f = x**2 + y**2
151+
r = f**0.5
152+
rtarget = self.getradius(f)
153+
154+
mass = self.discfraction * numpy.ones_like(x) / self.actualN
155+
internal_energy = self.u(rtarget)
156+
zscale = self.zscale(rtarget)
157+
r = r.clip(1.e-8, 2.)
158+
x= x / r
159+
y= y / r
160+
161+
vx = -y * self.vcirc(rtarget)
162+
vy = x * self.vcirc(rtarget)
163+
vz = numpy.zeros_like(x)
164+
165+
x = rtarget * x
166+
y = rtarget * y
167+
z = approximate_inverse_error_function(z) * zscale * 2.**0.5
156168

157169
return (mass, x, y, z, vx, vy, vz, internal_energy)
158-
170+
159171
@property
160172
def result(self):
161-
masses, x,y,z, vx,vy,vz, internal_energies = self.new_model()
173+
masses, x, y, z, vx, vy, vz, internal_energies = self.new_model()
162174
result = Particles(self.actualN)
163175
result.mass = nbody_system.mass.new_quantity(masses)
164176
result.x = nbody_system.length.new_quantity(x)
@@ -168,10 +180,10 @@ def result(self):
168180
result.vy = nbody_system.speed.new_quantity(vy)
169181
result.vz = nbody_system.speed.new_quantity(vz)
170182
result.u = nbody_system.specific_energy.new_quantity(internal_energies)
171-
183+
172184
if not self.convert_nbody is None:
173-
result = ParticlesWithUnitsConverted(result, self.convert_nbody.as_converter_from_si_to_generic())
185+
result = ParticlesWithUnitsConverted(
186+
result, self.convert_nbody.as_converter_from_si_to_generic())
174187
result = result.copy()
175-
176-
return result
177188

189+
return result

0 commit comments

Comments
 (0)