Observing Artefacts along azimuth direction after applying RMA to FMCW radar data

Question

I have simulated 3 scatterers for FMCW radar signal processing using Range Migration Algorithm.

Background:
The Range Migration algorithm steps are as follows:

1. Applying Azimuth FFT on the raw data
2. Apply "Matched Filter"
3. Apply "Stolt Interpolation"
4. Apply 2D Inverse Fourier Transform on the resulting signal

Signal Model:

Starting Signal dechirped signal model with RVP term removed,

After Azimuth Fourier Transform,

What I am observing

What I should be observing:

Y-axis : Azimuth direction X-axis : Range direction

Question:

So it seems the scatterer in the final image is diffused along the azimuth direction? But I am not sure what might be causing this diffusion along the azimuth axis. Could someone suggest what I should be looking out for? In terms of system parameters ? Or it might be due a step I am performing incorrectly? A few potential directions to look at is very much appreciated.

The following code is meant for someone who wants to advice on how to tackle the signal processing problem without actually having to implement the code from scratch

Data generation code for reference:

c = 3e8 # speed of light constant 
B = 3.6e9 # Bandwidth
pr = c/2/B # range resolution

PRF = 2500/3 # Pulse repetition frequency
PRI = 1/PRF # Pulse repetition Interval
Tp = PRI # Pulse width
Cr = B/Tp # Chirp rate

Fc = 3.48*1e10 # Center frequency
Fs = 25e6 # system sampling rate 
ts = 1/Fs # sampling interval

Nrg = round(Fs * Tp) # Number of range samples
t = np.linspace(-Nrg/2,Nrg/2,Nrg)*ts # Fast time 

radarVel = np.array([4,0,0]) # speed of aircraft in m/s
pa =  0.054 # azimuth resolution

# Derieved Parameters
Rs = 118 # slant distance to scene center
lamC = c/Fc # wavelength
delTh = (lamC)/(2*pa) # integration angle
L = Rs*delTh/sin(np.deg2rad(90)) # Synthetic Aperture Length (2.21)
Ta = L/eucDist(radarVel) # Synthetic Aperture time (2.22)
Naz = round(PRF * Ta /2) * 2  #Number of Azimuth samples (matter of convenience)

tSeq = np.arange(0,Naz) * PRI # timings

H = 34
radarPosc = np.array([0,np.sqrt((Rs**2)-(H**2)),H])# Radar Coords at Aperture center along the path
radarPoss = radarPosc - (radarVel * Naz//2 * PRI) # Radar's Starting Position


radarSeq = np.zeros((3,len(tSeq)))
radarSeq[0,:] = radarPoss[0] + tSeq * radarVel[0]
radarSeq[1,:] = radarPoss[1] + tSeq * radarVel[1]
radarSeq[2,:] = radarPoss[2] + tSeq * radarVel[2]

# Data generation
raw = np.zeros((Naz,Nrg), dtype=np.complex_)

KR = (4*pi/c)*(Fc+Cr*t) # Kr Spatial frequency (3.11)

for n in range(Naz):
    for scat in scattererCoords:

        R  = np.sqrt((radarSeq[1,n]-scat[1])**2 + (radarSeq[2,n]-scat[2])**2)
        x_comp = radarSeq[0,n]-scat[0] + norm(radarVel)*t # includes the motion between the pulses
        phase = np.exp(1j*KR*np.sqrt(R**2 + x_comp**2)) # without RVP term
        tau = 2*np.sqrt(R**2 + x_comp**2)/c 
        raw[n,:] += rect(t,Tp,tau) * phase

Range Migration Algorithm Procedure

# # 1D FT along azimuth track
S = fftshift(fft(raw, axis=0),axes=0)

# # Matched Filtering
phi_mf = -Rs*KYY   # Matched Filter
S_mf = S * np.exp(1j*phi_mf)

# Stolt interpolation
kstart = KYY[Naz//2,:].min()
kend   = KYY[Naz//2,:].max()

KYq = np.linspace(kstart, kend, int(Nrg)) # query points, create a uniformly sampled points

S_st = np.zeros((len(KYY), len(KYq)), dtype = 'complex_')

for i in range(Naz):
    S_st[i,:] = np.interp(KYq, KYY[i,:], S_mf[i,:],
                          left = 0, right = 0)
    
S_final =ifft2((S_st))

Edit: Have included the 1-d Azimuth plots which are in the range-time and azimuth frequency domain.

Close-up of the final processed image: