Given the coordinates of A, B and C, and knowing that the angle in D is straight, we do not need to know the distance between B and D. We could calculate it, but in truth it is not necessary to know the coordinates of D.
Consider the straight spherical triangle formed by the vertices A, B, and D, with right angle in D. Let
alpha be the angle in A. Knowing the distance between A and B (let's call it
d), the distance between A and D (let's call it
b) can be calculated as:
b = ATAN( TAN( d) * COS( alpha))
This formula is developed on the unit sphere. We will do all the calculations on it.
b, will always be included in the first quadrant for your use case (they will be in the interval between
2*k*PI multiples of that interval), the results should always be positive and they would not have ambiguities.
d can be calculated by developing the haversine formula:
d = 2 * ASIN( SQRT( SIN(( phi_b - phi_a) / 2)^2 + COS( phi_a) * COS( phi_b) * SIN(( lambda_b - lambda_a) / 2)^2))
lambda_b are the latitudes and longitudes of points A and B respectively, all in radians.
alpha is the absolute difference between the initial azimuth from A to B and the initial azimuth from A to C:
alpha = ABS( azimuth_ab - azimuth_ac)
And both azimuths can be calculated by developing the spherical law of sines and the spherical law of cosines for a spherical triangle with a vertex at the North pole (allow me in this opportunity to omit the full development):
azimuth_ab = ATAN2( SIN( lamda_b - lambda_a) * COS( phi_b), COS( phi_a) * SIN( phi_b) - SIN( phi_a) * COS( phi_b) * COS( lambda_b - lambda_a)) azimuth_ac = ATAN2( SIN( lamda_c - lambda_a) * COS( phi_c), COS( phi_a) * SIN( phi_c) - SIN( phi_a) * COS( phi_c) * COS( lambda_c - lambda_a))
Note: Sines and cosines can be positive or negative depending of the global position of the points, so we use the ATAN2() function instead of the ATAN() for a division.
azimuth_ac, we can find the coordinates of D. The development of these formulas was already exposed in the following answer: Get lat/long given current point, distance and bearing
phi_d = ASIN( SIN( phi_a) * COS( b) + COS( phi_a) * SIN( b) * COS( azimuth_ac)) lambda_d = lambda_a + ATAN( SIN( azimuth_ac) * SIN( b) * COS( phi_a) / (COS( b) - SIN( phi_a) * SIN( phi_d)))
Other notes and some reference links:
dis the length of the side opposite the angle at D in the spherical triangle, in the same way that the distance
bis the length of the side opposite the angle at B. Both lengths correspond to the angle of the arc segment, in radians, on the unit sphere. At no time is it necessary to use a certain radius for the sphere, but if you wish to know lengths on a sphere of radius
1, you must multiply the distances calculated by that radius.
I know that there are libraries for solving the inverse problem that allow us to calculate the distance between two points of known coordinates on the ellipsoid more precisely than the formula of the haversine (although probably it would not make sense to use them to reduce then that distance to the unit sphere), through iterative methods or reductions of series to their first terms. I'm not really familiar with them but if links are proposed in the comments I can add them to my answer.
Many of these formulas can be found on the Ed Williams' Aviation Formulary
Whuber and cffk have written great answers about the error produced by calculations on the sphere instead of on the ellipsoid: How accurate is approximating the Earth as a sphere?.
Sines of very small angles are very small too (and cosines are very close to
1). Use as many decimal places as you can.
Please verify the results with known coordinate data before making productive use of these formulas.