These formulations of the post-Newtonian equations for the phase and frequency are simple to implement: each pass through the loop increments the time by the sample time (Sample_Time in the example) and computes the phase and frequency using Eqs. () and (). However, there is an alternative formulation. In deriving these equations the ``natural'' equation that arises is of the form . [See e.g.  Eq.(3).] This in turn can be integrated to give an equation of the form . In our formulation this equation has been inverted - throwing away higher-order post-Newtonian terms as you go - to give Eq.(). However the equation in the form can also be implemented directly. In this type of formulation one would again increment the time, but then use a root-finding routine to find the frequency at each time step. Our chosen method has the advantage of avoiding a time-consuming root-finder at each time step; however the alternative formulation has undergone fewer damaging post-Newtonian transformations, and may therefore be more accurate.
In our formulation we only need to call a root-finding routine at the start of the chirp to find the value of when the system is at the initial frequency. In order to insure that we find the correct root for the starting time we begin a search at a time when the leading order prediction of the frequency is well below the desired starting frequency. We step forward in time until we bracket the root; we then call the Numerical Recipes root-finder rtbis() to compute the root precisely. This is depicted in the lower right corner of figure where we show the value of the ``time'' coordinate that corresponds to an initial frequency of Hz. This method is virtually assured of finding the correct root in that it will find the first solution as we proceed from right to left in figure . The primary problem in finding this root is that there may actually be no meaningful start-time for the specified chirp. For example, if you you were to specify a chirp with two objects with an initial frequency of 1000Hz, you can see from the figure that there is no value of (i.e. ) that corresponds to this frequency. In this case phase_frequency() will search from right to left for the start time. It will notice that it is passing over the peak in the graph and out of the regime of post-Newtonian viability. It will then terminate the search and notify the caller that there is no solution for the requested chirp.
The behavior of the frequency equation is shown in figure . As time increases the frequency rises to a maximum and then begins to decrease dramatically. Notice that the maximum occurs when the dimensionless time parameter is approximately unity; this feature is only weakly dependent on the mass ratio. The fact that means the post-Newtonian corrections in Eq.() are comparable to the leading order term. Therefore, this peak is a natural place to terminate the post-Newtonian chirp approximation. In the example the code terminated the chirp for precisely this reason. [See the warning message.]
Although it is not shown in the figure the behavior of as nears zero is very abrupt; the function goes sharply negative and then turns around and diverges to as (i.e. ). This abrupt behavior will happen on a time scale of order (a few microseconds). Typical sample times are likely to be on the order of a tenth of a millisecond, and therefore the iterative loop may step right over this maximum-minimum-divergence behavior of the frequency function altogether. Don't worry. The routine phase_frequency() handles this case gracefully. The routine will stop the chirp calculation and warn the caller if the time stepper goes beyond the coalescence time. It will also stop the chirp calculation if it senses that the time has stepped over the dip in frequency and is on the strongly divergent part of the frequency curve near the axis.