Another problem with EquatorialCylindricalShape
I think I discovered another problem with the EquatorialCylindricalShape class in the intersectSurface method. The main algorithm in intersectSurface is trying to walk along a segment (d0 to dm) of the vector beginning at the spacecraft and in the direction of lookB (the spacecraft pointing direction). It walks by an increment defined by the variable dalpha. There appears to be a unit conversion problem in the dalpha computation. In particular,
double dalpha = MAX(cos(g1lat * DEG2RAD), cmin) / (2.0 * demScale() * DEG2RAD);
and later in the code
dalpha = MAX(cos(g2lat * DEG2RAD), cmin) / (2.0 * demScale() * DEG2RAD);
The problem is that demScale() returns units of pixels per degree and DEG2RAD is in units of radians per degree so demscale() * DEG2RAD results in pixel radians per degree squared.
The correct computation for dalpha should be
dalpha = MAX(cos(g2lat * DEG2RAD), cmin) / (2.0 * demScale() / DEG2RAD);
This will result in the desired units which are pixels per radian.
This error was discovered at LROC using a high resolution SOCET DEM on a high emission angle image (oblique view). In this case the algorithm was failing inside craters or when there was an abrupt change in the radius of the DTM. Programs like phocube, campt, and qview can not compute the lat/lon/phase/ema/inc angles. We made the change to the code as described and it now correctly computes those values. Not sure how it will impact other instruments such as Dawn. Certainly it needs to be tested to make sure it doesn't fail with too many iterations and see how it might improve intersections at the limb. Please contact me for test data if you need it. I can not upload the cubes as they exceed the size limits of 29MB.