calc_mass_potential.py: update method gen_radius()

1 files changed, 23 insertions, 3 deletions

diff --git a/calc_mass_potential.py b/calc_mass_potential.py
index 22c3ccd..48a20d4 100755
--- a/calc_mass_potential.py
+++ b/calc_mass_potential.py
@@ -4,6 +4,10 @@
 # Created: 2016-06-24
 # Updated: 2016-06-24
 #
+# Change logs:
+# 2016-06-24:
+#   * Update method 'gen_radius()'
+#
 
 """
 Calculate the (gas and gravitational) mass profile and gravitational
@@ -207,15 +211,30 @@ class DensityProfile:
             bounds_error=False, fill_value=self.t_value[-1],
             assume_sorted=True)
 
-    def gen_radius(self, num=100):
+    def gen_radius(self, num=1000):
         """
         Generate radial points for following mass and potential calculation.
 
         The generated radial points are logarithmic-evenly distributed.
+
+        NOTE:
+        The radii are first generated to determine the inner-most bin width,
+        which is used to further divide the original inner-most bin (i.e.,
+        radius 0 - r_out[0]), and then the other radii are generated with
+        the constraint of given total number of points.
         """
         rout = self.r + self.r_err
-        rout_new = np.logspace(np.log10(rout[0]), np.log10(rout[-1]),
+        # first pass to determine the inner-most bin width
+        rout_tmp = np.logspace(np.log10(rout[0]), np.log10(rout[-1]),
                                num=num, base=10.0)
+        bw = rout_tmp[1] - rout_tmp[0]
+        # linear-evenly divide the first original bin (0 - rout[0])
+        nbin = int(np.ceil(rout[0] / bw))
+        rout_new1 = np.linspace(0.0, rout[0], num=nbin, endpoint=False)[1:]
+        # second pass to generate the other radii
+        rout_new2 = np.logspace(np.log10(rout[0]), np.log10(rout[-1]),
+                                num=(num-nbin+1), base=10.0)
+        rout_new = np.concatenate([rout_new1, rout_new2])
         rin_new = np.concatenate([[0.0], rout_new[:-1]])
         self.radius = (rout_new + rin_new) / 2.0
         self.radius_err = (rout_new - rin_new) / 2.0
@@ -232,7 +251,8 @@ class DensityProfile:
         #
         m_gas = np.zeros(self.radius.shape)
         if verbose:
-            print("Calculating the gas mass profile (#%d) ... " % len(m_gas))
+            print("Calculating the gas mass profile (#%d): ..." % len(m_gas),
+                  end="", flush=True)
         for i, r in enumerate(self.radius):
             if verbose and (i+1) % 10 == 0:
                 print("%d..." % (i+1), end="", flush=True)