Calculate Total Harmonic Distortion (THD) from FFT Plot

I am new to FFT and currently studying it for harmonic analysis. I have a set of data of inrush currents from the results of simulation involving equipment energization. The data has $$N = 40,001$$ sample points evenly spaced at $$t = 50 \mu_s$$. I generated the FFT spectrum through a python code below.

import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
from numpy.fft import fft, ifft, fftfreq, rfftfreq, rfft

df = pd.read_csv('inrush2.csv', delimiter = ',')

plt.plot(df['Time'], df['Phase A'])
plt.plot(df['Time'], df['Phase B'])
plt.plot(df['Time'], df['Phase C'])
plt.xlabel('Time (s)')
plt.ylabel('Current (kA)')
plt.show()

# Sampling rate
sr = np.size(df['Phase A'])
N = len(df['Phase A'])

frequency = rfftfreq(N, d=1.0/sr)
plt.plot(frequency, np.abs(rfft(df['Phase A'])))
plt.xlim(0.0, 480)
plt.show()


I have the following questions: What is the unit of the amplitude in the FFT spectrum if the input signal is in Amperes? I know the THD is $$100\frac{\sqrt{\sum_{k=2}^{\infty}V_k^2}}{V_1}$$ where $$V_k$$ are the harmonics and $$V_1$$ is the fundamental, but this formula doesn't seem to be applicable here. How do I calculate the THD from the FFT spectrum?

The output of a single-sided DFT such as this is scaled by the number of samples. Therefore, for an input time history that is $$N$$ samples in length, the DFT at each frequency point returns the product $$N\cdot A$$ where $$A$$ is the complex amplitude of the signal component at that frequency.