istat logic is broken; means photons aren't destroyed when istat = P_ERROR_MATOM, P_LOFREQ_FF, P_ADIABATIC

[jhmatthews]

Packets are supposed to be destroyed if kpkt() chooses low-frequency free-free or adiabatic cooling as a destruction process. Currently, it's not clear to me that this ever actually happens in the code.

In the original implementation (see e.g. #66 -- 11 years ago!!) istat was set to P_ADIABATIC in kpkt() and then this caused an exit out of the trans_phot loop fairly early on. Currently, it looks to me like this won't happen -- here's the bit of code, with irrelevant parts removed:

  while (istat == P_INWIND)
  {
    istat = translate (w, &pp, tau_scat, &tau, &current_nres);

    if (istat == P_ERROR)
    {
      Error ("trans_phot: abnormal return from translate on photon %d\n", p->np);
      break;
    }

    if (pp.w < weight_min)
    {
      pp.istat = P_ABSORB;
      stuff_phot (&pp, p);
      break;
    }

    if (istat == P_SCAT)
    {
      /* The photon has scattered, as either a resonance or continuum scatter. */
      pp.grid = n_grid = where_in_grid (wmain[pp.grid].ndom, pp.x);

      if ((ierr = scatter (&pp, &current_nres, &nnscat)))       // pp is modified
      {
        Error ("trans_phot_single: photon %d returned error code %d whilst scattering \n", pp.np, ierr);
      }

      /* Reinitialize parameters for the scattered photon so it can can continue through the wind
       */
      tau = 0;
      tau_scat = -log (1. - random_number (0.0, 1.0));
      pp.istat = P_INWIND;
      pp.ds = 0;

      stuff_phot (&pp, p);
      istat = p->istat;
    }

    if (pp.nscat == MAXSCAT)
    {
      pp.istat = P_TOO_MANY_SCATTERS;
      stuff_phot (&pp, p);
      break;
    }

    if (pp.istat == P_ERROR_MATOM || pp.istat == P_LOFREQ_FF || pp.istat == P_ADIABATIC)
    {
      p->istat = pp.istat;
      stuff_phot (&pp, p);
      break;
    }

    /* This is an insurance policy but it is not obvious that, for example nscat
     * and nrscat, need to be updated */

    p->istat = istat;
    p->nscat = pp.nscat;
    p->nrscat = pp.nrscat;
    p->w = pp.w;
  }

What happens here is that pp.istat gets updated by scatter, and while this used to cause the loop to be broken out of, in this case pp.istat gets set to P_INWIND before the relevant if statement. I'm finding it quite hard to follow what's intended here. I think I know how to fix the error -- the if statement involving P_ERROR_MATOM and so on should be moved up to before istat is set to P_INWIND -- but there is plenty of potential for confusion here.

Note this means we haven't been destroying k-packets properly with adiabatic cooling for years, which is a bit worrying given what we write in papers. This must have happened during some reorganisation of trans_phot which changed the behaviour from what was originally intended.

[jhmatthews]

Fixed by PR above.