newtonian mechanics - How to calculate a ballistic trajectory for a suborbital flight?

I'm trying to calculate the initial launch angles and velocity of a projectile (atmosphere's effects can be neglected), assuming that I know the lat/lon coordinates of both the launch and the destination points. Can you recommend any software libraries or reading materials that could help? Thank you.

Answer

If you don't consider atmospheric drag and other dissipative effects, assuming that the Earth is non-rotating, and ignore general relativity (which is a lot of assumptions to make, and later I can try to edit this to comment on how they might be relaxed), a ballistic trajectory is just a segment of an orbit. So finding the launch parameters just entails finding an orbit that intersects both the launch point and the destination.

Orbital trajectories are given by the equation

$$r_o = \frac{\rho}{1 + \epsilon\cos\phi}$$

where $\rho$ is a length parameter related to the size of the orbit and $\epsilon$ is the eccentricity. At the points where the orbit intersects the Earth's surface, you'll have $r_o$ equaling the radius of the Earth at that point; let's say $r_o = R_L$ for the launch point and $r_o = R_D$ for the destination. This gives you the following two equations,

$$R_L = \frac{\rho}{1 + \epsilon\cos\phi_L}$$

$$R_D = \frac{\rho}{1 + \epsilon\cos\phi_D}$$

along with the constraint that $\phi_D - \phi_L$ (or vice-versa, depending on how you define coordinates) has to be equal to the angular separation $\Delta$ between the source and destination. That angular separation can be calculated using the law of haversines.

This system of equations can be solved for the orbital parameters as follows, assuming $R_D \neq R_L$ (see below for the other case),

$$\epsilon = \frac{R_L - R_D}{R_D\cos(\phi_L + \Delta) - R_L\cos\phi_L}$$

$$\rho = \frac{R_D R_L [\cos(\phi_L + \Delta) - \cos\phi_L]}{R_D\cos(\phi_L + \Delta) - R_L\cos\phi_L}$$

Notice that they are dependent on $\phi_L$; this is because only the difference $\phi_D - \phi_L = \Delta$ is physically relevant, so you're free to choose the actual values to be whatever you like. Your choice of $\phi_L$ will influence the shape of the trajectory.

Having calculated $\epsilon$ and $\rho$, you can determine the launch angle $\alpha$ above the horizon by calculating the slope of the orbit at the launch point:

$$\tan\alpha = \left.\frac{1}{r_o}\frac{\mathrm{d}r_o}{\mathrm{d}\phi}\right|_{\phi_L} = \frac{\epsilon\sin\phi_L}{1 + \epsilon\cos\phi_L} = \frac{(R_L - R_D)\sin\phi_L}{R_D[\cos(\phi_L + \Delta) - \cos\phi_L]}$$

To determine the speed, you can use the fact that (according to my notes) the total energy of the projectile is given by

$$E = \frac{GMm(\epsilon^2 - 1)}{2\rho}$$

which is equal to the sum of kinetic and potential energies, $\frac{1}{2}mv^2 - \frac{GMm}{r}$. ($m$ is the mass of the projectile, $M$ is that of the Earth) Plugging in $r = R_L$, I get

$$v = \sqrt{\frac{2GM}{R_L} + \frac{GM(\epsilon^2 - 1)}{\rho}}$$

So the bottom line is that you plug $R_D$, $R_L$, $\Delta$, and some choice of $\phi_L$ into the formulas for $v$ and $\alpha$ to get the launch parameters.

---

I mentioned that the procedure above hits a snag if $R_D = R_L$. You'd wind up calculating $\epsilon = 0$, which on a spherical planet entails rolling your projectile along the ground ;-) which doesn't make sense.

If $R_D = R_L$, you can go back to the orbital equations and find that $\cos\phi_D = \cos\phi_L$. (Alternatively, this could come from the condition that the denominator of $\epsilon$ also be zero.) The only way to satisfy this is by setting $\phi_L = -\phi_D = \Delta/2$. This gives you the condition

$$\rho = R_L\left(1 + \epsilon\cos\frac{\Delta}{2}\right)$$

Again, you have a degree of freedom: you can choose any value of $\epsilon$, and plugging into this equation will give you the corresponding value of $\rho$. Once you have those, you can use the same procedures to calculate the launch angle and speed:

$$\tan\alpha = \left.\frac{1}{r_o}\frac{\mathrm{d}r_o}{\mathrm{d}\phi}\right|_{\phi_L} = \frac{\epsilon\sin\frac{\Delta}{2}}{1 + \epsilon\cos\frac{\Delta}{2}}$$

$$v = \sqrt{\frac{2GM}{R_L} + \frac{GM(\epsilon^2 - 1)}{R_L(1 + \epsilon\cos\frac{\Delta}{2})}}$$