can a great circle route be predicted from initial condition? Assuming the world is a sphere with no wind, can the great circle route of a vessel be predicted from the current position $\{\phi_i,\lambda_i\}$ and the current true course $\theta_i$?
Presently, I'm not concerning myself with speed, just the ability generate a formula for the GC so that I can draw it.
These are the useful formulae I've discovered.
$$ \begin{array}{ll}
d &= \cos^{-1}\left(\sin\phi_i \sin\phi_j + \cos\phi_i \cos\phi_j \cos (\lambda_i-\lambda_j) \right) \\
b &= \cos^{-1}\left(\frac{\sin\phi_j-\sin\phi_i\cos(d)}{\sin(d)\cos\phi_i}\right) \\
s &= \sin(\lambda_j-\lambda_i) \\
\theta_i &= \left\{
\begin{array}{ll}
                b  & \mbox{if } s < 0 \\
                2\pi-b & \mbox{if } s > 0
        \end{array}
\right.
\end{array}
$$
I am not presently concerning myself with the poles.
It is true to say that any great circle will cross all meridians so it does not seem unreasonable to determine a formula of the form:
$$ \phi = f(\lambda,\phi_0,\lambda_0,\theta_0) $$
I note that $\cos(2\pi - \theta) = \cos(\theta)$ so:
$$\begin{array}{ll}
 \cos(\theta_i) &= \frac{\sin\phi_j-\sin\phi_i\cos(d)}{\sin(d)\cos\phi_i} \\
&= \frac{\sin\phi_j-\sin\phi_i(\sin\phi_i \sin\phi_j + \cos\phi_i \cos\phi_j \cos (\lambda_i-\lambda_j))}{\sin\left(\cos^{-1}\left(\sin\phi_i \sin\phi_j + \cos\phi_i \cos\phi_j \cos (\lambda_i-\lambda_j) \right)\right)\cos\phi_i}
\end{array}$$
This is an equation that has all the right variables, but I can't rearrange it.
Is there a better approach I should be taking?
UPDATE
It's worth noting that if we use colatitude ($\phi'$), the distance equation looks just like the law of cosines:
$$cos(d) = \cos\phi'_i \cos\phi'_j + \sin\phi'_i \sin\phi'_j \cos (\lambda_i-\lambda_j)$$
 A: You’re in an angle-side-angle situation, and you want to know the length of a side. One angle is given to you by the original heading, I’ll call it $\eta$, say measured clockwise from north. The other angle is the angle at the pole, difference between the original longitude and the longitude of the meridian you’re thinking of as variable. Let’s call that $\theta$. The side is the “colatitude”, the distance, in degrees (or radians) from the pole to the original position. Let’s call this $\lambda$. And you want to know the colatitude of your new point, since you already know its longitude. Call this $\Lambda$. Then I would solve your problem by using two standards from trigonometry, the alternative form of Law of Cosines, which reads
$$
\cos a=-\cos b\cos c + \sin b \sin c \cos A\,,
$$
where the angles of the triangle are $a,b,c$ and the corresponding opposite sides are $A,B,C$. The other is the Law of Sines, which reads,
$$
\frac{\sin a}{\sin A}=\frac{\sin b}{\sin B}=\frac{\sin c}{\sin C}\,,
$$
very much like the corresponding one in plane trig.
Now, we need to name the angle opposite the side that measures your original colatitude, call it $\alpha$. It becomes the $a$ of the Law of Cosines, and the longitude-difference $\theta$ becomes your $b$ there, while the heading angle is $\eta$, your $c$, and the (original) given colatitude is your $A$. Now apply L of C to get $\alpha$, and plug this into Law of Sines, since you have
$$
\frac{\sin\Lambda}{\sin\eta}=\frac{\sin\lambda}{\sin\alpha}\,.
$$
You wanted $\Lambda$, and there it is.
A: Attempt at @Lubin's solution
Let A be the colatitude of the initial position
Let B be the colatitude of the final position 
Let C be the distance travelled along the GC
Let c be the longitude difference
Let b be the initial course 
Let a be the ill defined angle
We wish to find B given b, c and A
$$\cos a = -\cos b\cos c + \sin b \sin c \cos A$$
$$\frac{\sin a}{\sin A}=\frac{\sin b}{\sin B}=\frac{\sin c}{\sin C}$$
$$
\begin{array}{ll}
\sin B&=\frac{\sin b \sin A}{\sin a} \\
&=\frac{\sin b \sin A}{\sin (\cos^{-1}(-\cos b\cos c + \sin b \sin c \cos A)}\\
B &=\sin^{-1}\left(\frac{\sin b \sin A}{\sin (\cos^{-1}(-\cos b\cos c + \sin b \sin c \cos A)}\right)\\
\end{array}
$$
Now converting to latitudes from colatitudes
$$
B' =\cos^{-1}\left(\frac{\sin b \cos A'}{\sin (\cos^{-1}(-\cos b\cos c + \sin b \sin c \sin A')}\right)
$$
And into lat and lon
$$
\phi =\cos^{-1}\left(\frac{\sin \theta_0 \cos \phi_0}{\sin (\cos^{-1}(-\cos \theta_0 \cos (\lambda - \lambda_0) + \sin \theta_0 \sin (\lambda - \lambda_0) \sin \phi_0)}\right)
$$
And that work other than some obvious multiple possible inverses issues. I guess there is some conditional I can apply.
I wonder if it can be simplified.
There is an oddity for some paths where all the lats are over 90, so I take them from 180.
In fact, it's also not handling negative latitudes at all.
UPDATE
$$B =\sin^{-1}\left(\frac{\sin b \sin A}{\sqrt{1-(-\cos b\cos c + \sin b \sin c \cos A)^2}}\right)$$
Lost an inverse, but got nowhere.
