% FUNCTION:  [Ez,r] = MakeWaves( L, f, thetaR, N )
%
% Makes sample fading plots as a function of space.
%
% Author:  Greg Durgin, 3/13/2004
%
% Inputs:
%  L - the total distance (in meters) to plot
%  f - carrier frequency (in Hz)
%  thetaR - azimuth orientation of movement in space (in degrees)
%  N - total number of plane waves (default is 100)
%
% Outputs:
%  Ez - vector of total z-component of E-field as function of space
%  r - vector of corresponding positions (same size as Ez)
%
% Sample Usage:
% To make a fading profile over 3 meters in space at a carrier frequency of
% 2.0 GHz as observed with a 45-degree orientation in space:
% >> [Ez,r] = MakeWaves( 3, 2.4e9, 45 );
%
% Notes:
% The code provides a useful example of how to simulate the small-scale
% fading that occurs when many plane waves superimpose in space.  Feel free
% to adapt this code for your ECE3065 project.  The output of this function
% is the z-component of electric field (in complex phasor form) as a
% function of space.  The function automatically graphs Ez as a function of
% r on a dB-scale graph at the end of the routine, so that you can see the
% types of outputs required for the project.
%
% The function "p" at the end of this file is your azimuth spectrum of
% multipath power.  The default spectrum is omnidirectional (so the graphs
% will not depend on thetaR at all!) but this can be changed fairly
% quickly.

function [Ez,r] = MakeWaves( L, f, thetaR, N )

% Set a default value of 100 for N
if ~exist('N')  
    N = 100;
end;

% Initialize Variables
lambda = 3e8/f;                 % wavelength
M = 16*L/lambda;                % good number of space samples
k = 2*pi/lambda;                % freespace wavenumber
r = linspace(0,L*(1-1/M),M);    % make an array of sequential positions
Ez = zeros(size(r));            % initialize Ez component to zero
theta = linspace(0,2*pi*(1-1/N),N); % vector azimuth angles


% Step through N azimuth angles and add z-components of electric
% field from each direction.  The expression below that is added to Ez each
% loop is a single plane wave.  Its amplitude is related to the multipath
% angle spectrum and its phase is random, but tapers "k" radians/meter in
% the direction of wave travel.  By the end of this loop, Ez will be the
% sum of N different plane waves, creating the strange
% constructive-destructive interference pattern as a function of space.
%
for n=theta,
    Ez = Ez + (p(n))^.5 * exp( j*( 2*pi*rand(1) + k*r*cos(n-thetaR*pi/180) ) );
end;

% Normalize Ez so that average power is 1
Erms = (sum(abs(Ez).^2)/length(Ez))^.5;     % root-mean-square
Ez = Ez/Erms;                               % normalize to average power 1

% Plot a nice graph of power in dB as a function of space
figure(1);
plot(r,20*log10(abs(Ez)),'b-');     % plot dB-scale Ez as a function of r
title(sprintf('Small-Scale Fading in Received Power at %1.0f MHz, \\theta_R=%i',f/1e6,thetaR*180/pi));
xlabel( 'Position (m)' );
ylabel( 'Normalized Received Power (dB)' );

return; % finish


% Here is the function defining the azimuth spectrum, the distribution of
% multipath power as a function of angle-of-arrival.
%
% The example below is for an omnidirectional spectrum (a constant).
% Change it to see other types of angle spectra.
%
function P = p(theta)
    P = 1/2/pi;
% Here is another example, the cosine-squared spectrum
% P = cos(theta)^2;
return;