Calculating THD using MATLABs fft function

I have a sampled sine wave of a 10kHz signal. It was sampled at 1Msps (1Mhz) so there is 100 samples per period. The total vector is 1024 samples long so its just over 10 periods total.

I have measured the THD of the analog signal using the fft function of a tektronix scope and I know the correct THD to be less than 1%. However, when I try to do this using MATLAB I am not getting the correct results and I expect this is because I am doing it incorrectly. I have pasted my current script below in hopes that somebody would be able to tell me what I am doing wrong. Thanks in advance!

%Use fft to find harmonics
harmonics = abs(real(fft(data))).^2;
%Assume the fundamental frequency is the highest
fundamental_bin = min(find(harmonics==max(harmonics)));
fundamental_power = harmonics(power_bin);
%Find the sum of the rest of the harmonics
harmonics_power = sum((harmonics((2:numHarmonics).*fundamental_bin)));
%Calculate THD as a ratio of square roots of powers
THD = 100*sqrt(harmonics_power)/sqrt(fundamental_power);


When I run this script I get a THD of over 8%, which as I explained earlier is far too high. Can anyone point me in the right direction or give me an example of someone who has done this correctly?

Hi,
If you have found the correction for this code, can you please provide me with it.....

Thank you..