ds_to_plane inconsistent with other ds routines: results in errors in imported models

[jhmatthews]

In issue #1025 I've fixed one problem, but I'm now finding that some photons launched from underneath the wind don't pass through it in an imported model. Look at the purple photon here:

For imported models, what should happen in this case is that on the first call to translate in space, ds_to_wind should calculate the following distances:

• a) the distance to the edge of the domain defined by rmax
• b) the distance to a plane at the lowest z of the wind defined by zdom[ndom].windplane[0]
• c) the distance to a plane at the highest z the wind defined by zdom[ndom].windplane[1]
• d) the distance to cylinders at rho_min and rho_max.

It should then find whichever of these distances is smallest and return that value. In this case, b) should be the smallest distance, but the code instead finds a). Why does this happen?

First note that the initial photon position is: 8.2082e+09 4.8195e+10 -1.0000e-06 Secondly note this means the photon is going "downwards" i.e. is emerging from below the disk plane - that's why the z position is negative EPSILON (see line 173 of disk_photon_gen.c). Thirdly note that at the start of ds_to_wind we initially set ds to the distance to the edge of the domain and the following bit of code is just for imported models.

Here's the relevant bit of code in ds_to_wind:

      x = ds_to_plane (&zdom[ndom].windplane[0], &ptest);
      if (x > 0 && x < ds)
      {
        stuff_phot (pp, &qtest);
        move_phot (&qtest, x);
        rho = sqrt (qtest.x[0] * qtest.x[0] + qtest.x[1] * qtest.x[1]);
        if (zdom[ndom].wind_rhomin_at_disk <= rho && rho <= zdom[ndom].wind_rhomax_at_disk)
        {
          ds = x;
          *ndom_current = ndom;
          xxxbound = BOUND_ZMIN;
        }
      }
      x = ds_to_plane (&zdom[ndom].windplane[1], &ptest);
      if (x > 0 && x < ds)
      {
        stuff_phot (pp, &qtest);
        move_phot (&qtest, x);
        rho = sqrt (qtest.x[0] * qtest.x[0] + qtest.x[1] * qtest.x[1]);
        if (zdom[ndom].wind_rhomin_at_disk <= rho && rho <= zdom[ndom].wind_rhomax_at_disk)
        {
          ds = x;
          *ndom_current = ndom;
          xxxbound = BOUND_ZMAX;
        }
      }

I've discovered that ds_to_plane returns a negative number in this case. This means the x>0 check is not passed and so the ds is not set to this number.

Note in other cases (ds_to_cone, ds_to_cylinder, ds_to_sphere) we seem to only return a positive number from a root finding algorithm, so ds_to_plane seems to have inconsistent behaviour with these. I could just ensure it returns a positive number, but this could break other places where this is used. Additionally, I don't yet know whether the negative number would decrease in magnitude if I moved the photon a little bit further from the disk in the bottom hemisphere.

I think we should probably add a check function to Ed's unit tests that checks all dsto* routines.

[jhmatthews]

I've now checked and moving the photon in the southern hemisphere by ds increases the distance to the windplane. This means, I think, that the ds_to_plane routine is calculating the 3D distance to a plane in the (+x,+z) quadrant without putting the photon in the same quadrant.

[jhmatthews]

Fixed by #1029 with unit test also written.