src/gray_core.f90: refine realtime termination check

In some very rare cases the residual power check can stop the
integration too early.

For example, this extremely unlikely chain of events was observed in a
few cases while performing a parameter sweep (~3000 simulations):

  1. n⊥ converges on wrong branch ⇒ α jumps ⇒ dP/ds jumps (common)
  2. dP/ds jumps close to the peak (unlikely)
  3. the peak of dP/ds is close to 0.6P₀ (unlikely)
  4. dP/ds jumps to dP/ds_max + 10⁻¹² ⇒ wrong peak found (very unlikely)
  5. peak is at exactly the current point
     ⇒ three-point interpolation fails (unlikely)

src/gray_core.f90

@@ -529,24 +529,36 @@ contains
             end do
           end if
 
-          if (params%raytracing%realtime .and. pow <= 0.4_wp_*p0ray(1)) then
+          if (params%raytracing%realtime) then
-            ! Compute the approximate position of the absorption peak
+            ! Check whether we are past the absorption peak
             block
-              use utils, only : vmax, parabola_vertex
+              logical :: past_peak
-              real(wp_) :: power, peak(2)
+              integer :: tail
-              integer   :: index
+              ! We assume so if dP/ds has been decreasing in the last 8 steps
+              tail = max(2, i-8)
+              past_peak = i > 8 .and. all(dpds(tail:i) - dpds(tail-1:i-1) < 0)
 
-              ! Find maximum in dP/ds
+              if (past_peak .or. pow <= 0.1_wp_*p0ray(1)) then
-              call vmax(dpds, i, power, index)
+                ! Compute the approximate position of the absorption peak
+                block
+                  use utils, only : vmax, parabola_vertex
+                  real(wp_) :: power, peak(2)
+                  integer   :: index
 
-              ! Interpolate the peak with a parabola
+                  ! Find maximum in dP/ds
-              peak = parabola_vertex(psjki(1, index - 1:index + 1), &
+                  call vmax(dpds, i, power, index)
-                                     dpds(index - 1:index + 1))
+
+                  ! Interpolate the peak with a parabola
+                  peak = parabola_vertex(psjki(1, index - 1:index + 1), &
+                                         dpds(index - 1:index + 1))
+
+                  results%rho_peak = sqrt(peak(1))
+                end block
+                ! Stop propagation loop
+                exit
+              end if
 
-              results%rho_peak = sqrt(peak(1))
             end block
-            ! Stop propagation loop
-            exit
           end if
 
           ! print ray positions for j=nrayr in local reference system