The goal of this study is to develop algorithms and programs to compute wave travel times and directions at the nodes of a given 3D grid. We consider a velocity model that can be described analytically or presented by defining velocity values at the nodes of a 3D regular point grid. In the latter case, the velocity model is smoothed by a spatial low-pass filtering. An algorithm of 3D shooting is developed. The algorithm can be used in case of a point source and more generally for an initial wave front position defined as a parametric surface Phi(theta, phi). In the latter case, the initial points of rays on this surface are not known beforehand. In the case of a source point, the ray direction is determined by the variables theta, phi. We consider equations for partial derivatives with respect to the variables theta,phi together with the differential ray system. The entire system is solved by a fifth-order Runge-Kutta method with step-size and precision control. With known derivatives at ray points we can use Newton's iterative method, which guarantees quadratic convergence. Moreover, with these derivatives we can calculate geometrical spreading and thus the amplitudes. Different types of velocity structures such as homogeneous, gradient, etc. and different types of surface Phi(theta, phi) are considered. The results of numerical experiments are presented.