Correction for the OccultationCalculator::isObservationViable function

[lrbusquets]

Last week, I noticed some weird behaviour in the gaps of my Doppler observations caused by planetary occultations. Reviewing the code, I found what seemed to be the likely cause: here (line 175) the position of the occulting body is taken at the middle instant between the signal is sent and received, which causes a significant offset when one of the two ends is near the occulting body.

Although we've discussed with @DominicDirkx that the correct way to proceed would be with a convergence loop for the position of the occulting body, I moved forward with the following workaround, which uses a weighted linear average where the weights are determined by the distance from the occulting body to each of the two observation ends (note: I don't know why the first paragraph does not have the correct markdown format):

`bool OccultationCalculator::isObservationViable( const std::vector< Eigen::Vector6d >& linkEndStates, const std::vector< double >& linkEndTimes ) { bool isObservationPossible = 1; Eigen::Vector3d positionOfOccultingBodyAtFirstInstant; Eigen::Vector3d positionOfOccultingBodyAtSecondInstant; Eigen::Vector3d positionOfFirstEnd; Eigen::Vector3d positionOfSecondEnd;

double FirstEndEpoch;
double SecondEndEpoch;
double WeightedAvgEpoch;

double distanceToFirstEnd;
double distanceToSecondEnd;

// Iterate over all sets of entries of input vector for which occultation is to be checked.
for( unsigned int i = 0; i < linkEndIndices_.size( ); i++ )
{
    FirstEndEpoch = linkEndTimes.at( linkEndIndices_.at( i ).first );
    SecondEndEpoch = linkEndTimes.at( linkEndIndices_.at( i ).second );

    positionOfFirstEnd = linkEndStates.at( linkEndIndices_.at( i ).first  ).segment( 0, 3 );
    positionOfSecondEnd = linkEndStates.at( linkEndIndices_.at( i ).second ).segment( 0, 3 );

    // Get position of occulting body
    positionOfOccultingBodyAtFirstInstant = stateFunctionOfOccultingBody_(
                linkEndTimes.at( linkEndIndices_.at( i ).first )).segment( 0, 3 );
    positionOfOccultingBodyAtSecondInstant = stateFunctionOfOccultingBody_(
                linkEndTimes.at( linkEndIndices_.at( i ).second )).segment( 0, 3 );

    distanceToFirstEnd = (positionOfOccultingBodyAtFirstInstant - positionOfFirstEnd).norm();
    distanceToSecondEnd = (positionOfOccultingBodyAtSecondInstant - positionOfSecondEnd).norm();

    WeightedAvgEpoch = FirstEndEpoch * (distanceToSecondEnd / (distanceToFirstEnd + distanceToSecondEnd)) + 
                       SecondEndEpoch * (distanceToFirstEnd / (distanceToFirstEnd + distanceToSecondEnd));

    // Check if observing link end is occulted by body.
    if( mission_geometry::computeShadowFunction(
                positionOfFirstEnd, 0.0,
                stateFunctionOfOccultingBody_( WeightedAvgEpoch ).segment( 0, 3 ),
                radiusOfOccultingBody_,
                positionOfSecondEnd ) < 1.0E-10 )
    {
        isObservationPossible = 0;
        std::cout << "Breaking..." << std::endl;
        break;
    }
}

return isObservationPossible;

}`

This seems to correct the odd behaviour, at least, with one-way Doppler observations. I'll be happy to do more tests and/or push these changes if you wish.

[DominicDirkx]

A PR is (finally) open to address this here: https://github.com/tudat-team/tudat/pull/263/files I've taken your solution of calculating the occulting body evaluation epoch, and have also changed the function to calculate the occultation. It had some poor numerical behaviour in edge cases. A unit test has been added, which is passing for the 1-way case, N-way case needs to be added and then we can merge. Thanks for finding this bug and proposing a fix!

[lrbusquets]

Thank you for the update! I'm happy to know I made my tiny bit to make tudat better :)

[DominicDirkx]

N-way test added, and code is merged! Thanks a lot for finding this issue!