function [sl,sr]=Rickard(A,B)

N=200; % number of discretizations
minX=min(A(1),B(1)); % the range
maxX=max(A(4),B(4));
X=linspace(minX,maxX,N);

lowerA=mg(X,A(5:8),[0 A([9 9]) 0]);
upperA=mg(X,A(1:4));
lowerB=mg(X,B(5:8),[0 B([9 9]) 0]);
upperB=mg(X,B(1:4));

Il=zeros(1,N);
num=0;
den=0;
L=0;
for i=1:N
    if lowerB(i)<=lowerA(i)
        num=num+lowerB(i);
        den=den+upperA(i);
    elseif lowerB(i)>=upperA(i)
        num=num+lowerA(i);
        den=den+lowerA(i);
    else
        L=L+1;
        Il(L)=i;
        num=num+lowerB(i);
        den=den+upperA(i);
    end
end
Il(L+1:end)=[];
sl=num/den;
nums=[lowerA(Il); lowerB(Il)];
dens=[lowerA(Il); upperA(Il)];
sl0=1;
index=2*ones(1,L);
while sl0>sl
    sl0=sl;
    for i=1:L
        tempNum=num-nums(index(i),i);
        tempDen=den-dens(index(i),i);
        if sl>(tempNum+nums(3-index(i),i))/(tempDen+dens(3-index(i),i))
            index(i)=3-index(i);
            num=tempNum+nums(index(i),i);
            den=tempDen+dens(index(i),i);
            sl=num/den;
        end
    end
end

%% Compute Sr
Ir0=find(upperB<=lowerA);
Ir1=find(upperB>upperA);
Ir2=1:N;
Ir2([Ir0 Ir1])=[];
sr=sum([upperB([Ir0 Ir2]) upperA(Ir1)])/sum([lowerA(Ir0) upperA(Ir1) upperB(Ir2)]);

