function friction = friction_fun(sFLx, sFLy, sFL, sFRx, sFRy, sFR, sRLx, sRLy, sRL, sRRx, sRRy, sRR, z, zdot, theta, thetadot, phi, phidot)

global B C D W lf lr h tf tr Kf Cf Kr Cr


%FRONT LEFT
if sFL==0;
    muFL = 0; muFLx = 0; muFLy = 0;
else
    muFL = D*sin(C*atan(B*sFL));
    muFLx = -(sFLx/sFL)*muFL;
    muFLy = -(sFLy/sFL)*muFL;
end

%FRONT RIGHT
if sFR==0;
    muFR = 0; muFRx = 0; muFRy = 0;
else
    muFR = D*sin(C*atan(B*sFR));
    muFRx = -(sFRx/sFR)*muFR;
    muFRy = -(sFRy/sFR)*muFR;
end

%REAR LEFT
if sRL==0;
    muRL = 0; muRLx = 0; muRLy = 0;
else
    muRL = D*sin(C*atan(B*sRL));
    muRLx = -(sRLx/sRL)*muRL;
    muRLy = -(sRLy/sRL)*muRL;
end

%REAR RIGHT
if sRR==0;
    muRR = 0; muRRx = 0; muRRy = 0;
else
    muRR = D*sin(C*atan(B*sRR));
    muRRx = -(sRRx/sRR)*muRR;
    muRRy = -(sRRy/sRR)*muRR;
end


%CALCULATION OF LOAD (f_(i,j)z)

%AZ on each wheel
AZFL = (z-h) + lf*sin(theta) + tf*0.5*sin(phi);
AZFR = (z-h) + lf*sin(theta) - tf*0.5*sin(phi);
AZRL = (z-h) - lr*sin(theta) + tr*0.5*sin(phi);
AZRR = (z-h) - lr*sin(theta) - tr*0.5*sin(phi);

% AZdot on each wheel
AZFLdot = zdot + lf*thetadot*cos(theta) + tf*0.5*phidot*cos(phi);
AZFRdot = zdot + lf*thetadot*cos(theta) - tf*0.5*phidot*cos(phi);
AZRLdot = zdot - lr*thetadot*cos(theta) + tr*0.5*phidot*cos(phi);
AZRRdot = zdot - lr*thetadot*cos(theta) - tr*0.5*phidot*cos(phi);

%fFi,j Z
fFLz =  0.5*(W*lr/(lf+lr)) - AZFL*Kf - AZFLdot*Cf;
fFRz =  0.5*(W*lr/(lf+lr)) - AZFR*Kf - AZFRdot*Cf;
fRLz =  0.5*(W*lf/(lf+lr)) - AZRL*Kr - AZRLdot*Cr;
fRRz =  0.5*(W*lf/(lf+lr)) - AZRR*Kr - AZRRdot*Cr;

%Forces fij, x,y on each wheel
fFL = fFLz*muFL; fFLx = fFLz*muFLx; fFLy = fFLz*muFLy;
fFR = fFRz*muFR; fFRx = fFRz*muFRx; fFRy = fFRz*muFRy;

fRL = fRLz*muRL; fRLx = fRLz*muRLx; fRLy = fRLz*muRLy;
fRR = fRRz*muRR; fRRx = fRRz*muRRx; fRRy = fRRz*muRRy;

friction = [fFLx fFLy fFL fFRx fFRy fFR fRLx fRLy fRL fRRx fRRy fRR fFLz fFRz fRLz fRRz]';