diff --git a/gammapkt.cc b/gammapkt.cc index 3f49c1726..d66b989a5 100644 --- a/gammapkt.cc +++ b/gammapkt.cc @@ -1000,15 +1000,15 @@ void wollaeger_thermalisation(Packet &pkt) { bool end_packet = false; while (!end_packet) { // distance to the next cell - const auto [sdist, snext] = - grid::boundary_distance(pkt_copy.dir, pkt_copy.pos, pkt_copy.prop_time, pkt_copy.where, &pkt_copy.last_cross); + const auto [sdist, snext] = grid::boundary_distance(pkt_copy.dir, pkt_copy.pos, pkt_copy.prop_time, pkt_copy.where, + &pkt_copy.last_cross, true); const double s_cont = sdist * t_current * t_current * t_current / std::pow(pkt_copy.prop_time, 3); const int mgi = grid::get_cell_modelgridindex(pkt_copy.where); if (mgi != grid::get_npts_model()) { tau += grid::get_rho(mgi) * s_cont * mean_gamma_opac; // contribution to the integral } // move packet copy now - move_pkt_withtime(pkt_copy, sdist); + move_pkt_withtime<100 * CLIGHT_PROP>(pkt_copy, sdist); grid::change_cell(pkt_copy, snext); end_packet = (pkt_copy.type == TYPE_ESCAPE); @@ -1061,15 +1061,16 @@ void guttman_thermalisation(Packet &pkt) { bool end_packet = false; while (!end_packet) { // distance to the next cell - const auto [sdist, snext] = - grid::boundary_distance(pkt_copy.dir, pkt_copy.pos, pkt_copy.prop_time, pkt_copy.where, &pkt_copy.last_cross); + const auto [sdist, snext] = grid::boundary_distance(pkt_copy.dir, pkt_copy.pos, pkt_copy.prop_time, + pkt_copy.where, &pkt_copy.last_cross, true); const double s_cont = sdist * std::pow(t, 3.) / std::pow(pkt_copy.prop_time, 3.); + std::cout << "t: " << t << " s_cont: " << s_cont << "\n"; const int mgi = grid::get_cell_modelgridindex(pkt_copy.where); if (mgi != grid::get_npts_model()) { column_densities[i] += grid::get_rho_tmin(mgi) * s_cont; // contribution to the integral } // move packet copy now - move_pkt_withtime(pkt_copy, sdist); + move_pkt_withtime<100 * CLIGHT_PROP>(pkt_copy, sdist); grid::change_cell(pkt_copy, snext); end_packet = (pkt_copy.type == TYPE_ESCAPE); diff --git a/grid.cc b/grid.cc index 3ddc27b4c..132be6ed6 100644 --- a/grid.cc +++ b/grid.cc @@ -2465,10 +2465,11 @@ template static auto get_coordboundary_distances_cylindrical2d( const std::array pkt_pos, const std::array pkt_dir, const std::array pktposgridcoord, const std::array pktvelgridcoord, int cellindex, - const double tstart, - const std::array cellcoordmax) -> std::tuple, std::array> { + const double tstart, const std::array cellcoordmax, + bool neglect_exp) -> std::tuple, std::array> { // to get the cylindrical intersection, get the spherical intersection with Z components set to zero, and the // propagation speed set to the xy component of the 3-velocity + double CLIGHT_PROP_DIST_CALC = (neglect_exp) ? 100 * CLIGHT_PROP : CLIGHT_PROP; std::array d_coordminboundary{}; std::array d_coordmaxboundary{}; @@ -2476,7 +2477,7 @@ static auto get_coordboundary_distances_cylindrical2d( const std::array posnoz = {pkt_pos[0], pkt_pos[1]}; const double dirxylen = std::sqrt((pkt_dir[0] * pkt_dir[0]) + (pkt_dir[1] * pkt_dir[1])); - const double xyspeed = dirxylen * CLIGHT_PROP; // r_cyl component of velocity + const double xyspeed = dirxylen * CLIGHT_PROP_DIST_CALC; // r_cyl component of velocity // make a normalised direction vector in the xy plane const std::array dirnoz = {pkt_dir[0] / dirxylen, pkt_dir[1] / dirxylen}; @@ -2487,7 +2488,7 @@ static auto get_coordboundary_distances_cylindrical2d( if (r_inner > 0) { const double d_rcyl_coordminboundary = expanding_shell_intersection(posnoz, dirnoz, xyspeed, r_inner, true, tstart); if (d_rcyl_coordminboundary >= 0) { - const double d_z_coordminboundary = d_rcyl_coordminboundary / xyspeed * pkt_dir[2] * CLIGHT_PROP; + const double d_z_coordminboundary = d_rcyl_coordminboundary / xyspeed * pkt_dir[2] * CLIGHT_PROP_DIST_CALC; d_coordminboundary[0] = std::sqrt(d_rcyl_coordminboundary * d_rcyl_coordminboundary + d_z_coordminboundary * d_z_coordminboundary); } @@ -2497,7 +2498,7 @@ static auto get_coordboundary_distances_cylindrical2d( const double d_rcyl_coordmaxboundary = expanding_shell_intersection(posnoz, dirnoz, xyspeed, r_outer, false, tstart); d_coordmaxboundary[0] = -1; if (d_rcyl_coordmaxboundary >= 0) { - const double d_z_coordmaxboundary = d_rcyl_coordmaxboundary / xyspeed * pkt_dir[2] * CLIGHT_PROP; + const double d_z_coordmaxboundary = d_rcyl_coordmaxboundary / xyspeed * pkt_dir[2] * CLIGHT_PROP_DIST_CALC; d_coordmaxboundary[0] = std::sqrt(d_rcyl_coordmaxboundary * d_rcyl_coordmaxboundary + d_z_coordmaxboundary * d_z_coordmaxboundary); } @@ -2507,21 +2508,23 @@ static auto get_coordboundary_distances_cylindrical2d( ((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_coordminboundary[1] = CLIGHT_PROP_DIST_CALC * 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; + d_coordmaxboundary[1] = CLIGHT_PROP_DIST_CALC * t_zcoordmaxboundary; return {d_coordminboundary, d_coordmaxboundary}; } [[nodiscard]] auto boundary_distance(const std::array dir, const std::array pos, - const double tstart, const int cellindex, - enum cell_boundary *pkt_last_cross) -> std::tuple + const double tstart, const int cellindex, enum cell_boundary *pkt_last_cross, + bool neglect_exp) -> std::tuple /// Basic routine to compute distance to a cell boundary. { + double CLIGHT_PROP_DIST_CALC = (neglect_exp) ? 100 * CLIGHT_PROP : CLIGHT_PROP; + if constexpr (FORCE_SPHERICAL_ESCAPE_SURFACE) { if (get_cell_r_inner(cellindex) > globals::vmax * globals::tmin) { return {0., -99}; @@ -2545,18 +2548,18 @@ static auto get_coordboundary_distances_cylindrical2d( if constexpr (GRID_TYPE == GRID_CARTESIAN3D) { // keep xyz Cartesian coordinates for (int d = 0; d < ndim; d++) { - pktvelgridcoord[d] = dir[d] * CLIGHT_PROP; + pktvelgridcoord[d] = dir[d] * CLIGHT_PROP_DIST_CALC; } } else if constexpr (GRID_TYPE == GRID_CYLINDRICAL2D) { // xy plane radial velocity - pktvelgridcoord[0] = (pos[0] * dir[0] + pos[1] * dir[1]) / pktposgridcoord[0] * CLIGHT_PROP; + pktvelgridcoord[0] = (pos[0] * dir[0] + pos[1] * dir[1]) / pktposgridcoord[0] * CLIGHT_PROP_DIST_CALC; // second cylindrical coordinate is z - pktvelgridcoord[1] = dir[2] * CLIGHT_PROP; + pktvelgridcoord[1] = dir[2] * CLIGHT_PROP_DIST_CALC; } else if constexpr (GRID_TYPE == GRID_SPHERICAL1D) { // the only coordinate is radius from the origin - pktvelgridcoord[0] = dot(pos, dir) / pktposgridcoord[0] * CLIGHT_PROP; // radial velocity + pktvelgridcoord[0] = dot(pos, dir) / pktposgridcoord[0] * CLIGHT_PROP_DIST_CALC; // radial velocity } else { assert_always(false); } @@ -2634,7 +2637,7 @@ static auto get_coordboundary_distances_cylindrical2d( std::array{-1}; // distance to reach the cell's lower boundary on each coordinate if constexpr (GRID_TYPE == GRID_SPHERICAL1D) { last_cross = BOUNDARY_NONE; // handle this separately by setting d_inner and d_outer negative for invalid direction - const double speed = vec_len(dir) * CLIGHT_PROP; // just in case dir is not normalised + const double speed = vec_len(dir) * CLIGHT_PROP_DIST_CALC; // just in case dir is not normalised const double r_inner = grid::get_cellcoordmin(cellindex, 0) * tstart / globals::tmin; d_coordminboundary[0] = (r_inner > 0.) ? expanding_shell_intersection(pos, dir, speed, r_inner, true, tstart) : -1.; @@ -2650,7 +2653,7 @@ static auto get_coordboundary_distances_cylindrical2d( } std::tie(d_coordminboundary, d_coordmaxboundary) = get_coordboundary_distances_cylindrical2d( - pos, dir, pktposgridcoord, pktvelgridcoord, cellindex, tstart, cellcoordmax); + pos, dir, pktposgridcoord, pktvelgridcoord, cellindex, tstart, cellcoordmax, neglect_exp); } else if constexpr (GRID_TYPE == GRID_CARTESIAN3D) { // There are six possible boundary crossings. Each of the three @@ -2670,13 +2673,13 @@ static auto get_coordboundary_distances_cylindrical2d( ((pktposgridcoord[d] - (pktvelgridcoord[d] * tstart)) / (grid::get_cellcoordmin(cellindex, d) - (pktvelgridcoord[d] * globals::tmin)) * globals::tmin) - tstart; - d_coordminboundary[d] = CLIGHT_PROP * t_coordminboundary; + d_coordminboundary[d] = CLIGHT_PROP_DIST_CALC * t_coordminboundary; const double t_coordmaxboundary = ((pktposgridcoord[d] - (pktvelgridcoord[d] * tstart)) / (cellcoordmax[d] - (pktvelgridcoord[d] * globals::tmin)) * globals::tmin) - tstart; - d_coordmaxboundary[d] = CLIGHT_PROP * t_coordmaxboundary; + d_coordmaxboundary[d] = CLIGHT_PROP_DIST_CALC * t_coordmaxboundary; } } else { assert_always(false); diff --git a/grid.h b/grid.h index 72bba2793..f09a6f019 100644 --- a/grid.h +++ b/grid.h @@ -130,7 +130,8 @@ void write_grid_restart_data(int timestep); [[nodiscard]] auto get_ndo_nonempty(int rank) -> int; [[nodiscard]] auto get_totmassradionuclide(int z, int a) -> double; [[nodiscard]] auto boundary_distance(std::array dir, std::array pos, double tstart, int cellindex, - enum cell_boundary *pkt_last_cross) -> std::tuple; + enum cell_boundary *pkt_last_cross, + bool neglect_exp = false) -> std::tuple; [[nodiscard]] inline auto get_elem_abundance(int modelgridindex, int element) -> float // mass fraction of an element (all isotopes combined) diff --git a/vectors.h b/vectors.h index df854cd74..b21332c49 100644 --- a/vectors.h +++ b/vectors.h @@ -141,6 +141,7 @@ template return doppler_nucmf_on_nurf(dir_rf, get_velocity(pos_rf, prop_time)); } +template constexpr auto move_pkt_withtime(std::span pos_rf, const std::array dir_rf, double &prop_time, const double nu_rf, double &nu_cmf, const double e_rf, double &e_cmf, const double distance) -> double @@ -150,7 +151,7 @@ constexpr auto move_pkt_withtime(std::span pos_rf, const std::array= 0); const double nu_cmf_old = nu_cmf; - prop_time += distance / CLIGHT_PROP; + prop_time += distance / CLIGHT_ASSUMED; pos_rf[0] += (dir_rf[0] * distance); pos_rf[1] += (dir_rf[1] * distance); @@ -170,8 +171,10 @@ constexpr auto move_pkt_withtime(std::span pos_rf, const std::array constexpr auto move_pkt_withtime(Packet &pkt, const double distance) -> double { - return move_pkt_withtime(pkt.pos, pkt.dir, pkt.prop_time, pkt.nu_rf, pkt.nu_cmf, pkt.e_rf, pkt.e_cmf, distance); + return move_pkt_withtime(pkt.pos, pkt.dir, pkt.prop_time, pkt.nu_rf, pkt.nu_cmf, pkt.e_rf, pkt.e_cmf, + distance); } [[nodiscard]] [[gnu::const]] constexpr auto get_arrive_time(const Packet &pkt) -> double