Various authors have presented the aberration function of an optical system as a power series expansion with respect to the ray coordinates in the exit pupil and the coordinates of the intersection point with the image field of the optical system. In practical applications, for reasons of efficiency and accuracy, an expansion with the aid of orthogonal polynomials is preferred for which, since the 1980s, orthogonal Zernike polynomials have become the reference. In the literature, some conversion schemes of power series coefficients to coefficients for the corresponding Zernike polynomial expansion have been given. In this paper we present an analytic solution for the conversion problem from a power series expansion in three or four dimensions to a double Zernike polynomial expansion. The solution pertains to a general optical system with four independent pupil and field coordinates and to a system with rotational symmetry in which case three independent coordinate combinations have to be considered. The conversion of the coefficients is analytically in closed form and the result is independent of a specific sampling scheme or sampling density as this is the case for the commonly used least squares fitting techniques. Computation schemes are given that allow the evaluation of coefficients of arbitrarily high order in pupil and field coordinates.