function backscatter(par,type)
% The input parameter "type" is a character a,b,c,d,or e.
% The input parameter "par" depends on the type, look at
% the function makeSurf() to see what the "par" is.

N     = 2000;   % the number of grid points (1-D)
angle = 0.01;   % the antenna beam width from center
h     = 717e3;  % altitude of satellite
dr    = 0.5;    % range resolution.

% lay out the spacial grid on the sea surface.
x_max = 1.25*h*angle;                   % Size of grid
x = x_max * 2*((0:(N-1))/(N-1)-0.5);    % The grid array
y = x';
rho2 = x.^2 + y.^2;                     % Distance from the center squared.

alpha = 0.5/(h*angle)^2;                % A temporary constant
G = exp(-alpha*rho2);                   % The antenna gain

z = makeSurf(x,y,par,type);             % Compute the surface height.
r = sqrt((h-z).^2 + rho2);              % Range to satellite.
d = r-h;                                % Subtract altitude.

d_min = -20*dr;                         % Location of first bin.

Np = 128;
dgrid = d_min+(0:(Np-1))*dr;            % Create range bins.

P = NaN(1,Np);                          % Allocate the power array.

for n = 1:Np                            % Loop over range bins.
    d_target = dgrid(n);                % The range bin in this round.
                                        % The next line does the magic.
    inBin = abs(d-d_target)<(dr/2);     % Array of 1's and 0's.
    P(n)  = sum(sum( G.*inBin ));       % Power in bin n.
end

figure(2)
plot(dgrid,P,'.')                       % plot the power versus range
ylim([0,1.1*max(P)])                    % set the range of the plot

end

% This function creates the surface height array.
function z = makeSurf(x,y,par,type)
    H = par(1);
    Nx = size(x,2);
    Ny = size(y,1);
    switch type
        case 'a'
            z = 0;
        case 'b'
            T = par(2);
            theta = par(3);
            k = 4*pi^2/(9.8*T^2);
            kx = k*cosd(theta);
            ky = k*sind(theta);
            z = 0.5*H*cos(kx*x +ky*y);
        case 'c'
            z = H*(rand(Ny,Nx)-0.5);
        case 'd'
            z = 0.25*H*randn(Ny,Nx);
        case 'e'
            Hs = par(2);
            T = par(3);
            theta = par(4);
            k = 4*pi^2/(9.8*T^2);
            kx = k*cosd(theta);
            ky = k*sind(theta);
            z = 0.5*Hs*cos(kx*x +ky*y);
            figure(10)
            imagesc(z)
            z = z + 0.25*H*randn(Ny,Nx);
    end
            figure(11)
            imagesc(z)
end
        