@@ -258,6 +258,7 @@ void ICPhoton<Electrons, Photons>::generate_grid() {
258258 logspace_boundary_center (std::log2 (gamma_min), std::log2 (gamma_max), gamma_size, gamma, dgamma);
259259
260260 IC_tab = MeshGrid ({gamma_size, nu_size}, -1 );
261+
261262 generated = true ;
262263
263264 I_nu_dnu = Array::from_shape ({nu_size});
@@ -287,36 +288,36 @@ Real ICPhoton<Electrons, Photons>::compute_I_nu(Real nu) {
287288
288289 for (size_t i = gamma_size; i-- > 0 ;) {
289290 const Real gamma_i = gamma (i);
290- const Real upscatter = 4 * IC_x0 * gamma_i * gamma_i;
291+ const Real nu_seed = nu / ( 4 * IC_x0 * gamma_i * gamma_i) ;
291292 const Real Ndgamma = dN_e (i);
292293
293- if (nu > upscatter * nu0.back ())
294+ if (nu_seed > nu0.back ())
294295 break ;
295296
296297 bool extrapolate = true ;
297298 for (size_t j = nu_size; j-- > 0 ;) {
298299 const Real nu0_j = nu0 (j);
299300
300301 if (IC_tab (i, j) < 0 ) { // integral at (gamma(i), nu(j)) has not been evaluated
301- Real nu_comv = gamma_i * nu0_j;
302- Real inv = 1 / nu_comv;
302+ const Real nu_comv = gamma_i * nu0_j;
303+ const Real inv = 1 / nu_comv;
303304
304305 const Real grid_value = I_nu_dnu (j) * sigma (nu_comv) * inv * inv;
305306 IC_tab (i, j) = (j != nu_size - 1 ) ? (IC_tab (i, j + 1 ) + grid_value) : grid_value;
306307 }
307308
308- if (upscatter * nu0_j < nu ) {
309+ if (nu0_j < nu_seed ) {
309310 IC_I_nu += Ndgamma * (IC_tab (i, j + 1 ) + (IC_tab (i, j) - IC_tab (i, j + 1 )) / (nu0 (j + 1 ) - nu0 (j)) *
310- (nu0 (j + 1 ) - nu / upscatter ));
311+ (nu0 (j + 1 ) - nu_seed ));
311312
312313 extrapolate = false ;
313314 break ;
314315 }
315316 }
316317
317318 if (extrapolate) {
318- IC_I_nu += Ndgamma *
319- (IC_tab (i, 0 ) + (IC_tab (i, 0 ) - IC_tab (i, 1 )) / (nu0 (1 ) - nu0 (0 )) * (nu0 (0 ) - nu / upscatter ));
319+ IC_I_nu +=
320+ Ndgamma * (IC_tab (i, 0 ) + (IC_tab (i, 0 ) - IC_tab (i, 1 )) / (nu0 (1 ) - nu0 (0 )) * (nu0 (0 ) - nu_seed ));
320321 }
321322 }
322323
0 commit comments