%% f_bSim
%
% Calculate simluated bunching, using either uniform or lognormal density

function bSim = calcBFromZ(lZ,uZ, data)

    if strcmp(data.density, 'uniform')
                
        bSim = max(uZ - lZ,0);
           
    elseif strcmp(data.density, 'lognormal')        
        
        if (min([uZ ; lZ]) < 0 || ~isreal([uZ ; lZ]))
            bSim = 1e8;
        else
            
            probSim = max(logncdf(uZ,data.mu,data.sigma)-logncdf(lZ,data.mu,data.sigma),0);
            
            % scale by density at kink
            bSim = probSim ./ lognpdf(data.zStar,data.mu,data.sigma);
            
        end
        
    elseif strcmp(data.density, 'flexible')
        bSim = calcBunchFlex(uZ, lZ, data) ./ data.densAK;
        
    else
        error('data.density must be uniform, lognormal, or flexible')
    end

end