Here are two old scripts that find individual harmonics algebraically:
Least -squares fit time series to C + Asin(wt) + Bcos(wt), where C is the mean, w is omega.
The first one (harm-6b.go) works for regular data with no blanks (then C is the simple time average).
I think you can prove that the result of the least-squares fit is exactly the same as that of a Fourier transform for a regularly-spaced data.
The least-squares fit is a useful method when there are missing data or when the data interval isn't regular.
Cheers,
Ryo