Summary: New methods for the fast, accurate and efficient calculation of large classes of seismic rays joining two points \(x_ S\) and \(x_ R\) in very general two-dimensional configurations are presented. The medium is piecewise homogeneous with arbitrary interfaces separating regions of different elastic properties (i.e., differing wave speeds \(c_ P\) and \(c_ S)\). In general there are \(2^{N+1}\) rays joining \(x_ S\) to \(x_ R\) while making contact with N interfaces. Our methods find essentially all such rays for a given N by using continuation or homotopy methods on the wave speeds to solve the ray equations determined by Snell's law. In addition travel times, ray amplitudes and caustic locations are obtained. When several receiver positions \(x_ R^{(j)}\) are to be included, as in a gather, our techniques easily yield all the rays for the entire gather by employing continuation in the receiver location. The applications, mainly to geophysical inverse problems, are reported elsewhere.