The method of Väisälä can be extremely useful in many cases. It gets a lot of use when you have a very short arc of observation, not enough to determine a "real" orbit, and just want to get an orbit good enough to predict positions for the next week or so. It can also provide a good starting orbit for an iterative orbital improvement scheme.
I've found very little documentation of this method, and I've also come up with what I consider to be two nice modifications for implementing it. They are used for the Väisälä option in my Find_Orb orbit determination software. Thus, these comments.
"Normally", the method of Väisälä assumes you have two observations of an object, at times t1 and t2. A Väisälä orbit is one that satisfies both observations exactly, and has the object at an apsis at time t1 (perihelion or aphelion, normally). Mathematically, this means that dr/dt = 0 at t=t1, where r=distance between the sun and the object.
The first way in which I've modified this is to "re-define" this sort of orbit to be one in which, at times t1 and t2, the object is at a specific distance from the Sun. As a consequence, there must be an apsis at the mid-time, (t1 + t2) / 2. If t1 and t2 are close, the distance at apsis will be close to the selected distance.
The reason I did this is that it simplifies the code immensely. You know the RA/dec to the object at times t1 and t2, and you know it was at a particular distance from the sun on those dates (the same distance on both dates). Given that, determining the distances of the object from the observer and the (x, y, z) Cartesian coordinates for the objects at times t1 and t2 is quite simple. And from there, getting an orbit is straightforward. Many orbital determination systems have, somewhere near their code, a function that tackles the problem "the object was at (x1, y1, z1) at time t1 and at (x2, y2, z2) at time t2; find the orbit specified by this". That's the only function you need. (In my case, it was required for the method of Herget, and was therefore already in place.)
The second modification is a "linearization" trick. It addresses the fact that Väisälä orbital methods are usually limited to inclusion of only two observations. You generally have at least four observations to play with after two nights, and it would seem logical to use all the data. (As will be noted below, this trick also can improve Herget orbits.)
My first thought was to proceed as follows: Compute the "usual" sort of Väisälä orbit using the first and last observations. Use this as a starting point in doing a least-squares fit in which the orbit is constrained to a given apsis distance and date (q and Tperihelion). Since a general-solution orbit has six free parameters, the above two constraints leave four parameters to solve for. (If you're using "classical" orbital elements, these would be the eccentricity, argument of perihelion ω, inclination i, longitude of the ascending node Ω.)
The above would have worked quite well. However, I found a method which was easier to implement and less susceptible to instabilities and mathematical singularities. It's also faster. It works as follows. Suppose you've computed the "usual" two-point Väisälä orbit. Your residuals at times t1 and t2 are zero; you have non-zero residuals everywhere in between.
Consider the residuals in, say, RA. If you do a linear fit to them, you would get a function such as
delta_RA = A + B * t
You can do another linear fit to get a similar function in declination. Furthermore, such fits can be done analytically, in a single step, and are not at all computationally demanding.
The core of my second modification was this. Assume you've computed the above A and B coefficients, as well as C and D coefficients for
delta_dec = C + D * t
Now modify your first and last observations, (ra1, dec1) and (ra2, dec2), as follows:
ra1new = ra1 - (A + B * t1)
dec1new = dec1 - (C + D * t1)
ra2new = ra2 - (A + B * t2)
dec2new = dec2 - (C + D * t2)
Now re-compute the Väisälä orbit, except this time, instead of computing an orbit connecting (ra1, dec1) to (ra2, dec2), compute one connecting the "new" values.
What you've done here is a first-order correction: you've computed the average residuals and the trends in those residuals, and adjusted the first and last observations to get an orbit in which those trends are canceled out and the average residual is zero. The result is very similar to what you would have gotten with the four-parameter fit, but requires a fraction of the math and computation time.
Another nice aspect of this trick is that it can be used with Herget orbits, which also assume zero residuals for two observations (again, usually the first and last). Apply this linearization method, and the errors become distributed among all observations.
Sonia Keys has pointed out one slight flaw in this scheme: it breaks down as one approaches a celestial pole. If one end is, say, half as far from a celestial pole as is the other, then a second of RA is about twice as big at one end as it is at the other.
Usually, this is not much of a problem. This method is normally used for short arcs only. If the arc is, say, five degrees long, then you'd have to be near declination +/- 80 degrees to see significant damage from this effect. But as Sonia points out, one can get around the issue almost completely by rotating to a different coordinate system, one in which the end observations have "declinations" of zero. Then do the above linearization trick, and rotate back. Sonia has implemented this; I've not bothered.