function varargout = stepest(varargin)
    
// Estimate step response and plot of idpoly type model
// 
// Calling Sequence
// stepest(sys)
// stepData = stepest(sys,flag)
// Parameters
// sys : idpoly type polynomial 
// flag : boolean type variable,default value is false(%F)
// stepData : stores step response if the flag value is true
// Description
// stepest function estimate and plot the impulse response of idpoly type function.
// Examples
// a = [1 0.2];b = [0 0.2 0.3];
// sys = idpoly(a,b,'Ts',0.1)
// stepest(sys);
// 
// Examples
// a = [1 0.2];b = [0 0.2 0.3];
// sys = idpoly(a,b,'Ts',0.1)
// flag = %T
// stepData = stepest(sys,flag)
// Authors
// Ashutosh Kumar Bhargava  

    [lhs,rhs] = argn(0)
    // checking the number of inputs
    if rhs > 2  then
        error(msprintf(gettext("%s: Unexpected number of input arguments "),"stepest"))
    end
    modelData = varargin(1)
    if typeof(modelData) <> "idpoly" then
        error(msprintf(gettext("%s: Plant model must be ""idpoly"" type. "),"stepest"))
    end
    // adding noise
    if rhs == 2 then
        noiseFlag = varargin(2)
        if typeof(noiseFlag) <> 'boolean' then
            error(msprintf(gettext("%s: Last input data must be ""boolean"".type "),"stepest"))
        end
    else
        noiseFlag = %F
    end
    z = poly(0,'z')
    aPoly = poly(modelData.a(length(modelData.a):-1:1),'z','coeff')
    bPoly = poly(modelData.b,'z','coeff')
    fPoly = poly(modelData.f(length(modelData.f):-1:1),'z','coeff')
    afPoly = aPoly*fPoly
    
    bCoeff = modelData.b
    extra = 1
    if ~bCoeff(1,1) then
        afCoeff = coeff(afPoly)
        bLength = length(bCoeff);afLength = length(afCoeff)
        if bLength == afLength then
            extra = 1
        else
            extra = z^-(bLength-afLength)
        end
    end
    
    bPoly = poly(modelData.b(length(modelData.b):-1:1),'z','coeff')
    sys = syslin('d',bPoly,afPoly)*extra
    if size(find(gsort(abs(roots(sys.den)))>1),'*') then
        tempRoot = clean(roots(sys.den))
        [sorted index] = gsort(abs(real(tempRoot)))
        tempRoot = tempRoot(index)
        n = 26/log10(abs(real(tempRoot(1))))
         
    else
        finalValue = horner(sys,1)*1
        n = 10
        tempValue = sum(ldiv(sys.num,sys.den,n))
        if finalValue >= 0 then
            while finalValue*0.99 >= tempValue || n >=1000
                n = n+10
                tempValue = sum(ldiv(sys.num,sys.den,n))
            end
        elseif finalValue < 0 then
            while finalValue*0.99 <tempValue || n >=1000
                n = n+10
                tempValue = sum(ldiv(sys.num,sys.den,n))
            end  
        end
    end
    uData = ones(1,n)
    yData = flts(uData,sys)
    timeData = (0:(n-1))*modelData.Ts
    if noiseFlag then
        varargout(1) = yData'
    else
        plot2d2(timeData,yData,2)
        varargout(1) = []
    end
endfunction