A method of determining an object's initial orbit about an astronomical body, comprising: from known observation positions o 1, o 2,... o n, measuring n angle-only observation directions d o1, d o2,... d on of the object at times t 1, t 2,.... t n where n is three or more; characterised in that the method comprises: generating a search grid in r i and Δ r = r j - r i with a sampling resolution where r i and r j are a range-to-object of the different i th and j th observations; for different combinations of r i and Δ r, solving a boundary-value problem as a function of ( o i, o j), ( d oi, d oj) ( r i, Δ r) to determine a state vector of a three-dimensional position and a velocity of the object at time ti; propagating the position and velocity of the object to estimate observation directions d e1, d e2,... d en at times t k for k = 1 to n, k not equal i or j; and computing an error metric between the measured and estimated observation directions at times t k; identifying the ( r i, Δ r) pair on the search grid at which the error metric represents a best fit to the measured observation directions; and reporting an orbital element set corresponding to the identified ( r i, Δ r) pair that describes the object's initial orbit. from known observation positions o 1, o 2,... o n, measuring n angle-only observation directions d o1, d o2,... d on of the object at times t 1, t 2,.... t n where n is three or more; characterised in that the method comprises: generating a search grid in r i and Δ r = r j - r i with a sampling resolution where r i and r j are a range-to-object of the different i th and j th observations; for different combinations of r i and Δ r, solving a boundary-value problem as a function of ( o i, o j), ( d oi, d oj) ( r i, Δ r) to determine a state vector of a three-dimensional position and a velocity of the object at time ti; propagating the position and velocity of the object to estimate observation directions d e1, d e2,... d en at times t k for k = 1 to n, k not equal i or j; and computing an error metric between the measured and estimated observation directions at times t k; identifying the ( r i, Δ r) pair on the search grid at which the error metric represents a best fit to the measured observation directions; and reporting an orbital element set corresponding to the identified ( r i, Δ r) pair that describes the object's initial orbit. generating a search grid in r i and Δ r = r j - r i with a sampling resolution where r i and r j are a range-to-object of the different i th and j th observations; for different combinations of r i and Δ r, solving a boundary-value problem as a function of ( o i, o j), ( d oi, d oj) ( r i, Δ r) to determine a state vector of a three-dimensional position and a velocity of the object at time ti; propagating the position and velocity of the object to estimate observation directions d e1, d e2,... d en at times t k for k = 1 to n, k not equal i or j; and computing an error metric between the measured and estimated observation directions at times t k; solving a boundary-value problem as a function of ( o i, o j), ( d oi, d oj) ( r i, Δ r) to determine a state vector of a three-dimensional position and a velocity of the object at time ti; propagating the position and velocity of the object to estimate observation directions d e1, d e2,... d en at times t k for k = 1 to n, k not equal i or j; and computing an error metric between the measured and estimated observation directions at times t k; identifying the ( r i, Δ r) pair on the search grid at which the error metric represents a best fit to the measured observation directions; and reporting an orbital element set corresponding to the identified ( r i, Δ r) pair that describes the object's initial orbit.