I've been using the VSOP87 theory to calculate transits, such as transits of Venus as seen from Earth. To detect a transit at the closest approach of a sun and planet, the method I use is to compute the sums of the radii of the sun and planet, and if it is less than their separation then a transit occurs. In this discussion for the sake of simplicity I assume an observer on the Earth but the method that I am using is good for any transit involving two terrestrial planets during the time range where the VSOP87 theory is sufficiently accurate.
So far it is working well but at present my code is not able to detect a grazing transit. A grazing transit is one that cannot be observed by a geocentric observer but may be observed by an observer on the earth's surface that is closer to the transit centreline.
The version of the VSOP87 theory that I am using provides the heliocentric coordinates of the planets, so the sun is at the origin. Given the geocentric observer O, the sun S and the planet P, the VSOP87 theory provides the vectors SO> and SP>. (Please excuse my adhoc vector notation.)
The vectors OS> = -SO> and OP> = SP> - SO> form the basis of the transit calculation. The actual calculation uses the coordinates but at its core the algorithm that I use (from "Astronomical Algorithms" by Meeus) converts these to vectors (represented as rectangular xyz coordinates in code).
To calculate a grazing transit, the idea that I had is to find an additional point to represent the observer that is closest to the centreline for the transit. Call this closest point C. The vector that is derived from this point is CO>. Then we can find the vectors CS> = CO> + OS> and CP> = CO> + OP> and perform the transit calculation again using CS> and CP> instead of OS> and OP>. Doing it this way is intended to take the parallax of the sun into account; this is not negligible.
Where I'm a little stuck is the necessary mathematics to find the vector CO>. It's been many years since I studied vectors and my mathematics here is a little deficient. So far I have determined that the magnitude of the vector is the radius of the earth (I am neglecting oblateness here for the time being) and its direction is the projection of the vector SP> onto the plane that is perpendicular to the vector SO>.
It's also possible that the method itself is deficient or someone else can suggest an easier way, which is why I have posted it here.