function [Fx,Fy,Fz,Mx,My,Mz,sigma,tau_t,tau_l] = Force_quad_fun(att,btt,ctt, ...
                                  atl,btl,ctl, ...
                                  as,bs,cs, ...
                                  theta0,theta1,r,b)
    Fx=0;
    Fy=0;
    Fz=0;
    
    Mx=0;
    My=0;
    Mz=0;
    
    for i =1:2
        th0=theta0(i);
        th1=theta1(i);
    
        attt=att(i);  
        bttt=btt(i);
        cttt=ctt(i);
    
        atll=atl(i);
        btll=btl(i);
        ctll=ctl(i);
    
        ass=as(i);
        bss=bs(i);
        css=cs(i);

        Fx_part= @(theta) ( attt.*(theta.^2 - 2) - 2*ass.*theta + bttt.*theta - bss+ cttt ) .* (sin(theta)) + ...
                ( ass.*(theta.^2 - 2) + 2*attt.*theta + bss.*theta + bttt + css ) .* (cos( theta));
        Fx_i= Fx_part(th1)-Fx_part(th0);
        Fx= Fx + r*b*Fx_i;

        Fy_part= @(theta)  atll.*theta.^3/3 + btll.*theta.^2/2 + ctll.*theta;
        Fy_i= Fy_part(th1)-Fy_part(th0);
        Fy= Fy- r*b*Fy_i;

        Fz_part= @(theta) ( ass.*(theta.^2 - 2) + 2*attt.*theta + bss.*theta + bttt+ css ) .* (sin(theta)) + ...
                ( -attt.*(theta.^2 - 2) + 2*ass.*theta - bttt.*theta + bss- cttt ) .* (cos(theta));
        Fz_i= Fz_part(th1)-Fz_part(th0);
        Fz=  Fz +r*b*Fz_i;

        Mx_part = @(theta) ...
        ( atll.*(theta.^2 - 2) .* (sin(theta)) + 2*atll.*theta .* (cos( theta) )) + ...
        ( btll.*theta .* (sin(theta)) + btll.* (cos( theta)) + ctll .* (sin( theta)) );
        Mx_i = Mx_part(th1) - Mx_part(th0);
        Mx = Mx - r^2 * b * Mx_i;

        My_part = @(theta) ...
                attt.*theta.^3/3 + bttt.*theta.^2/2 + cttt.*theta;
        My_i = My_part(th1) - My_part(th0);
        My = My - r^2 * b * My_i;

        Mz_part = @(theta) ...
                ( -atll.*(theta.^2 - 2) .* (cos(theta)) + 2*atll.*theta .* (sin(theta)) ) + ...
                ( -btll.*theta .* (cos(theta)) + btll.*sin( theta) - ctll .* (cos(theta)));
        Mz_i = Mz_part(th1) - Mz_part(th0);
        Mz = Mz - r^2 * b * Mz_i;


        sigma = @(th) ass(i).*th.^2 + bss(i).*th + css(i);
        tau_t   = @(th) attt(i).*th.^2  + bttt(i).*th  + cttt(i);
        tau_l = @(th) atll(i).*th.^2  + btll(i).*th  + ctll(i);
        % Fx = Fx + r*b * integral(@(th) ...
        %     tau(th).*cos(th) - sigma(th).*sin(th), th0, th1);
        % 
        % Fz = Fz + r*b * integral(@(th) ...
        %     sigma(th).*cos(th) + tau(th).*sin(th), th0, th1);
        % 
        % My = My + r^2*b * integral(@(th) tau(th), th0, th1);
        % 
    
    end 

end