% Besmellah...
clearvars;
close all

% parameter definition:
c = 3e8;
lambda = 0.025;
La = 0.5;
B = 200e6;
Tp = 5e-6;

tet_bw = lambda/La;
f0 = c/lambda;
slant_renge_res = c/2/B;
TBP = B*Tp;
minR = c*Tp/2;
Kr = B/Tp;
osf = 1.2;
Fs = B*osf;
ts = 1/Fs;


% baseband signal:
tb = -Tp/2:ts:Tp/2;
Ntb = length(tb);
sb = exp(1i*pi*Kr*tb.^2);

figure;plot(tb,real(sb)); 
title('baseband signal');
grid on;
xlabel('t(s)')

freqb = (-Ntb/2:Ntb/2-1)/Ntb*Fs;
figure; plot(freqb,db(nrml(fftshift(fft(sb)))));
title('FFT of baseband signal');
xlabel('f(Hz)');
grid on;

smf = conv(sb,conj(fliplr(sb)),'same');
figure;plot(tb,20*log10(nrml(abs(smf))));
title('auto-correlation of baseband signal');
grid on;
xlabel('t(s)')


% SAR and scene Geometry:
h = 8000;
xt = 8000;
yt = 0;
zt = 0;
R0 = sqrt(h^2+xt^2);
Ls = tet_bw*R0;
Yp1 = -Ls/2*1.2;
Yp2 = -Yp1;
az_os = 1.2;
dy = La/2/az_os;
yp = Yp1:dy:Yp2;
Np = length(yp);
xp = zeros(1,Np);
zp = ones(1,Np)*h;
Rs = sqrt((xp-xt).^2+(yp-yt).^2+(zp-zt).^2);
Rmin = min(Rs);
Rmax = max(Rs);

% RAW data simulation:
t = 2*Rmin/c-Tp/2:ts:2*Rmax/c+Tp/2;
Nt = length(t);
RawData = zeros(Np,Nt);
Range_Compressed_Data = zeros(Np,Nt);
for np=1:Np
    R = Rs(np);
    delay = 2*R/c;
    teta = asin((yp(np)-yt)/R);
    Inbeam = abs(teta)<=tet_bw/2;
    RawData(np,:) = exp(1i*pi*Kr*(t-delay).^2)*exp(-1i*4*pi/lambda*R)*...
        Inbeam.*(t>=delay-Tp/2).*(t<=delay+Tp/2);
    Range_Compressed_Data(np,:)=conv(RawData(np,:),conj(fliplr(sb)),'same');
end
figure;imagesc(real(RawData));
title('RAW Data')
xlabel('Range')
ylabel('Azimuth')
figure;imagesc(abs(Range_Compressed_Data));
title('Range-Compressed Data');
xlabel('Range')
ylabel('Azimuth')


% RCMC:
Range_Doppler_Data = fftshift(fft(Range_Compressed_Data,[],1),1);
figure;imagesc(abs(Range_Doppler_Data));
title('Range-Compressed Data in Range-Doppler Domain');
xlabel('Range')
ylabel('Doppler')

PRF = 1000;
Vr = dy*PRF;
feta = (-Np/2:Np/2-1)'/Np*PRF;
RCMC_Data = zeros(Np,Nt);
dR = c*ts/2;

for ii=1:Np
    dRrd = lambda^2*R0*feta(ii).^2/8/Vr^2;
    smpl_shft = round(dRrd/dR);
    fft_rng = fft(Range_Doppler_Data(ii,:));
    shft_fft_rng = fft_rng.*exp(1i*2*pi/Nt*(0:Nt-1)*smpl_shft);
    RCMC_Data(ii,:) = ifft(shft_fft_rng);
end
figure;imagesc(abs(RCMC_Data));
title('After RCMC');
xlabel('Range')
ylabel('Doppler')

% Azimuth Compression:
R0near = sqrt(h^2+(xt-100)^2);
R0far = sqrt(h^2+(xt+100)^2);
R0s = R0near:dR:R0far;
indR0s = floor(Tp/2/ts)+floor((R0s-Rmin)/dR)+1;
Ka = 2*Vr^2/lambda./R0s;
[KA,Feta] = meshgrid(Ka,feta);
Img = zeros(Np,Nt);
Img(:,indR0s) = RCMC_Data(:,indR0s).*exp(-1i*pi*Feta.^2./KA);
Img(:,indR0s) = fft(Img(:,indR0s));
figure;imagesc((nrml(abs(Img))));
title('Image');
xlabel('Range')
ylabel('Azimuth')