From e6b727785fab339103a103194317ffdc6ba37684 Mon Sep 17 00:00:00 2001 From: Aaron LI Date: Tue, 25 Jul 2017 20:38:32 +0800 Subject: clusters/solver.py: Use multiple data points to fit the power-law Signed-off-by: Aaron LI --- fg21sim/extragalactic/clusters/solver.py | 36 +++++++++++++++++--------------- 1 file changed, 19 insertions(+), 17 deletions(-) (limited to 'fg21sim/extragalactic/clusters') diff --git a/fg21sim/extragalactic/clusters/solver.py b/fg21sim/extragalactic/clusters/solver.py index d7660b0..67bda58 100644 --- a/fg21sim/extragalactic/clusters/solver.py +++ b/fg21sim/extragalactic/clusters/solver.py @@ -322,30 +322,32 @@ class FokkerPlanckSolver: def fix_boundary(self, uc): """ - Truncate the lower end (i.e., near the lower boundary) of the - distribution/spectrum and then extrapolate as a power law, in order - to avoid the unphysical pile-up of electrons at the lower regime. - - TODO: - Fit a power law to the same number (``buffer_np``) of data points, - then extrapolate it to fix the lower buffer region. + Due to the no-flux boundary condition adopted, particles may + unphysically pile up near the lower boundary. Therefore, a + buffer region spanning ``self.buffer_np`` cells is chosen, within + which the densities are replaced by extrapolating from the upper + density distribution as a power law, and the power-law index + is determined by fitting to the data points of ``self.buffer_np`` + cells on the upper side of the buffer region. References: Ref.[2],Sec.(3.3) """ if self.buffer_np is None: return uc + if (uc <= 0.0).sum() > 0: + logger.warning("solved density has negative values!") + return uc - uc = np.asarray(uc) x = self.x - # Calculate the power-law index - xa = x[self.buffer_np] - xb = x[self.buffer_np+1] - ya = uc[self.buffer_np] - yb = uc[self.buffer_np+1] - if ya > 0 and yb > 0: - # Truncate and extrapolate as a power law - s = np.log(yb/ya) / np.log(xb/xa) - uc[:self.buffer_np] = ya * (x[:self.buffer_np] / xa) ** s + # Calculate the power-law index from ``self.buffer_np`` data + # points just right of the buffer region. + xp = x[self.buffer_np:(self.buffer_np*2)] + yp = uc[self.buffer_np:(self.buffer_np*2)] + # Power-law fit + pfit = np.polyfit(np.log10(xp), np.log10(yp), deg=1) + xbuf = x[:self.buffer_np] + ybuf = 10 ** np.polyval(pfit, np.log10(xbuf)) + uc[:self.buffer_np] = ybuf return uc def time_step(self): -- cgit v1.2.2