diff --git a/tardis/montecarlo/src/cmontecarlo.c b/tardis/montecarlo/src/cmontecarlo.c index 357a4c32506..714cb284aa2 100644 --- a/tardis/montecarlo/src/cmontecarlo.c +++ b/tardis/montecarlo/src/cmontecarlo.c @@ -137,67 +137,62 @@ void calculate_chi_bf(rpacket_t * packet, storage_model_t * storage) rpacket_set_chi_boundfree(packet, bf_helper * doppler_factor); } -double +void compute_distance2boundary (rpacket_t * packet, const storage_model_t * storage) { double r = rpacket_get_r (packet); double mu = rpacket_get_mu (packet); double r_outer = storage->r_outer[rpacket_get_current_shell_id (packet)]; double r_inner = storage->r_inner[rpacket_get_current_shell_id (packet)]; - double d_outer = - sqrt (r_outer * r_outer + ((mu * mu - 1.0) * r * r)) - (r * mu); - double d_inner; + double check, distance; if (rpacket_get_recently_crossed_boundary (packet) == 1) { rpacket_set_next_shell_id (packet, 1); - return d_outer; + distance = sqrt (r_outer * r_outer + ((mu * mu - 1.0) * r * r)) - (r * mu); + } + else if (mu > 0.0) + { // direction outward + rpacket_set_next_shell_id (packet, 1); + distance = sqrt (r_outer * r_outer + ((mu * mu - 1.0) * r * r)) - (r * mu); } else - { - double check = r_inner * r_inner + (r * r * (mu * mu - 1.0)); - if (check < 0.0) - { - rpacket_set_next_shell_id (packet, 1); - return d_outer; + { // going inward + if ( (check = r_inner * r_inner + (r * r * (mu * mu - 1.0)) )>= 0.0) + { // hit inner boundary + rpacket_set_next_shell_id (packet, -1); + distance = - r * mu - sqrt (check); } else - { - d_inner = mu < 0.0 ? -r * mu - sqrt (check) : MISS_DISTANCE; + { // miss inner boundary + rpacket_set_next_shell_id (packet, 1); + distance = sqrt (r_outer * r_outer + ((mu * mu - 1.0) * r * r)) - (r * mu); } } - if (d_inner < d_outer) - { - rpacket_set_next_shell_id (packet, -1); - return d_inner; - } - else - { - rpacket_set_next_shell_id (packet, 1); - return d_outer; - } + rpacket_set_d_boundary (packet, distance); } tardis_error_t -compute_distance2line (const rpacket_t * packet, const storage_model_t * storage, - double *result) +compute_distance2line (rpacket_t * packet, const storage_model_t * storage) { - tardis_error_t ret_val = TARDIS_ERROR_OK; - if (rpacket_get_last_line (packet)) - { - *result = MISS_DISTANCE; - } - else + if (!rpacket_get_last_line (packet)) { double r = rpacket_get_r (packet); double mu = rpacket_get_mu (packet); double nu = rpacket_get_nu (packet); double nu_line = rpacket_get_nu_line (packet); + double distance, nu_diff; double t_exp = storage->time_explosion; double inverse_t_exp = storage->inverse_time_explosion; int64_t cur_zone_id = rpacket_get_current_shell_id (packet); double doppler_factor = 1.0 - mu * r * inverse_t_exp * INVERSE_C; double comov_nu = nu * doppler_factor; - if (comov_nu < nu_line) + if ( (nu_diff = comov_nu - nu_line) >= 0) + { + distance = (nu_diff / nu) * C * t_exp; + rpacket_set_d_line (packet, distance); + return TARDIS_ERROR_OK; + } + else { if (rpacket_get_next_line_id (packet) == storage->no_of_lines - 1) { @@ -232,14 +227,14 @@ compute_distance2line (const rpacket_t * packet, const storage_model_t * storage fprintf (stderr, "nu = %f\n", nu); fprintf (stderr, "doppler_factor = %f\n", doppler_factor); fprintf (stderr, "cur_zone_id = %" PRIi64 "\n", cur_zone_id); - ret_val = TARDIS_ERROR_COMOV_NU_LESS_THAN_NU_LINE; - } - else - { - *result = ((comov_nu - nu_line) / nu) * C * t_exp; + return TARDIS_ERROR_COMOV_NU_LESS_THAN_NU_LINE; } } - return ret_val; + else + { + rpacket_set_d_line (packet, MISS_DISTANCE); + return TARDIS_ERROR_OK; + } } void @@ -712,12 +707,9 @@ montecarlo_compute_distances (rpacket_t * packet, storage_model_t * storage) } else { - rpacket_set_d_boundary (packet, - compute_distance2boundary (packet, storage)); - double d_line; - compute_distance2line (packet, storage, &d_line); + compute_distance2boundary(packet, storage); + compute_distance2line (packet, storage); // FIXME MR: return status of compute_distance2line() is ignored - rpacket_set_d_line (packet, d_line); compute_distance2continuum (packet, storage); } } diff --git a/tardis/montecarlo/src/cmontecarlo.h b/tardis/montecarlo/src/cmontecarlo.h index b22b658e813..8f2b5676671 100644 --- a/tardis/montecarlo/src/cmontecarlo.h +++ b/tardis/montecarlo/src/cmontecarlo.h @@ -33,7 +33,7 @@ double rpacket_doppler_factor(const rpacket_t *packet, const storage_model_t *st * * @return distance to shell boundary */ -double compute_distance2boundary (rpacket_t * packet, +void compute_distance2boundary (rpacket_t * packet, const storage_model_t * storage); /** Calculate the distance the packet has to travel until it redshifts to the first spectral line. @@ -43,9 +43,8 @@ double compute_distance2boundary (rpacket_t * packet, * * @return distance to the next spectral line */ -tardis_error_t compute_distance2line (const rpacket_t * packet, - const storage_model_t * storage, - double *result); +tardis_error_t compute_distance2line (rpacket_t * packet, + const storage_model_t * storage); /** Calculate the distance to the next continuum event, which can be a Thomson scattering, bound-free absorption or free-free transition. diff --git a/tardis/montecarlo/src/test_cmontecarlo.c b/tardis/montecarlo/src/test_cmontecarlo.c index 80fa06fb190..6657c14392c 100644 --- a/tardis/montecarlo/src/test_cmontecarlo.c +++ b/tardis/montecarlo/src/test_cmontecarlo.c @@ -220,10 +220,10 @@ test_compute_distance2boundary(void){ storage_model_t sm; init_rpacket(&rp); init_storage_model(&sm); - double D_BOUNDARY = compute_distance2boundary(&rp, &sm); - rpacket_set_d_boundary(&rp, D_BOUNDARY); + compute_distance2boundary(&rp, &sm); + double D_BOUNDARY = rpacket_get_d_boundary(&rp); dealloc_storage_model(&sm); - return D_BOUNDARY; + return D_BOUNDARY; } double test_compute_distance2line(void){ @@ -231,12 +231,11 @@ test_compute_distance2line(void){ storage_model_t sm; init_rpacket(&rp); init_storage_model(&sm); - double D_LINE; // FIXME MR: return status of compute_distance2line() is ignored - compute_distance2line(&rp, &sm, &D_LINE); - rpacket_set_d_line(&rp, D_LINE); + compute_distance2line(&rp, &sm); + double D_LINE = rpacket_get_d_line(&rp); dealloc_storage_model(&sm); - return D_LINE; + return D_LINE; } double @@ -279,17 +278,12 @@ test_increment_j_blue_estimator(void){ storage_model_t sm; init_rpacket(&rp); init_storage_model(&sm); - int64_t j_blue_idx = 0; - double D_BOUNDARY = compute_distance2boundary(&rp, &sm); - rpacket_set_d_boundary(&rp, D_BOUNDARY); - double D_LINE; - // FIXME MR: return status of compute_distance2line() is ignored - compute_distance2line(&rp, &sm, &D_LINE); - rpacket_set_d_line(&rp, D_LINE); + int64_t j_blue_idx = 0; + compute_distance2line(&rp, &sm); move_packet(&rp, &sm, 1e13); - double d_line = rpacket_get_d_line(&rp); - increment_j_blue_estimator(&rp, &sm, d_line, j_blue_idx); - double res = sm.line_lists_j_blues[j_blue_idx]; + double d_line = rpacket_get_d_line(&rp); + increment_j_blue_estimator(&rp, &sm, d_line, j_blue_idx); + double res = sm.line_lists_j_blues[j_blue_idx]; dealloc_storage_model(&sm); return res; } @@ -379,12 +373,7 @@ test_calculate_chi_bf(void){ rk_state mt_state; irandom(&mt_state); int64_t j_blue_idx = 0; - double D_BOUNDARY = compute_distance2boundary(&rp, &sm); - rpacket_set_d_boundary(&rp, D_BOUNDARY); - double D_LINE; - // FIXME MR: return status of compute_distance2line() is ignored - compute_distance2line(&rp, &sm, &D_LINE); - rpacket_set_d_line(&rp, D_LINE); + compute_distance2line(&rp, &sm); move_packet(&rp, &sm, 1e13); double d_line = rpacket_get_d_line(&rp); increment_j_blue_estimator(&rp, &sm, d_line, j_blue_idx);