%% Effeciency here is the size of D, very ineffecient, how to search fast?
function d=density3(P,D)
if (P<1)||(log10(P)<1)
    d=D(1,1);
elseif(P<D(2001,2))
    i=fix((log10(P)-1)*100)+1;
    pace=0.01;    %%²½³¤
    x1=log10(P)-log10(D(i,2));
    x2=pace-x1;
    d=D(i,1)*x2/pace+D(i+1,1)*x1/pace;
else
    d=D(2001,1);
end