function [sin_out, cos_out] = cordicsincos(theta, niters)
% [y, x] = CORDICCSINCOS(theta, niters) Compute the CORDIC-based approximation to [SIN(theta), COS(theta)].

% Copyright 2009-2014 The MathWorks, Inc.
%#codegen

coder.allowpcode('plain');

if eml_ambiguous_types
    sin_out = coder.ignoreConst(reshape(zeros(size(theta)),size(theta)));
    cos_out = coder.ignoreConst(reshape(zeros(size(theta)),size(theta)));
else
    coder.internal.prefer_const(theta);

    % Input argument checking
    eml_cordic_check_theta_arg(theta, 'cordicsincos');

    % Input-output dimensions
    sz            = size(theta);
    needToReshape = ~isvector(theta) || (length(sz) > 2);

    if isfi(theta) && isfloat(theta)
        % Call CORDICSINCOS with the equivalent builtin floating-point type
        if issingle(theta)
            theta_float = single(theta);
        else
            theta_float = double(theta);
        end
        if nargin > 1
            [sin_out_flt, cos_out_flt] = cordicsincos(theta_float, niters);
        else
            [sin_out_flt, cos_out_flt] = cordicsincos(theta_float);
        end

        % Allocate memory for FI Double/Single outputs
        NT_FI_Float = numerictype(theta);
        sin_out = coder.nullcopy(fi(zeros(sz), NT_FI_Float));
        cos_out = coder.nullcopy(fi(zeros(sz), NT_FI_Float));
        
        % Cast builtin values to FI values
        for idx=1:length(theta)
            sin_out(idx) = sin_out_flt(idx);
        end
        for idx=1:length(theta)
            cos_out(idx) = cos_out_flt(idx);
        end
        return; % EARLY RETURN
    end
    
    % Compute maximum number of iterations constant
    if isa(theta, 'double')
        maxNITERS = 52;
    elseif isa(theta, 'single')
        maxNITERS = 23;
    elseif isfi(theta)
        thetaFPNT = numerictype(theta);
        maxNITERS = thetaFPNT.WordLength - 1;
    else
        % Builtin integer THETA input
        switch class(theta)
          case { 'int8',  'uint8'}
            maxNITERS = 7;
          case {'int16', 'uint16'}
            maxNITERS = 15;
          case {'int32', 'uint32'}
            maxNITERS = 31;
          otherwise
            maxNITERS = 63;
        end
    end
        
    if nargin > 1
        % NITERS is specified
        coder.internal.prefer_const(niters);
        eml_cordic_check_niters_arg(niters);
        if ~isempty(niters) && isfinite(niters)
            if (double(niters) <= maxNITERS)
                numIters_dbl = double(niters); % use specified constant
            else
                numIters_dbl = maxNITERS; % use maximum value
            end
        else
            % NITERS empty, inf, or nan
            numIters_dbl = maxNITERS; % use maximum value
        end
    else
        % NITERS not specified
        numIters_dbl = maxNITERS; % use maximum value
    end
    
    % Compute "real-world-values" for constant initializations
    intArray_dbl = (0:(numIters_dbl-1))';
    inputLUT_dbl = atan(2 .^ (-intArray_dbl));
    preNrmIC_dbl = 1 / fixed.internal.cordic_compute_gain(numIters_dbl);
    
    % Initialize constants
    if isa(theta, 'double')
        needTempXYOut = needToReshape;
        lut           = inputLUT_dbl;
        K             = preNrmIC_dbl;
        x0            = double(K);
        y0            = 0.0;
        if needTempXYOut
            xout = coder.nullcopy(zeros(sz));
            yout = coder.nullcopy(zeros(sz));
        end
        cos_out = coder.nullcopy(zeros(sz));
        sin_out = coder.nullcopy(zeros(sz));
    elseif isa(theta, 'single')
        needTempXYOut = needToReshape;
        lut           = single(inputLUT_dbl);
        K             = single(preNrmIC_dbl);
        x0            = single(K);
        y0            = single(0.0);
        if needTempXYOut
            xout = coder.nullcopy(single(zeros(sz)));
            yout = coder.nullcopy(single(zeros(sz)));
        end
        cos_out = coder.nullcopy(single(zeros(sz)));
        sin_out = coder.nullcopy(single(zeros(sz)));
    else
        % Fixed-point or integer THETA
        if isfi(theta)
            thetaNT      = numerictype(theta);
            ioWordLength = thetaNT.WordLength;
        else
            % Builtin integer THETA input
            switch class(theta)
              case { 'int8',  'uint8'}
                ioWordLength =  8;
              case {'int16', 'uint16'}
                ioWordLength = 16;
              case {'int32', 'uint32'}
                ioWordLength = 32;
              otherwise
                ioWordLength = 64;
            end
        end
        
        eml_invariant(ioWordLength > 1, ...
            'fixed:cordic:inputWordLengthNotGTOne');
    
        needTempXYOut = true; % Due to use of LOCAL FIMATH for CORDIC internals
        ioFracLength  = ioWordLength - 2;
        if isfi(theta)
            ioNT = numerictype(numerictype(1, ioWordLength, ioFracLength),'DataType',thetaNT.DataType);
        else
            ioNT = numerictype(1, ioWordLength, ioFracLength);
        end
        
        % First initialize all values using the "fimathless FI" rules
        % (i.e. cast "real-world-values" using NEAREST and SATURATE)
        lutTmp = fi(inputLUT_dbl, ioNT);
        kTmp   = fi(preNrmIC_dbl, ioNT);
        
        % Use same LOCAL FIMATH for CORDIC arithmetic
        Fm      = eml_al_cordic_fimath(theta);
        lut     = fi(lutTmp, ioNT, Fm);
        K       = fi(kTmp,   ioNT, Fm);
        x0      = fi(K,      ioNT, Fm);
        y0      = fi(0,      ioNT, Fm);
        xout    = coder.nullcopy(fi(zeros(sz), ioNT, Fm));
        yout    = coder.nullcopy(fi(zeros(sz), ioNT, Fm));
        cos_out = coder.nullcopy(fi(zeros(sz), ioNT)); % NOT using local Fm
        sin_out = coder.nullcopy(fi(zeros(sz), ioNT)); % NOT using local Fm
    end
    
    % Type for numerical operation and intermediate results
    if isfloat(theta)
        lutValues = lut;
    else
        lutValues = fi(lut, ioNT, Fm);
    end
    
    if needToReshape
        % FULL MATRIX or N-D cases: Form column vector and reshape afterward.
        % (NOTE: this will likely generate extra input and output copy code...)
        angle = theta(:);
        
        for idx=1:length(angle)
            % Correct input angle to [-pi/2, pi/2]
            if isfloat(angle)
                [negate, z0] = eml_al_cordic_quad_correction_before_float(angle(idx));
            else
                [negate, z0] = eml_al_cordic_quad_correction_before(angle(idx), K);
            end
            
            % CORDIC computing
            [xn, yn] = fixed.internal.cordic_rotation_kernel_private(x0, y0, z0, lutValues, numIters_dbl);
            
            % Restore signs of output sin and cos
            % (NOTE: using temp xout, yout since needToReshape == TRUE)
            [xout(idx), yout(idx)] = eml_al_cordic_quad_correction_after(xn, yn, negate);
        end
    else
        % No input-copy/reshape required
        for idx=1:length(theta)
            % Correct input angle to [-pi/2, pi/2]
            if isfloat(theta)
                [negate, z0] = eml_al_cordic_quad_correction_before_float(theta(idx));
            else
                [negate, z0] = eml_al_cordic_quad_correction_before(theta(idx), K);
            end
            % CORDIC computing
            [xn, yn] = fixed.internal.cordic_rotation_kernel_private(x0, y0, z0, lutValues, numIters_dbl);

            % Restore signs of output sin and cos
            if needTempXYOut
                % Using temp variables, e.g., for LOCAL FIMATH.
                % (NOTE: may generate extra output copy code below...)
                [xout(idx), yout(idx)] = eml_al_cordic_quad_correction_after(xn, yn, negate);
            else
                % No extra output copies required
                [cos_out(idx), sin_out(idx)] = eml_al_cordic_quad_correction_after(xn, yn, negate);
            end
        end
    end
    
    if needTempXYOut
        % Copy temp XOUT and YOUT vars to COS_OUT and SIN_OUT outputs
        [cos_out, sin_out] = ...
            eml_cordic_copy_reshape_xy(xout, yout, needToReshape, sz);
    end
end % if eml_ambiguous_types

end % function

