function [Rc,Rb,Rr,Rg,Rall] = Resistence_fun(as,bs,cs,rw, ...
                                  final_rover,lun_soil, contact,theta_slope,h_f)   

% a) soil compaction resistance

Rc= 0;

for i =1:2
    th0=contact.theta0(i);
    th1=contact.theta1(i);

    ass=as(i);
    bss=bs(i);
    css=cs(i);

    Rc_part= @(theta) (- ass.*(theta.^2) + 2*ass - bss.*theta -  css ) .* (cos(theta)) + ...
        ( 2*ass.*theta +  bss) .* (sin( theta));
    Rc_i= Rc_part(th1)-Rc_part(th0);
    Rc= -Rc + rw*final_rover.bw*Rc_i;
end


% b) bulldozing resistance [Theory of bearing capacity of soil-Terzaghi]

Kc= 33.37; % (lun_soil.Nc-tan(lun_soil.phi))*cos(lun_soil.phi)^2; % modulus of cohesiom 
Kgamma=72.77; % (2*lun_soil.Ngamma/tan(lun_soil.phi)+1)*cos(lun_soil.phi)^2; % modulus of densitty of soil deformation 
alfa= acos(1-h_f/rw); % angle of approach [rad] == theta_f (entry contact angle)
lr=h_f*tan(45*pi/180-lun_soil.phi/2)^2; % distance of rupture [m]

Rb= ((final_rover.bw*sin(alfa+lun_soil.phi))/(2*sin(alfa)*cos(lun_soil.phi)))*(2*lun_soil.c*Kc*h_f+lun_soil.gamma*Kgamma*h_f^2)+...
    pi*lun_soil.gamma*lr^3*(90-lun_soil.phi*180/pi)/540+pi*lun_soil.c*lr^2/180+lun_soil.c*lr^2*tan(45*pi/180+lun_soil.phi/2);

% c) Rolling resistance
Rr= lun_soil.fr.*final_rover.Ww; % [N]

% d) Gravitational resistance !!! non so se è lo stesso teta di prima
% (controlla)
Rg= final_rover.Ww.*sind(theta_slope);
                                
% e) Obstacle resistance [Wong] 

Ro=0; % per ora lo ignoro

Rall= Rc+Rb+Rr+Rg+Ro;


end