Skip to content

Commit ea09412

Browse files
committed
RprofsAveraged: support variable dt
1 parent 83c9f7d commit ea09412

File tree

1 file changed

+20
-12
lines changed

1 file changed

+20
-12
lines changed

src/stagpy/stagyydata.py

Lines changed: 20 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -234,6 +234,8 @@ class RprofsAveraged(step.Rprofs):
234234
235235
The [`StepsView.rprofs_averaged`][stagpy.stagyydata.StepsView.rprofs_averaged]
236236
attribute is an instance of this class.
237+
238+
The current implementation does not take into account time-changing geometry.
237239
"""
238240

239241
steps: StepsView
@@ -246,21 +248,27 @@ def _steps_with_rprofs(self) -> StepsView:
246248
def _cached_data(self) -> dict[str, dt.Rprof]:
247249
return {}
248250

251+
@cached_property
252+
def _times(self) -> NDArray[np.float64]:
253+
return np.fromiter((s.time for s in self._steps_with_rprofs), dtype=float)
254+
255+
@cached_property
256+
def _dtimes(self) -> NDArray[np.float64]:
257+
midpoints = (self._times[:-1] + self._times[1:]) / 2
258+
return np.diff(midpoints, prepend=self._times[0], append=self._times[-1])
259+
249260
def __getitem__(self, name: str) -> dt.Rprof:
250-
# the averaging method has two shortcomings:
251-
# - does not take into account time changing geometry;
252-
# - does not take into account time changing timestep.
253261
if name in self._cached_data:
254262
return self._cached_data[name]
255-
steps_iter = iter(self._steps_with_rprofs)
256-
rpf = next(steps_iter).rprofs[name]
257-
rprof = np.copy(rpf.values)
258-
nprofs = 1
259-
for step_ in steps_iter:
260-
nprofs += 1
261-
rprof += step_.rprofs[name].values
262-
rprof /= nprofs
263-
self._cached_data[name] = dt.Rprof(rprof, rpf.rad, rpf.meta)
263+
integral_prof = sum(
264+
dtime * s.rprofs[name].values
265+
for dtime, s in zip(self._dtimes, self._steps_with_rprofs, strict=True)
266+
)
267+
self._cached_data[name] = dt.Rprof(
268+
values=integral_prof / (self._times[-1] - self._times[0]),
269+
rad=self._first_rprofs[name].rad,
270+
meta=self._first_rprofs[name].meta,
271+
)
264272
return self._cached_data[name]
265273

266274
@property

0 commit comments

Comments
 (0)