function [ret,model] = sprop(Cell_Info_X)
%propegation model
% usage: 
%       [ret] = 
%       Cell_Info_X = input struct
%       
%

%warning off MATLAB:divideByZero

[M,N] = size(Cell_Info_X.terrain);                     %dimentions M rows(hight) by N columns(width)
azimuth = Cell_Info_X.cellSite.azimuth;             %transmitter azimuth angle azimuth=0 on y+ axis
eirp = Cell_Info_X.cellSite.eirp;                   %transmitter EIRP
BSx = Cell_Info_X.cellSite.BSx;                     %transmitter x position
BSy = Cell_Info_X.cellSite.BSy;                     %transmitter y position
freq = 1e6*Cell_Info_X.cellSite.freq;               %transmitter frequency in mhz * 1e6 Hz/MHz
cellSize = Cell_Info_X.cellSite.cellSize;           %size of each pixel in meters
AGL = Cell_Info_X.cellSite.height;                  %antenna height AGL
terrain = double(Cell_Info_X.terrain);              %terrain in meters ASL double format
ASL = AGL + terrain(BSy,BSx);                       %antenna height ASL
[FX,FY] = gradient(terrain);                        %X (L to R) and Y (T to B) gradients of terrain (slope in meters of elevation per pixel)


%Contants
c = 3e8;                 %speed of light
lambda = c/freq;

%variation of parameters data
%-----------------------------
Gr = 78;                           %Receiver gain
n = 1.86;                          %path loss exponent

gamma0 = 1;                        %azmuthial radiation pattern fourier coefficient F0
gamma1 =  0.022;                   %F1
gamma2 =  0.005;                   %F2
gamma3 = -0.017;                   %F3
gamma4 =  0.002;                   %F4
gamma5 =  0.0075;                  %F5
gamma6 = -0.002;                   %F6
gamma7 = -0.007;                   %F7
gamma8 = -0.009;                   %F8
gamma9 =  0.003;                   %F9

delta0 = 1;                        %elevation radietion pattern fourier coefficient F0
delta1 = 0.075;                    %F1
delta2 = 0.007;                    %F2
delta3 = 0.000;                    %F3
delta4 = 0.000;                    %F4

%terrain coefficients
terrv0 = -0.17;

%initilize prediction matrix
pred=zeros([M,N]);



for I = 1:N,     %(width X(to right))
    for J = 1:M, %(hight Y(top down))
        
        r=sqrt((cellSize^2)*((I-BSx)^2+(J-BSy)^2)+(ASL-terrain(M,N))^2);    %distance in meters to observer
        phi=atan2(-(J-BSy),I-BSx)-pi/2;                                     %angle with respect to phi=0 on y-axis -phi = cw, phi=-pi/2 for x, phi=pi/2 for -x
        theta=atan2(terrain(J,I)-ASL,cellSize*sqrt((I-BSx)^2+(J-BSy)^2));   %elevation angle of incident wave

        %azimuth patern factor
        az=(gamma0-gamma1*cos(azimuth-phi)-gamma2*cos(2*(azimuth-phi))-gamma3*cos(3*(azimuth-phi))-gamma4*cos(4*(azimuth-phi))-gamma5*cos(5*(azimuth-phi))-gamma6*cos(6*(azimuth-phi))-gamma7*cos(7*(azimuth-phi))-gamma8*cos(8*(azimuth-phi))-gamma9*cos(9*(azimuth-phi)) );
        el=(delta0-delta1*cos(theta)-delta2*cos(2*theta) );                 %elevation pattern factor
        
        Sy=abs(cos(phi));           %signal y component cos since phi=0 on y axis
        Sx=abs(sin(phi));           %signal x component
        %                   vert comp on x slope                          vert comp on y slope  
        terrv= ( Sx*cos(-sin(phi)*atan2(FX(J,I),30)-theta) + Sy*cos(-cos(phi)*atan2(FY(J,I),30)-theta));        %terrain diffraction factor

        point = eirp + terrv0*terrv + Gr - 20*log(4*pi/lambda) - 10*n*log(r)*az*el; %calculate recieved signal strength
        
        if point < -114,            %limit value to above -114dB
            pred(J,I)= -114;
        else
            pred(J,I)= point;
        end
        
    end
end


%subplot(2,1,1);
%imagesc(Cell_Info_X.meas);
%colorbar;
%subplot(2,1,2);
%imagesc(pred);
%TITLE('Prediction for Cell Info C');
%colorbar;

model=int8(pred);
ret=Cell_Info_X.cellSite;      %return transmitter attributes
