The scaling of FFT magnitude and phase

Question

In Matlab, I simulated and graphed out the FFT of a signal composed of three pure sinusoids of distinct magnitudes and phases. Figuring out the scaling of the frequency axis is fairly straightforward - the number of bins in the frequency domain is equivalent to the sample size in the time domain. So bin #31 corresponds to $ 31 \times F_s/N = 31 \times 1000/1500 \approx 20 \operatorname{Hz}$.

What I don't understand is the scaling of the y-axis (magnitude or the phase). Given the definition of DFT and the math derivation from FS to DFT, there has to be a scaling factor of $\frac{1}{NT_s}$ between the digital amplitude and the analog amplitude:

$$ C_k \,=\, \frac{1}{T_0} \int_{-\frac{T_0}{2}}^\frac{T_0}{2} x(t) \, e^{-j \omega_0 kt} \operatorname{dt} \,=\, \frac{1}{T_0} \int_{T_0} x(t) \, e^{-j \frac{2\pi}{T_0} kt} \operatorname{dt} \,=\, \frac{1}{NT_s} \sum_{n=0}^{N-1} x(nT_s) \, e^{-j \frac{2\pi}{NT_s} k nT_s} \,=\, \frac{1}{NT_s} \underbrace{\sum_{n=0}^{N-1} x(nT_s) \, e^{-j \frac{2\pi}{N} k n}}_{\text{DFT}} $$

Take bin #31 as an example, the analog amplitude that corresponds to $ 20\operatorname{Hz} $ is $3$ (or two spectral lines each with height $1.5$ at $ \pm 20 \operatorname{Hz} $ because $ \cos(\theta) = 0.5e^{j\theta}+0.5^{-j\theta} $). If I scale the unscaled magnitude of $x[31]$ I get $ 2250 \times \frac{1}{NT_s} = 2250 \times \frac{1}{1500\times 1000} = 0.0015 $, which is obviously wrong. If I want to get to the value of $ 1.5 $, I have to divide the digital magnitude by $ 1500 $, which is exactly my sampling rate .... why?

1. Why do I rescale the magnitude by $ 1/N $ whereas the math derivation suggests a scaling factor of $ 1/(NT_s) $?

2. If I have to scale the magnitude by $ 1/N $ I would have thought the phase will need the same scaling - after all, both the magnitude and phase originate from the same complex spectrum equation. But to my surprise, the FFT gives me the exact phase value without any rescaling - at bin #31, the phase is $0.2$! How is that possible?

Matlab code:

fs=1000;                    % sampling @1000Hz
t=0:1/fs:1.5-1/fs;          % from 0-1.499s in 0.001s interval so we have 1500 samples

f1=20;
f2=30;
f3=40;

x1=3*cos(2*pi*f1*t+0.2);    % @f1=20Hz, amplitude=3, phase=0.2rad
x2=1*cos(2*pi*f2*t-0.3);    % @f2=30Hz, amplitude=1, phase=-0.3rad
x3=2*cos(2*pi*f3*t+2.4);    % @f3=40Hz, amplitude=2, phase=2.4rad

x=x1+x2+x3;                 % x has 1500 samples; x is a real number
X=fft(x);                   % X has 1500 "bins" of frequencies 
X_mag = abs(X);             % get magnitude spectrum from the complex spectrum; X_mag is a real number

figure(1)                   % unscaled magnitude spectrum
plot(X_mag);

figure(2)                   % scaled magnitude spectrum
plot(X_mag/1500);           

figure(3)                   % unscaled phase spectrum
plot(angle(X));

FFT spectrum graphs:

Answer 1

Answering part 2:

The phase is the arctangent of the two components of the complex result of a FFT, or atan2(imaginary_part,real_part). If both components are scaled equally by any value (other than zero), the atan2() will remain unchanged.