Toroidal magnetic configurations are widely exploited in industry and scientific research, involving a vast spectrum of applications, such as thermonuclear fusion, particle detectors, SMES systems and medical devices. To properly design and analyse these sys-tems, it is crucial to determine the magnetic field generated by dif-ferent configurations. The multipole expansion theory can be ap-plied to the analysis of toroidal configurations, by solving the La-place equation for the magnetic scalar potential in toroidal coor-dinates. Contrarily to the case of accelerator magnets with straight axis, in this case the correlation between the current dis-tribution and the field harmonics cannot easily be identified. This paper proposes a methodology for the computation of field har-monics in toroidal coordinates, which is validated by comparison with the results obtained through the Biot-Savart law. This work was carried out in the frame of the GaToroid project ongoing at CERN.