diff --git a/src/gray_core.f90 b/src/gray_core.f90 index ff3cba1..3ec99b0 100644 --- a/src/gray_core.f90 +++ b/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 - ! Compute the approximate position of the absorption peak + if (params%raytracing%realtime) then + ! Check whether we are past the absorption peak block - use utils, only : vmax, parabola_vertex - real(wp_) :: power, peak(2) - integer :: index + logical :: past_peak + integer :: tail + ! 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 - call vmax(dpds, i, power, index) + if (past_peak .or. pow <= 0.1_wp_*p0ray(1)) then + ! 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 - peak = parabola_vertex(psjki(1, index - 1:index + 1), & - dpds(index - 1:index + 1)) + ! Find maximum in dP/ds + call vmax(dpds, i, power, index) + + ! 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