Jigsaw's solution for target body parameters (triaxial shape and mean radius) are inconsistent with radii triangulated during the bundle

[blandoplanet]

ISIS version(s) affected: 3.10.* (haven't tried older/newer versions)

Description
Brief Summary:
When we solve for target body parameters (in this cases for a global Europa control network) in jigsaw the resulting triaxial shape (a, b, c radii) and the mean radius are inconsistent with the triangulated radii of the tie points derived during bundle adjustment of the same control network.

Specifics:
Before solving for target body, we bundle adjust the final Europa control network (Galileo and Voyager images) including solving for radius. The result is a set of lat/lon/radius values for each of our 53,000+ points. See /work/projects/EuropaBasemap/lweller/GLL_VGR_Network/TargetBody_FinalSolutions/2020EUROPA_PM/RadAngTwist_GlobalMerged_2020_CilixFree_NewPM_bundleout_points.csv
This bundle has an updated prime meridian and the images have updated (though not quit final) SPICE. The image list can be found here: /work/projects/EuropaBasemap/lweller/GLL_VGR_Network/TargetBody_FinalSolutions/2020EUROPA_PM/
NewPM_GalileoVoyager_Europa_Merged_2020.lis

These radii are "good." A quick and dirty dem shows that the point cloud forms a triaxial shape with the long axis near 180 deg, and the short axis aligned with the pole. Here I just shoved the point cloud through ASP point2dem. Nothing fancy. The cloud is very sparse in some areas, thus a lot of no-data. The point is simply that the cloud is triaxial (yellows high, green/light blue intermediate, dark blue low).

An interrogation of the points themselves indicates that they are consistent with IAU values. Only ~30 points out of 53,000+ have corrections larger than 500 m. The largest point values cluster around 1562.6 km, which is the IAU/pck a-axis and are all near 180 deg. The smallest values cluster at 1559.5 km (again consistent with the IAU/pck) and are all at high latitudes, as expected. If we just take the mean (a terrible way to calculate the mean radius) we get 1560.99 km, which is pretty darn close to the IAU mean radius of 1560.8 km - especially considering a lot of our points happen to be near the longest axis (very few points near the poles). The largest radii is 1563.37 km, and this appears to be a bit of an outlier (the next largest is 1562.97 km). So a quick look at the data suggests the points are well behaved and consistent with IAU values. Finally the difference between a-axis and c-axis (which can be interpreted geophysically) is about 3.1 km.

Now, we take the same network and use it to solve for the triaxial shape (one jigsaw run) and the mean radius (a second jigsaw run). These solutions can be found here: /work/projects/EuropaBasemap/lweller/GLL_VGR_Network/TargetBody_FinalSolutions/2020EUROPA_TRIRADIUS/
The results are radii that are significantly larger than IAU values: e.g., 1563.9 km, 1562.8 km, and 1561.6 km for a, b, and c respectively (using no constraints). Note that in the target body solution we can't actually solve for radii to have a direct comparison. A few red flags. 1) the value for the a-axis is larger than even the largest outlier from the points themselves (1563.9 vs 1563.4 km). 2) The mean radius of 1562.7 km is 1.7 km larger than the simple average of the points, and similar to the largest values in our cloud (i.e., the a-axis). Recall that our "average" is already biased toward points on the long axis. 3) the short axis of 1561.6 is 0.8 km longer then where most of the points on the c-axis seem to be. 4) a-c is only 2.3 km. So it isn't just a shift from the IAU, the relative length of the axis has changed.

@lwellerastro spent a lot of effort assessing whether the sigma values used made any difference. They do tweak the result, but not in a meaningful way. These results can be found in the dir mentioned above. Each "sigma" value has its own network and bundleout (labeling should be clear). The pvl files are also in that directory.

Additional context
We used jigsaw to solve for target body parameters for Enceladus following an identical method. In that case the triaxial solution fall right in the radii point cloud from the bundle solution (I just double checked!). The values we found are consistent with the IAU values, and the mean radius we derived was reasonable.

Possible Solution
It's not entirely clear that this is a bug. I'll note that there is essentially no documentation for this application (see #4111) so it's possible that it is pilot error. However, at this point that seems unlikely to me. We've used this successfully before. It's also possible that the point cloud is just too sparse for whatever method is being used. If that's the case, the uncertainty values reported should be much larger (again, it's not clear how the uncertainty is determined). It's also possible that this is a bug (introduced between isis3.10 and isis3.8.0)

[lwellerastro]

We used isis3.8.1 for Enceladus.

[jessemapel]

There haven't been significant changes to Jigsaw between 3.8.1 and 3.10 so I doubt it is a bug. The target body code in particular hasn't been touched in many years.

Thank you for the extensive documentation and explanation on this. It will likely take some detective work to figure this one out.

[blandoplanet]

Some historical context from Brent Archinal (may or may not be useful):

"The original algorithm was developed by Tim Colvin and Mert Davies at Rand for their software, which after being transferred here became the ISIS2 randlsq program. (Although unfortunately they had not apparently not implemented it at the time of Tim’s publication of the software description at https://www.rand.org/pubs/notes/N3330.html, so details are not included there. The algorithm though basically includes an expansion of the observation equation to solve for a and b, or a, b, and c, instead of just R). That part of the program was (quite unfortunately) not included when randlsq was converted into the ISIS3 jigsaw program, but Ken Edmundson added it in a few years ago (also unfortunately after he had documented Jigsaw with his 2012 ISPRS publication and before he was ever allotted time to update the documentation) based on the randlsq algorithm/code. I do not recall what level of testing he did with the program, so cannot say whether he would have caught this problem or not. I can try to look through my records/e-mail to see what discussion we had about that."

[jessemapel]

@blandoplanet I looked at the document you posted and it doesn't appear to be solving for target radii. It seems to only solve for target rotation parameters, sensor parameters, and point parameters.