function [D] = D(flag,Kf,K0,theta1)
if Kf==0
   D = 1; 
else
    if flag==1
        D = (1/Kf-1/K0)*(theta1+1/Kf)^-1;
        if isnan(D) || isinf(D)
            disp('D IS NAN OR INF -> SETTING TO ZERO!')
            D = 0;
        end
    else 
        D = 0;
    end
end