function hd=ExoterDB(M0,sigmaMun,R0,sigmaRun)
warning off all
%%Constant:
Mearth=5.9736e24;
Rearth=6372797;
load Newdata3;
%%step:
step=0.002;

%find how much mass uncertainty is equivalent to radius uncertainty in the
%sense of composition
xx=[0.5,1,2,4,8,16,20];
yy=[3.3,3.4,3.6,3.9,4.1,4.4,4.5];
ratio=interp1(xx,yy,M0,'pchip','extrap');
correction=sqrt((sigmaMun/M0)^2+(ratio*(sigmaRun/R0))^2)/(abs(sigmaMun/M0)+ratio*abs(sigmaRun/R0));
sigmaM=correction*sigmaMun;
sigmaR=correction*sigmaRun;


%%create the picture
figure;
terplot(11);
terlabel('\beta (silicate)','\alpha (iron)','\gamma (water)');
title(['Mass:',num2str(M0),'M_\oplus',' \sigma_M:',num2str(sigmaMun),'M_\oplus',...
       '   Radius:',num2str(R0),'R_\oplus',' \sigma_R:', num2str(sigmaRun),'R_\oplus'],'FontSize',14,'FontWeight','bold');
hold on

%%Mass, Radius and their uncertainties we know:
w=0;

for v=-3:3
m0=M0*(1+v*sigmaM/M0);
r0=R0*(1-v*sigmaR/R0);

%%control parameters:
j=0;
k=0;
a=0;

while ((a>=0)&(a<=1))&(j*k<0.5)
    [avector,bvector,mvector]=meshgrid(a,0:step:(1-a),m0);
    Radiusesti=interp3(x,y,z,Radius,avector,bvector,mvector);
    if ((max(Radiusesti)>=r0) && (min(Radiusesti)<=r0))
        [f,g]=min(abs(Radiusesti-r0));
        w=w+1;
        iron(w)=a;
        silicate(w)=bvector(g);
        err(w)=abs(v);
        j=1;k=0;
    else
        k=1;
    end
    a=a+step;
end
end
hd=Exocontour(iron,silicate,err);
hold off;