Моделирование радарных снимков на основе модельного поля волнений

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
 
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
 
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
 
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

    % 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);