X = load('Cell_Info.mat');
Compare_master = load('Compare.mat');
Names = fieldnames(X);
close all;
Compare_master = Compare_master.Compare;
n=3;
c = 3E8;
for cell= 1:length(Names),
%for cell= 1:1,
    Cell_Info = getfield(X,Names{cell});
    freq = Cell_Info.freq;
    f = (freq)*10^6;
    APL=12;
    lambda = c/f;
  %  Image_dims = size(A);
    %Anew = zeros(Image_dims);
    lambdaloss = 20*log10((lambda)/(4*pi));


my_raster=Cell_Info;
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

for i=-length(Cell_Info.ReceivedPowerMap):length(Cell_Info.ReceivedPowerMap),
    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 = [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 = my_raster.ReceivedPowerMap(1:203, 1:203);
if cell == 1
    Cell_Info_A = Cell_Info;
    my_raster.ReceivedPowerMap = my_raster.ReceivedPowerMap - Compare_master;
end

if cell >1
    Image_dims = size(Compare_master);    
    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_master(im_x,im_y);
            end
        end
    end 
  
  my_raster.ReceivedPowerMap = my_raster.ReceivedPowerMap - newCompare;
end

if cell == 1
    Cell_Info_A = my_raster;
end
if cell == 2
    Cell_Info_B = my_raster;
end
if cell == 3
    Cell_Info_C = my_raster;
end
if cell == 4
    Cell_Info_D = my_raster;
end
if cell == 5
    Cell_Info_E = my_raster;
end
if cell == 6
    Cell_Info_F = my_raster;
end
if cell == 7
    Cell_Info_G = my_raster;
end
if cell == 8
    Cell_Info_H = my_raster;
end

        %Title = Names{cell}; Title( find(Title=='_') ) = ' ';
        %figure(cell); imagesc(my_raster.ReceivedPowerMap); axis equal; axis xy; colorbar; title(['Prediction for ' Title]);
end
save test4 Cell_Info_A Cell_Info_B Cell_Info_C Cell_Info_D;
save Berntsson_Leif Cell_Info_E Cell_Info_F Cell_Info_G Cell_Info_H;

%imagesc(my_raster.ReceivedPowerMap);
%axis equal;
%axis xy;
%colorbar;
