X = load('Cell_Info.mat');
Names = fieldnames(X);
close all;
n=3;
c = 3E8;
%A = Cell_Info_A.ReceivedPowerMap;

%for cell= 1:length(Names),
for cell= 1:4,
    Cell_Info = getfield(X,Names{cell});
    freq = Cell_Info.freq;
    f = (freq)*10^6;
    APL=22;
    lambda = c/f;
    lambdaloss = 20*log10((lambda)/(4*pi));
    my_raster=Cell_Info;
    Compare = zeros(203,203);
    A = Cell_Info.ReceivedPowerMap;
    
for i=1:length(Cell_Info.ReceivedPowerMap),
    for j=1:length(Cell_Info.ReceivedPowerMap),
      

           if (i~=101)|(j~=101)
               my_raster.ReceivedPowerMap(j,i)=0.7*Cell_Info.erp - 10*n*log10(sqrt((Cell_Info.BS_x - i)^2 + (Cell_Info.BS_y - j)^2)) + lambdaloss-APL;
           end
            
    end
end
n=2;
for i=-length(Cell_Info.ReceivedPowerMap)/2:length(Cell_Info.ReceivedPowerMap)/2,
    for j=1:length(Cell_Info.ReceivedPowerMap)-101,
        if j>abs((sqrt(3)/3)*i)
            C=[cos((Cell_Info.azimuth-90)*2*pi/360),-sin((Cell_Info.azimuth-90)*2*pi/360);sin((Cell_Info.azimuth-90)*2*pi/360),cos((Cell_Info.azimuth-90)*2*pi/360)]*[i;j];
            
            C(1) = C(1)+101;
            C(2) = C(2)+101;
            x = int16(C(1));
          
            y = int16(C(2));
            if x<=0
                x = 1;
            end
            if y<=0
                y=1;
            end
   
            if (x~=101)|(y~=101)
                my_raster.ReceivedPowerMap(x,y)=Cell_Info.erp - 10*n*log10(sqrt((Cell_Info.BS_x - double(x))^2 + (Cell_Info.BS_y - double(y))^2)) + lambdaloss-APL;
            end
            
            if my_raster.ReceivedPowerMap(x,y)>0
                my_raster.ReceivedPowerMap(x,y)=0;
            end
        end
    end
end
my_raster.ReceivedPowerMap(101,101) = 0;
    my_raster.ReceivedPowerMap = my_raster.ReceivedPowerMap(1:203,1:203);
    for i=1:203,
        for j=1:203,
            if A(i,j)~=0
                
                Compare(i,j) = Compare(i,j) + my_raster.ReceivedPowerMap(i,j) - A(i,j);
             
            end
        end
    end
    
if cell == 1
    Cell_Info_A = Cell_Info;
    Compare_master = Compare;
end
if cell >1
    Image_dims = size(Compare);    
    newCompare = zeros(Image_dims);
    
    for im_x = 1:Image_dims(1),
        for im_y = 1:Image_dims(2),
            x1 = double(im_x) +(Cell_Info.lon- Cell_Info_A.lon)*9263.6;
            y1 = double(im_y) + (Cell_Info.lat- Cell_Info_A.lat)*11091.8;
            if (x1<Image_dims(1))&(y1<Image_dims(2))&(int16(x1)>0)&(int16(y1)>0)
                 newCompare(int16(x1), int16(y1)) = Compare(im_x,im_y);
            end
        end
    end 
  
    Compare_master = Compare_master+newCompare;
end

        Title = Names{cell}; Title( find(Title=='_') ) = ' ';
        figure(cell); imagesc(my_raster.ReceivedPowerMap); axis equal; axis xy; colorbar; title(['Prediction for ' Title]);
end


Compare_master = Compare_master./4;
save Compare Compare;
figure(5); imagesc(Compare_master); axis equal; axis xy; colorbar; 
%imagesc(Compare_master);
%axis equal;
%axis xy;
%colorbar;


        