% In the name of Allah, the Most Merciful and the Beneficent

% Ripple Localization Algorithm

%{ 
M. Farooq-I-Azam; Q. Ni; E. Ansari, “Intelligent Energy Efficient 
Localization Using Variable Range Beacons in Industrial Wireless Sensor 
Networks,” in IEEE Transactions on Industrial Informatics , vol.PP, no.99, 
pp.1-1 
doi: 10.1109/TII.2016.2606084
%}

% This is the simulation engine for the Ripple Localization Algorithm
% (RLA). This engine can be run with loops and other wrappers and results
% for different scenarios can be obtained. For example, we can use an outer
% loop which changes the number of beacon nodes inside the loop so that
% results for different number of beacon nodes can be obtained.
% 
% Copyright: Muhammad Farooq-i-Azam
% Email: ripplelocalizationalgorithm@gmail.com

% Infolab21, School of Computing and Communications
% Lancaster University, Lancaster, UK.

% Department of Electrical Engineering
% COMSATS Institute of Information Technology, Lahore, Pakistan.

% License: GPL


FileName = 'BasicEngine001.mat';

DOI = 0.2;

TOLERANCE = 0.01;

NumDumbNodes = 100;
DumbNodeRange = 10;
NumBeaconNodes = 20; 
BeaconMinRange = 10;
BeaconSignalStep = 10;
BeaconMaxRange = 100;

DIMENSIONS = 2;
FieldLength = 100;
FieldWidth = 100;
FieldHeight = 0;

SensorFieldArea = FieldLength * FieldWidth;

% Array containing actual positions of dumb nodes.
% This array will be subsequently expanded and contain radii
% of inner and outer circles of beacon node for each dumb node.
% This array will also hold estimated positions of dumb nodes and other
% results.
DumbNodes = FieldLength .* rand(NumDumbNodes, DIMENSIONS);

% Array containing actual positions of beacon nodes
BeaconNodes = FieldLength .* rand(NumBeaconNodes, DIMENSIONS);

% Calculate actual distances of dumb nodes from each beacon node.
ActualDistance = pdist2(DumbNodes, BeaconNodes);

% Create temporary array of -1s to hold radii of inner and outer circles
% of the ring around a beacon node in which the dumb node lies. This matrix
% will be used to expand the DumbNodes Matrix.

TwoColumns = ones(NumDumbNodes, DIMENSIONS);
TwoColumns = - TwoColumns;

% Determine the radii of the rings in which the dumd nodes lie.
% For each beacon node, two columns are added to DumbNodes matrix to 
% hold radius of inner and outer circles

for CurrentBeaconNode = 1:NumBeaconNodes
   DumbNodes = [DumbNodes TwoColumns]; 

   for BeaconCurrentRange = ...
                BeaconMinRange : BeaconSignalStep : BeaconMaxRange 
      for CurrentDumbNode = 1 : NumDumbNodes
         Distance = ActualDistance(CurrentDumbNode, CurrentBeaconNode);

	 RBeaconCurrentRange = ...
                normrnd(BeaconCurrentRange, DOI*BeaconCurrentRange);

         if BeaconCurrentRange == BeaconMinRange

            if Distance < RBeaconCurrentRange
              
               DumbNodes(CurrentDumbNode, ...
                   DIMENSIONS*CurrentBeaconNode + 1) = BeaconMinRange;

            end
     
	 elseif Distance > (BeaconCurrentRange - BeaconSignalStep) ...
                    && Distance < RBeaconCurrentRange

	    DumbNodes(CurrentDumbNode, DIMENSIONS*CurrentBeaconNode + 1) ...
                    = BeaconCurrentRange - BeaconSignalStep;

	    DumbNodes(CurrentDumbNode, DIMENSIONS*CurrentBeaconNode + 2) ...
                    = BeaconCurrentRange;

	 end 

      end 

   end 	       

end 

% Matrix DumbNodes now has dimensions of 
% NumDumbNodes x DIMENSIONS + DIMENSIONS * NumBeaconNodes.
% Matrix DumbNodes now has the radii of the inner and outer circle rings 
% from each beacon node for every dumb node in the row in which actual 
% position coordinates of the dumb node are present.


BASE = DIMENSIONS + DIMENSIONS * NumBeaconNodes;

NodeDensity = NumDumbNodes / SensorFieldArea

NodesPerNodeArea = NodeDensity * pi*DumbNodeRange^2

NodesAreaToFieldArea = NumDumbNodes*pi*DumbNodeRange^2 / SensorFieldArea

UBR = NumDumbNodes / NumBeaconNodes

PercentBeaconNodes = NumBeaconNodes * 100 / NumDumbNodes

Frequency = 1+(BeaconMaxRange - BeaconMinRange) / BeaconSignalStep

ANR = BeaconMaxRange / DumbNodeRange

DumbNodes = [DumbNodes, TwoColumns];

DumbNodes = [DumbNodes, TwoColumns];

DumbNodes = [DumbNodes, TwoColumns];
DumbNodes = [DumbNodes, TwoColumns];

DumbNodes = [DumbNodes, TwoColumns];
DumbNodes = [DumbNodes, TwoColumns];


for CurrentDumbNode = 1 : NumDumbNodes
    Neighbors =[];      
    
    % Matrix r will hold estimated ranges from the beacon nodes
    r = [];
    
    for CurrentBeaconNode = 1 : NumBeaconNodes

        Column = DIMENSIONS * CurrentBeaconNode + 1;
        r1 = DumbNodes( CurrentDumbNode, Column ); % Inner Circle
        if ( r1 > 0 )
            r2 = DumbNodes( CurrentDumbNode, Column + 1 );
            Range = (r1 + r2) / 2;    % Estimated range
            r = [r ; Range];
            Neighbors = [Neighbors ; BeaconNodes(CurrentBeaconNode, :)];
        end    
    end    

    % After above for loop, Neighbors contains coordinates of all beacon 
    % nodes from which the dumb node receives signals. Range contains 
    % distances of the dumb node from these beacon nodes.
    
    % We can now use multilateration to estimate position
    
    % Number of neighbor beacon nodes
    
    k = size(Neighbors, 1);
    
    if ( k >= 3)
        
        xb = Neighbors(:, 1);
    
        yb = Neighbors(:, 2);
    
        A = [];
        b = [];
        for i = 1 : k - 1
            A = [A; Neighbors(k, :) - Neighbors(i, :)];
            b = [b; r(i)^2 - r(k)^2 + xb(k)^2 - xb(i)^2 + yb(k)^2 - ...
                                                            yb(i)^2];
        end 
    
        position = A \ b;
    
        position = position ./ 2;
        
        position = position';
    
    
        DumbNodes(CurrentDumbNode, BASE + 1) = position(1, 1);
        DumbNodes(CurrentDumbNode, BASE + 2) = position(1, 2);
        
    
        ActualX = DumbNodes(CurrentDumbNode, 1);
        ActualY = DumbNodes(CurrentDumbNode, 2);
        MLocalizationError = pdist([ActualX, ActualY; position], ...
                                                        'euclidean');
        DumbNodes(CurrentDumbNode, BASE + 3) = MLocalizationError;
        
        DumbNodes(CurrentDumbNode, BASE + 4) = 0; 
        

        
    end
    
end 

Results = [ DumbNodes(:, 1:2) DumbNodes(:, BASE+1 : end ) ];

NumBlindNodes2 = sum ( DumbNodes(:, BASE + 3) < 0 );

NumSettledNodes2 = sum ( DumbNodes(:, BASE + 3) >= 0 );

LocalizationEfficiency2 = NumSettledNodes2 * 100 / NumDumbNodes;
Errors2 = DumbNodes(:, BASE + 3); 

Negatives2 = Errors2 < 0;

Errors2 ( Negatives2 ) = 0;

MeanLocalizationError2 = mean ( Errors2 );

MedianLocalizationError2 = median ( Errors2 );

ModeLocalizationError2 = mode (Errors2 );

NormalError12 = Errors2 / BeaconSignalStep;

AverageNormalError12 = mean(NormalError12);

NormalError22 = Errors2 / BeaconMaxRange;

AverageNormalError22 = mean(NormalError22);

% Localization Error normalized to dumb node range. It is this 
% normal error which is usually used
NormalError32 = Errors2 / DumbNodeRange;

NormalError3NoBlinds2 = NormalError32 ( NormalError32 > 0 );

AverageNormalError32 = mean( NormalError3NoBlinds2 )

MedianNormalError32 = median( NormalError3NoBlinds2 )

ModeNormalError32 = mode ( NormalError3NoBlinds2 )

RangeNormalError32 = range ( NormalError3NoBlinds2 )



save(FileName);






