From da0342963dcdc4b7bc2d586894771a74a2a7c191 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Tue, 16 Apr 2024 19:25:10 +0200 Subject: [PATCH 1/2] Integral model: separate calculation at origin for spectral density --- src/gstools/covmodel/models.py | 19 ++++++++++++++----- 1 file changed, 14 insertions(+), 5 deletions(-) diff --git a/src/gstools/covmodel/models.py b/src/gstools/covmodel/models.py index ab94e279..7811085f 100644 --- a/src/gstools/covmodel/models.py +++ b/src/gstools/covmodel/models.py @@ -465,16 +465,25 @@ def spectral_density(self, k): # noqa: D102 if self.nu > 50.0: return ( (0.5 * self.len_rescaled / np.sqrt(np.pi)) ** self.dim + * (self.nu / (self.nu + self.dim)) * np.exp(-x) - * self.nu - / (self.nu + self.dim) * (1.0 + 2 * x / (self.nu + self.dim + 2)) ) - return ( + # separate calculation at origin + res = np.empty_like(k) + x_gz = np.logical_not(np.isclose(x, 0)) + x0 = x[x_gz] + k0 = k[x_gz] + # limit at k=0 + res[np.logical_not(x_gz)] = ( + 0.5 * self.len_rescaled / np.sqrt(np.pi) + ) ** self.dim * (self.nu / (self.nu + self.dim)) + res[x_gz] = ( self.nu - / (x ** (self.nu * 0.5) * 2 * (k * np.sqrt(np.pi)) ** self.dim) - * inc_gamma_low((self.nu + self.dim) / 2.0, x) + / (x0 ** (self.nu * 0.5) * 2 * (k0 * np.sqrt(np.pi)) ** self.dim) + * inc_gamma_low((self.nu + self.dim) / 2.0, x0) ) + return res def calc_integral_scale(self): # noqa: D102 return ( From cf2ae456688cfcf8783e5b2abf1ada44f962de85 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Wed, 17 Apr 2024 10:55:32 +0200 Subject: [PATCH 2/2] Integral model: make calculations more readable --- src/gstools/covmodel/models.py | 29 ++++++++++------------------- 1 file changed, 10 insertions(+), 19 deletions(-) diff --git a/src/gstools/covmodel/models.py b/src/gstools/covmodel/models.py index 7811085f..b1a9d68e 100644 --- a/src/gstools/covmodel/models.py +++ b/src/gstools/covmodel/models.py @@ -460,29 +460,20 @@ def cor(self, h): def spectral_density(self, k): # noqa: D102 k = np.asarray(k, dtype=np.double) - x = (k * self.len_rescaled / 2.0) ** 2 + fac = (0.5 * self.len_rescaled / np.sqrt(np.pi)) ** self.dim + lim = fac * self.nu / (self.nu + self.dim) # for nu > 50 we just use an approximation of the gaussian model if self.nu > 50.0: - return ( - (0.5 * self.len_rescaled / np.sqrt(np.pi)) ** self.dim - * (self.nu / (self.nu + self.dim)) - * np.exp(-x) - * (1.0 + 2 * x / (self.nu + self.dim + 2)) - ) + x = (k * self.len_rescaled / 2) ** 2 + return lim * np.exp(-x) * (1 + 2 * x / (self.nu + self.dim + 2)) # separate calculation at origin + s = (self.nu + self.dim) / 2 res = np.empty_like(k) - x_gz = np.logical_not(np.isclose(x, 0)) - x0 = x[x_gz] - k0 = k[x_gz] - # limit at k=0 - res[np.logical_not(x_gz)] = ( - 0.5 * self.len_rescaled / np.sqrt(np.pi) - ) ** self.dim * (self.nu / (self.nu + self.dim)) - res[x_gz] = ( - self.nu - / (x0 ** (self.nu * 0.5) * 2 * (k0 * np.sqrt(np.pi)) ** self.dim) - * inc_gamma_low((self.nu + self.dim) / 2.0, x0) - ) + k_gz = np.logical_not(np.isclose(k, 0)) + x = (k[k_gz] * self.len_rescaled / 2) ** 2 + # limit at k=0 (inc_gamma_low(s, x) / x**s -> 1/s for x -> 0) + res[np.logical_not(k_gz)] = lim + res[k_gz] = 0.5 * self.nu * fac / x**s * inc_gamma_low(s, x) return res def calc_integral_scale(self): # noqa: D102