Found the error in my derivation. My first derivation is not correct. The heaviside cannot be delayed by simply adding a delay term.
SInce $y(t)$ only has a value for $t>t_d$ the fourier transform yields $\int_{t_d}^{\infty}Ae^{-i\omega t}dt$.
The integral $\int_{t_d}^{\infty}Ae^{-i\omega t}dt = A\int_{t_d}^{\infty}cos(\omega t)dt + Ai\int_{t_d}^{\infty}sin(\omega t)dt $ is not defined. This also required for the derivation for the heaviside fourier transform where $\int_{0}^{\infty}Ae^{-i\omega t}dt$ is not defined
Using a briefer notation: Since $\int_{t_d}^{\infty}$ does not hold because of the infinity in the integral, it is possible to rewrite this to:
$\int_{t_d}^{\infty}Ae^{-i\omega t}dt = \int_{0}^{\infty}Ae^{-i\omega t}dt - \int_{0}^{t_d}Ae^{-i\omega t}dt$
The second term ($\int_{0}^{t_d}$) is finite, and therefore defined. The first term is can derived analogous to the fourier transform of the 'normal' (without delay) heaviside function. The first integral on the rhs is the correction term, if you will, which is a pulse in time domain and a sinc in the frequency domain.
The discrete Fourier transform implicitly regards the input data as periodic. The length of the input data is an integer multiple of the frequencies that the output data correspond to. If you generate the input data from a sinusoidal that has a different period that doesn't evenly divide the length of the data, then you're effectively creating a jump at the point where the first data point follows the last in the periodic continuation. A jump has appreciable Fourier components at all frequencies, and this is superimposed on your signal.
Nevertheless, you do get a signal with a clear maximum around the input frequency. I tried this using this online FFT calculator. In your case, as you calculated, we'd expect the maximum to be around the $21$st and $22$nd output data point (with the first data point corresponding to frequency $0$), and indeed, the first $40$ output values are:
30.322000,0.000000
30.391603,-0.332305
30.619103,-0.635705
30.983132,-0.983397
31.541176,-1.314554
32.258933,-1.685996
33.213385,-2.076468
34.383753,-2.509669
35.850294,-2.989665
37.679693,-3.538320
39.950279,-4.172808
42.781442,-4.897551
46.400390,-5.801462
51.088228,-6.910079
57.343876,-8.338977
66.017176,-10.282208
78.721064,-13.088380
98.974578,-17.487260
136.068697,-25.448952
225.318492,-44.421047
727.024048,-150.774758
-543.414483,118.240910
-191.999756,43.727082
-114.594042,27.275221
-80.661116,20.008091
-61.661208,15.918785
-49.536603,13.289953
-41.141647,11.459233
-35.002598,10.087933
-30.321248,9.041060
-26.630074,8.216618
-23.664950,7.534937
-21.230810,6.967688
-19.199079,6.488785
-17.492157,6.083581
-16.004951,5.723743
-14.733225,5.416510
-13.616822,5.138707
-12.640339,4.905611
-11.768659,4.665358
So it seems you must have been doing something wrong when you got spikes all over the place at roughly the same amplitude.
You can actually calculate the output in closed form in this case by writing the sinusoidals in terms of complex exponentials and summing the resulting partial geometric series.
Edit:
Here's that last suggestion written out.
You have the time-domain signal $\sin2\pi\nu t$ with $\nu=440\text{Hz}$, and you sample it at the frequency $f=44100\text{Hz}$, so you have the samples
$$
s_k=\sin\frac{2\pi\nu k}f=\frac1{2\mathrm i}\left(\exp\left(\frac{2\pi\mathrm i\nu k}f\right)-\exp\left(-\frac{2\pi\mathrm i\nu k}f\right)\right)
$$
for $0\le k\lt n$, with $n=2048$. Performing a discrete Fourier transform on these samples yields
\begin{eqnarray*}
\tilde s_j
&=&
\sum_{k=0}^{n-1}s_k\exp\left(\frac{2\pi\mathrm i jk}n\right)
\\
&=&
\sum_{k=0}^{n-1}\frac1{2\mathrm i}\left(\exp\left(\frac{2\pi\mathrm i\nu k}f\right)-\exp\left(-\frac{2\pi\mathrm i\nu k}f\right)\right)\exp\left(\frac{2\pi\mathrm i jk}n\right)
\\
&=&
\frac1{2\mathrm i}\sum_{k=0}^{n-1}\left(\exp\left(2\pi\mathrm i\left(\frac jn+\frac\nu f\right)k\right)-\exp\left(2\pi\mathrm i\left(\frac jn-\frac\nu f\right)k\right)\right)
\\
&=&
\frac1{2\mathrm i}\sum_{k=0}^{n-1}\left(\omega_{j+}^k-\omega_{j-}^k\right)
\\
&=&
\frac1{2\mathrm i}\left(\frac{1-\omega_{j+}^n}{1-\omega_{j+}}-\frac{1-\omega_{j-}^n}{1-\omega_{j-}}\right)\end{eqnarray*}
with $\omega_{j\pm}=2\pi\mathrm i\left(\frac jn\pm\frac\nu f\right)$. The second term dominates and is maximal at $\frac jn=\frac\nu f$ (which in your case is at the fractional value $j=\frac{n\nu}f$). It's a phase factor times
$$
\frac{\sin n\pi\left(\frac jn-\frac\nu f\right)}{\sin\pi\left(\frac jn-\frac\nu f\right)}\;,
$$
which is a discrete periodic analogue of the $\operatorname{sinc}$ function. Fitting to this functional form might improve your estimate of the frequency from the output data of the Fourier transform.
Best Answer
The Fourier transform of $\sin(\omega_0 t)$ is $$ \mathcal{F}\left[\sin(\omega_0 t)\right] = -\frac{i}{2}\delta(\omega -\omega_0) + \frac{i}{2}\delta(\omega + \omega_0) $$ where $\delta$ is the Dirac delta function. This is seen in your first image. The magnitude of this function is $$ |\mathcal{F}\left[\sin(\omega_0 t)\right]| = \frac{1}{2}\delta(\omega -\omega_0) + \frac{1}{2}\delta(\omega + \omega_0) $$ as seen in your second image. So then what about the phase? Well, the delta function is zero everywhere except where its argument is zero. That means that for all $\omega \ne \pm \omega_0$, $\mathcal{F}\left[\sin(\omega_0 t)\right](\omega) = 0$. What's the phase of zero? It's undefined. So your phase plot won't be much of a plot.
The reason the plot in the linked question is so noisy is that numerical error is causing the values of $\mathcal{F}\left[\sin(\omega_0 t)\right]$ to not compute to exactly zero. Instead it computes to an extremely small magnitude complex number whose phase is essentially random.