diff options
Diffstat (limited to 'fg21sim/extragalactic/clusters')
| -rw-r--r-- | fg21sim/extragalactic/clusters/halo.py | 33 | ||||
| -rw-r--r-- | fg21sim/extragalactic/clusters/solver.py | 15 |
2 files changed, 27 insertions, 21 deletions
diff --git a/fg21sim/extragalactic/clusters/halo.py b/fg21sim/extragalactic/clusters/halo.py index 75789b6..ce4bcaf 100644 --- a/fg21sim/extragalactic/clusters/halo.py +++ b/fg21sim/extragalactic/clusters/halo.py @@ -182,7 +182,7 @@ class RadioHalo: Returns ------- - electron_spec : `~numpy.ndarray` + electron_spec : float 1D `~numpy.ndarray` The solved electron spectrum at ``zend``. Unit: [cm^-3] """ @@ -196,8 +196,7 @@ class RadioHalo: tstop = COSMO.age(zend) if n0_e is None: # Accumulated constantly injected electrons until ``tstart``. - n_inj = np.array([self.fp_injection(gm) - for gm in self.gamma]) + n_inj = self.fp_injection(self.gamma) n0_e = n_inj * tstart self.electron_spec = self.fpsolver.solve(u0=n0_e, tstart=tstart, @@ -218,8 +217,8 @@ class RadioHalo: Parameters ---------- - gamma : float - Lorentz factor of electrons + gamma : float, or float 1D `~numpy.ndarray` + Lorentz factors of electrons t : None Currently a constant injection rate is assumed, therefore this parameter is not used. Keep it for the consistency @@ -227,8 +226,8 @@ class RadioHalo: Returns ------- - Qe : float - Current electron injection rate at specified energy (gamma). + Qe : float, or float 1D `~numpy.ndarray` + Current electron injection rate at specified energies (gamma). Unit: [cm^-3 Gyr^-1] References @@ -251,16 +250,16 @@ class RadioHalo: Parameters ---------- - gamma : float - The Lorentz factor of electrons + gamma : float, or float 1D `~numpy.ndarray` + The Lorentz factors of electrons t : float Current (cosmic) time when solving the equation Unit: [Gyr] Returns ------- - diffusion : float - Diffusion coefficient + diffusion : float, or float 1D `~numpy.ndarray` + Diffusion coefficients Unit: [Gyr^-1] References @@ -283,8 +282,8 @@ class RadioHalo: Returns ------- - advection : float - Advection coefficient, describing the energy loss/gain rate. + advection : float, or float 1D `~numpy.ndarray` + Advection coefficients, describing the energy loss/gain rates. Unit: [Gyr^-1] """ advection = (abs(self._loss_ion(gamma, t)) + @@ -394,16 +393,16 @@ class RadioHalo: Parameters ---------- - gamma : float - The Lorentz factor of electrons + gamma : float, or float 1D `~numpy.ndarray` + The Lorentz factors of electrons t : float The cosmic time/age Unit: [Gyr] Returns ------- - loss : float - The energy loss rate + loss : float, or float 1D `~numpy.ndarray` + The energy loss rates Unit: [Gyr^-1] References diff --git a/fg21sim/extragalactic/clusters/solver.py b/fg21sim/extragalactic/clusters/solver.py index dea7f3f..d7660b0 100644 --- a/fg21sim/extragalactic/clusters/solver.py +++ b/fg21sim/extragalactic/clusters/solver.py @@ -104,6 +104,13 @@ class FokkerPlanckSolver: NOTE ---- + All above functions should accept two parameters: ``(x, t)``, + where ``x`` is an 1D float `~numpy.ndarray` representing the adopted + logarithmic grid points along the spatial/energy axis, ``t`` is the + time of each solving step. + + NOTE + ---- The diffusion coefficients (i.e., calculated by ``f_diffusion()``) should be *positive* (i.e., C(x) > 0), otherwise unstable or wrong results may occur, due to the current numerical scheme/algorithm @@ -281,9 +288,9 @@ class FokkerPlanckSolver: dx_phalf = self.dx_phalf dx_mhalf = self.dx_mhalf dt = self.tstep - B = np.array([self.f_advection(x_, tc) for x_ in x]) - C = np.array([self.f_diffusion(x_, tc) for x_ in x]) - Q = np.array([self.f_injection(x_, tc) for x_ in x]) + B = self.f_advection(x, tc) + C = self.f_diffusion(x, tc) + Q = self.f_injection(x, tc) # B_phalf = self.X_phalf(B) B_mhalf = self.X_mhalf(B) @@ -307,7 +314,7 @@ class FokkerPlanckSolver: b[-1] = 1 + (dt/dx[-1]) * (C_mhalf[-1]/dx_mhalf[-1])*Wplus_mhalf[-1] # Escape from the system if self.f_escape is not None: - T = np.array([self.f_escape(x_, tc) for x_ in x]) + T = self.f_escape(x, tc) b += dt / T # Right-hand side r = dt * Q + uc |