function [ elementElementCoupling, elementCoupling, totalCoupling ] = couple( frequency, transmitGaindB, transmitBeamwidthDegrees, receiveGaindB, receiveBeamwidthDegrees, numTransmitColumns, numTransmitRows, numReceiveColumns, numReceiveRows, arraySpacingInches, receiveXIndex, receiveYIndex, transmitXIndex, transmitYIndex )
%COUPLE Summary of this function goes here
%   Detailed explanation goes here

impedance = 50;

lambda = 3e8/(frequency * 1e9);	%calculation and storage of wavelength
xElementSpacing = lambda * .6;  %Spacing between elements in the X direction
yElementSpacing = lambda * .6;  %Spacing between elements in the Y direction

transmitGain = 10^(transmitGaindB/10); %calculation of transmit gain
transmitBeamwidthRadians = transmitBeamwidthDegrees* pi/180; %calculation of transmit beamwidth in radians
receiveGain = 10^(receiveGaindB/10); %calculation of receive gain
receiveBeamwidthRadians = receiveBeamwidthDegrees* pi/180; %calculation of receive beamwidth in radians

arraySpacing = arraySpacingInches *  0.0254;    %Space between arrays in meters

power = 200; %arbitrary power value
radiatedPower = power * numTransmitColumns * numTransmitRows * transmitGain * receiveGain;   %radiated power calculation

%Variable Initialization
    totalVoltage = complex(0,0);	%total voltage (phasor) per receive element
    totalPower = 0;
    maxElementPower = 0;            %The highest power coupling to a single receive element
    fullArrayPower = 0;             %total coupling between the two arrays
    powerOnElements=zeros(numReceiveColumns, numReceiveRows);
    couplingBetweenElements=zeros(numReceiveColumns,numReceiveRows,numTransmitColumns,numTransmitRows);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%          Calculation of Power transferred to       %%%
%%%            Element of the Receive Array            %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=1:numReceiveColumns
    xReceive=xElementSpacing*(i-1); %x position of receive element
    for j=1:numReceiveRows
        yReceive=yElementSpacing*(j-1); %y position of receive element
        for k=1:numTransmitColumns
            xTransmit=-(arraySpacing+xElementSpacing*(k-1)); %x position of transmit element
            for l=1:numTransmitRows
                yTransmit=yElementSpacing*(l-1); %y position of transmit element
                
                distance=sqrt((xReceive-xTransmit)*(xReceive-xTransmit)+(yReceive-yTransmit)*(yReceive-yTransmit));    %distance between elements
				angle = atan((yTransmit-yReceive)/(xTransmit-xReceive));    %angle between transmit and receive elements
                angleDegrees = angle * 180 / pi;
                
                
                
                angleApproximate = 5 * round(angleDegrees / 5);
%                gainTransmit = 12*(angle/transmitBeamwidthRadians)*(angle/transmitBeamwidthRadians);
%                gainReceive = 12*(angle/receiveBeamwidthRadians)*(angle/receiveBeamwidthRadians);
                gainTransmit = transmitGain;
                gainReceive = receiveGain;

                elementPower=power * gainTransmit * gainReceive * lambda * lambda / (4 * pi * distance) / (4 * pi * distance); %calculation of power between these two elements (use of Friis Transmission Equation)
                
                couplingBetweenElements(i,j,k,l) = elementPower;
                
                voltageMagnitude=sqrt(elementPower/impedance);      %voltage magnitude
                phase=2 * pi * distance / lambda;           %voltage phase
                VR = voltageMagnitude * cos(phase);         %real part of voltage phasor
                VI = voltageMagnitude * sin(phase);         %imaginary part of voltage phasor
                V = complex(VR, VI);                        %voltage phasor
				totalVoltage = totalVoltage + V;            %add total voltage transferred to receive element (i,j) from transmit array
                totalPower = totalPower + elementPower;
            end
        end
        powerOnElements(i,j) = (real(totalVoltage)) * (real(totalVoltage)) / impedance;     %calculate and store the total power transferred from the transmit array for receive element (i,j)
        totalVoltage = 0 + 0i;                                   %reinitialize vTotal variable for use in next iteration
        %powerOnElements(i,j) = totalPower;
        totalPower = 0;
        
        if(powerOnElements(i,j) > maxElementPower)
            maxElementPower = powerOnElements(i,j);                       %if this is the element recieving the most power, set MAX as this value
        end
        fullArrayPower = fullArrayPower + powerOnElements(i,j);              %Add this received power to the total coupling value
    end
end

elementElementCoupling = 10 * log10(couplingBetweenElements(receiveXIndex,receiveYIndex,transmitXIndex,transmitYIndex)/power);
elementCoupling = 10 * log10(maxElementPower/radiatedPower);                    
totalCoupling = 10 * log10(fullArrayPower/radiatedPower);                    %calculate the total coupling between the full arrays

end