diff --git a/grid.cc b/grid.cc index 94f962cc6..0cff60990 100644 --- a/grid.cc +++ b/grid.cc @@ -1335,16 +1335,23 @@ auto get_coordboundary_distances_cylindrical2d(const std::array &pkt_ } // z boundaries are the same as Cartesian - const double t_zcoordminboundary = - ((pktposgridcoord[1] - (pktvelgridcoord[1] * tstart)) / - ((grid::get_cellcoordmin(cellindex, 1)) - (pktvelgridcoord[1] * globals::tmin)) * globals::tmin) - - tstart; - d_coordminboundary[1] = CLIGHT_PROP * t_zcoordminboundary; - const double t_zcoordmaxboundary = ((pktposgridcoord[1] - (pktvelgridcoord[1] * tstart)) / - ((cellcoordmax[1]) - (pktvelgridcoord[1] * globals::tmin)) * globals::tmin) - - tstart; - d_coordmaxboundary[1] = CLIGHT_PROP * t_zcoordmaxboundary; + if ((pktvelgridcoord[1] * tstart) > pktposgridcoord[1]) { + d_coordminboundary[1] = -1; + + const double t_zcoordmaxboundary = ((pktposgridcoord[1] - (pktvelgridcoord[1] * tstart)) / + ((cellcoordmax[1]) - (pktvelgridcoord[1] * globals::tmin)) * globals::tmin) - + tstart; + d_coordmaxboundary[1] = CLIGHT_PROP * t_zcoordmaxboundary; + } else { + const double t_zcoordminboundary = + ((pktposgridcoord[1] - (pktvelgridcoord[1] * tstart)) / + ((grid::get_cellcoordmin(cellindex, 1)) - (pktvelgridcoord[1] * globals::tmin)) * globals::tmin) - + tstart; + d_coordminboundary[1] = CLIGHT_PROP * t_zcoordminboundary; + + d_coordmaxboundary[1] = -1; + } return {d_coordminboundary, d_coordmaxboundary}; } @@ -2486,12 +2493,6 @@ auto get_totmassradionuclide(const int z, const int a) -> double { std::tie(d_coordminboundary, d_coordmaxboundary) = get_coordboundary_distances_cylindrical2d( pos, dir, pktposgridcoord, pktvelgridcoord, cellindex, tstart, cellcoordmax); - if ((pktvelgridcoord[1] * tstart) > pktposgridcoord[1]) { - d_coordminboundary[1] = -1; - } else { - d_coordmaxboundary[1] = -1; - } - } else if constexpr (GRID_TYPE == GridType::CARTESIAN3D) { // There are six possible boundary crossings. Each of the three // cartesian coordinates may be taken in turn. For x, the packet