x_coord = Cell_Info_D.cellSite.BSx;
y_coord = Cell_Info_D.cellSite.BSy;
cheight = Cell_Info_D.cellSite.height + Cell_Info_D.terrain(y_coord,x_coord);
termap = Cell_Info_D.terrain;
lambda = (3*10^8)/(Cell_Info_D.cellSite.freq * 10^6);

% Create Map and take mountain shadowing into account
s = MakeShadowMap(x_coord, y_coord, cheight, Cell_Info_D.terrain);
backup = s;

[n,m] = size(s);
for z = 1:n
    for v = 1:m
        ydist = 30*(z-y_coord);
        xdist = 30*(v-x_coord);
        dist  = sqrt(xdist^2 + ydist^2);

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %Determine which points lie in 110 degree directivity region
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  
	  angl = (atan2(ydist,xdist)*180/pi);
        if (((angl >= -90) & (angl <=-55)) | ((angl < -90) & (angl >= -125)))
            s(z,v) = s(z,v) + 20;
        elseif (((angl >-55) & (angl < -35)) | ((angl <=-125) & (angl >-145)))
            s(z,v) = s(z,v) + 10;
        end
        s(z,v) = -40 + s(z,v) - 20*log10(dist);

	  % -120 is minimum gain
        if (s(z,v) < -120)
            s(z,v) = -120;
        end
    end
end
imagesc(s)
axis xy
figure
imagesc(Cell_Info_D.meas)
axis xy


