sirocco-rt / sirocco

This is the repository for Sirocco, the radiative transfer code used to winds in AGN and other syatems
GNU General Public License v3.0
27 stars 24 forks source link

error in frame conversion in update_flux_estimators: affects rad_force_es #1030

Closed jhmatthews closed 8 months ago

jhmatthews commented 1 year ago

Currently in update_flux_estimators() we store the radiation force and flux estimators as "observer frame quantities", but we store volume and ne as "CMF quantities". Currently the volume and density are transformed in the same way, which must be wrong, and I believe it is the density which is incorrect, i.e. n_obs should be gamma * n_cmf.

  electron_density_obs = xplasma->ne / wmain[nwind].xgamma_cen; // Mihalas & Mihalas p146
  volume_obs = wmain[nwind].vol / wmain[nwind].xgamma_cen;

  for (i = 0; i < 4; i++)
  {
    xplasma->rad_force_es[i] *= (volume_obs * electron_density_obs) / (volume_obs * VLIGHT);
    xplasma->F_vis[i] /= volume_obs;
    xplasma->F_UV[i] /= volume_obs;
    xplasma->F_Xray[i] /= volume_obs;
  }

  for (i = 0; i < NFLUX_ANGLES; i++)
  {
    xplasma->F_UV_ang_x[i] /= volume_obs;
    xplasma->F_UV_ang_y[i] /= volume_obs;
    xplasma->F_UV_ang_z[i] /= volume_obs;
  }

Do you know if rad_force_es is used anywhere in the hydro work @nscepi ?

Also note that rad_force_es has a volume_obs in both the numerator and the denominator ??

nscepi commented 1 year ago

Nick used to use the force estimator from PYTHON but at some point he switched to only using the flux estimator, which we then use to reconstruct a force in PLUTO. The machinery to use the force is still here, you just need to set the flag PY_RAD_DRIV to ACCELERATIONS instead of FLUXES.

jhmatthews commented 1 year ago

Thanks! so I guess it doesn't affect any work in progress, but should still be fixed of course.