function ic = bjt(vce,ib,beta,vcesat)

%  This function generates the collector current values "ic" for
%  a bjt Transistor as a function of the collector-emitter voltages "vce".
%
%    ic = bjt(vce,ib,beta,vcesat)
%
%    where "vce" is a vector of collector-emitter values
%                "ib" is the base-current
%                "beta" forward current gain (ic/ib)
%                "vcesat" is the vce value for the boundary between the
%                saturation and active regions (usually .2V)
%
%   Written by Kevin D. Donohue (donohue@engr.uky.edu) Nov. 2005

ic = zeros(size(vce));  %  Initialize output array with all zeros


%  For non-cutoff operation:    
        %  Find points in vds that are in the saturation region
        ktri = find(vce<vcesat & vce >= 0);  %  Only take point up to the gate voltage minus the threshold.
        %  If points are found in the triode region compute ids with proper formula 
        % approximate saturation regions with simple line
        m = beta*ib/vcesat;  %  Slope
        if ~isempty(ktri)
            ic(ktri) = m*vce(ktri);
        end
        % Find points in active region
        ksat = find(vce>=vcesat);   %  Take points greater than the excess voltage
        %  if points are found in the saturation regions compute ids with proper formula 
        if ~isempty(ksat)
            ic(ksat) = beta*ib;
        end
    %  If points of vds are outside these ranges then the ids values remain zero
end