Simulation of radar images based on a simulated wave field

Question

I am modeling a radar image based on my already modeled field of sea surface waves. Faced with the problem that the final image does not work right. It turns out to be some kind of nonsense. You can see the resulting image and what the field of surface waving looks like.

The script in MATLAB looks like this (without taking into account main_prog from which all these functions are called).

Functions of parameters:

function [g, speed, dx, dy, x, y, df, f_start, f_end, dtt, theta] = common_data()   
    g = 9.81; % acceleration of gravity
    speed = 13; % wind speed

    dx = 1; % step along the x-axis
    dy = 1; % step along the y-axis
    x = 0:dx:1000; 
    y = 0:dy:1000; 

    df = 0.1; % frequency range step
    f_start = 0.016; % frequency range start
    f_end = 1.25; % frequency range end

    dtt = pi / 18; % angle step 
    theta = -pi/2:dtt:pi/2; % angle range
end

function [xd, alpha, wmaxd, wmax, gm] = JONSWAP_param(g, speed)
    xd = 2 * 10^(4); % acceleration parameter
    alpha = 0.076 * xd^(-0.22); % Phillips constant
    wmaxd = 3.5 * xd^(-0.33); % dimensionless frequency
    wmax = wmaxd * (g / speed); % max frequency of spectrum
    gm = 3.3; % gamma-factor
end

Spectrum calculations:

function [JSP] = JONSWAP_Spec(g, alpha, wmax, gm, f)
    sgm = zeros(size(f)); 
    EXP = zeros(size(f));  
    JSP = zeros(size(f)); 
    
    for i = 1:length(f)
        if wmax >= f(i) 
            sgm(i) = 0.07; 
        else
            sgm(i) = 0.09; 
        end
        
    EXP(i) = exp(- ( (f(i) - wmax)^2 ) / ( 2 * sgm(i)^2 * wmax^2) );  
    JSP(i) = alpha * g^(2) * (2*pi)^(-4) * ( (f(i))^(-5) ) * ( exp( (-5/4) * (f(i)/wmax)^(-4) ) ) * ( gm^(EXP(i)) ); % JONSWAP spectrum
    end 
end

function [St, Swt, eta0] = Freq_Ang_Spec(dtt, theta, df, JSP)
    St = cos(theta).^(4); % angular spectrum
    Norm = trapz(dtt, St); 

    Swt = JSP .* St'; % freq-angular spectrum
    
    eta0 = sqrt((2 * Swt .* df * dtt) ./ Norm); % amplitude of surface waving 
end

Surface Waving simulation:

function [phase, Kx, Ky, etaWA] = SurfaceWaving_simulator(g, x, y, theta, w, eta0)
    phase = 2*pi*rand(length(theta),length(w)); % initial phase

    t = 0; % time moment

    Kabs = (w.^2) / g; % wavenumber 
    Kx = Kabs .* cos(theta)'; % projection of the wavevector on the x-axis 
    Ky = Kabs .* sin(theta)'; % projection of the wavevector on the y-axis 
    
    etaWA = zeros(length(x),length(y)); 
    
    for i = 1:length(x)
        for j = 1:length(y)
            etaWA(i,j) = sum(eta0 .* cos(w * t - Kx * i - Ky * j + phase), 'all');
        end
    end
end

RADAR data and RADAR image simulation:

function [az_ang, inc_ang, slope_x, slope_y, Rfl_VV_s, Rfl_HH_s, K_rad, K_b, f_b] = Geom_Sig_Settings(g, x, y, etaWA)
    H_rad = 14; % radar height above sea level

    r = sqrt(x.^2 + y.^2); % horizontal distance 
    R = sqrt(r.^2 + H_rad^2); % the trajectory of the signal to each point of the surface

    az_ang = atan2(y, x); % azimuth angle of the radar
    inc_ang = acos(H_rad./R); % angle of incidence   
    graz_ang = pi/2 - inc_ang; % grazing angle  

    % Calculation of local surface slopes (geometric modulation of the returned signal)
    [sx, sy] = gradient(etaWA);
    slope_x = atan(sx); % slope on x
    slope_y = atan(sy); % slope on y
    % Converting slopes to the radar coordinate system
    slope_rad_x = slope_x .* cos(az_ang) + slope_y .* sin(az_ang); 
    slope_rad_y = -slope_x .* sin(az_ang) + slope_y .* cos(az_ang);

    % Calculation of geometric reflection coefficients 
    Rfl_VV = @(inc_ang) (cos(inc_ang)).^(4) .* (1 + (sin(inc_ang)).^(2)).^(2) / (cos(inc_ang) + 0.111).^(4); % the square of the reflection coefficient module (case VV)
    Rfl_HH = @(inc_ang) (cos(inc_ang)).^(4)  / (0.111 .* cos(inc_ang) + 1).^(4); % the square of the reflection coefficient module (case HH)

    Rfl_VV_s = sqrt(Rfl_VV(inc_ang + slope_rad_x)); % reflection coefficient module (case VV)
    Rfl_HH_s = sqrt(Rfl_HH(inc_ang + slope_rad_x)) + (slope_rad_y / sin(inc_ang)).^(2) .* sqrt(Rfl_VV(inc_ang)); % reflection coefficient module (case HH)

    %% Radar signal parameters
    c_light = 3 * 10^(8); % light speed
    f_rad = 10 * 10^(9); % the frequency of the radar signal emitted by the radar
    lambda_rad = c_light / f_rad; % the wavelength of the radar signal
    K_rad = (2*pi) / lambda_rad; % the wavenumber of the radar signal

    K_b = 2 * K_rad * sin(inc_ang); % the Bragg-Wolfe condition
    sigma = 0.073 / 997; 
    f_b = sqrt(g * K_b + sigma * (K_b).^(3)) / (2*pi); % frequency of Bragg waves
end

function [JSP_b, RCS_VV] = RadBackscat_Simulator(g, alpha, wmax, gm, az_ang, Rfl_VV_s, Rfl_HH_s, K_rad, K_b, f_b)
    % Calculation of JONSWAP spectrum using f_B
    JSP_b = JONSWAP_Spec(g, alpha, wmax, gm, f_b);

    RCS_VV = 16 * pi * (K_rad)^4 .* (abs(Rfl_VV_s))^2 .* JSP_b; 

    %RCS_HH = 16 * pi * (K_rad)^4 .* (abs(Rfl_HH_s))^2 .* JSP_b;

    %% Radar image visualization
    figure(3)
    imagesc(log10(RCS_VV));
    colormap gray;
    colorbar;
    title('RADAR image');
    xlabel('x');
    ylabel('y');
end

At the moment, I'm not sure what the problem might be, but nothing else comes to mind. I think the whole problem lies in the fact that my model does not take into account the position of the radar, because depending on the azimuthal angle az_ang and angle of incidence inc_ang at which the radar "looks" at a particular point on the ocean surface, the slope values slope_x slope_y will depend. That is, in fact, if we change the position of the radar, it will also change slopes. In my case, the slopes are counted without taking this fact into account.

I tried to take this fact into account when calculating the slopes, but it didn't give any results:

    % Converting slopes to the radar coordinate system
    slope_rad_x = slope_x .* cos(az_ang) + slope_y .* sin(az_ang); 
    slope_rad_y = -slope_x .* sin(az_ang) + slope_y .* cos(az_ang);

What am I doing wrong? What could be the problem?

I took the formulas for reflection coefficients here ( page 3 formulas (3),(4) & (6),(7) ):

https://www.researchgate.net/publication/47652627_A_semiempirical_model_of_the_normalized_radar_cross-section_of_the_sea_surface_-_1_Background_model