/*
 * MATLAB Compiler: 3.0
 * Date: Mon Nov 03 13:00:57 2003
 * Arguments: "-B" "macro_default" "-O" "all" "-O" "fold_scalar_mxarrays:on"
 * "-O" "fold_non_scalar_mxarrays:on" "-O" "optimize_integer_for_loops:on" "-O"
 * "array_indexing:on" "-O" "optimize_conditionals:on" "-m" "-W" "main" "-L"
 * "C" "-t" "-T" "link:exe" "-h" "libmmfile.mlib" "DancingMonkeys" 
 */
#include "bilinear.h"
#include "libmatlbm.h"
#include "libmmfile.h"
static mxArray * _mxarray0_;

static mxChar _array2_[50] = { 'N', 'u', 'm', 'e', 'r', 'a', 't', 'o', 'r',
                               ' ', 'c', 'a', 'n', 'n', 'o', 't', ' ', 'b',
                               'e', ' ', 'h', 'i', 'g', 'h', 'e', 'r', ' ',
                               'o', 'r', 'd', 'e', 'r', ' ', 't', 'h', 'a',
                               'n', ' ', 'd', 'e', 'n', 'o', 'm', 'i', 'n',
                               'a', 't', 'o', 'r', '.' };
static mxArray * _mxarray1_;
static mxArray * _mxarray3_;
static mxArray * _mxarray4_;

static mxChar _array6_[51] = { 'F', 'i', 'r', 's', 't', ' ', 't', 'w', 'o',
                               ' ', 'a', 'r', 'g', 'u', 'm', 'e', 'n', 't',
                               's', ' ', 'm', 'u', 's', 't', ' ', 'h', 'a',
                               'v', 'e', ' ', 't', 'h', 'e', ' ', 's', 'a',
                               'm', 'e', ' ', 'o', 'r', 'i', 'e', 'n', 't',
                               'a', 't', 'i', 'o', 'n', '.' };
static mxArray * _mxarray5_;

void InitializeModule_bilinear(void) {
    _mxarray0_ = mclInitializeDouble(2.0);
    _mxarray1_ = mclInitializeString(50, _array2_);
    _mxarray3_ = mclInitializeDouble(6.283185307179586);
    _mxarray4_ = mclInitializeDouble(1.0);
    _mxarray5_ = mclInitializeString(51, _array6_);
}

void TerminateModule_bilinear(void) {
    mxDestroyArray(_mxarray5_);
    mxDestroyArray(_mxarray4_);
    mxDestroyArray(_mxarray3_);
    mxDestroyArray(_mxarray1_);
    mxDestroyArray(_mxarray0_);
}

static mxArray * Mbilinear(mxArray * * pd,
                           mxArray * * kd,
                           mxArray * * dd,
                           int nargout_,
                           mxArray * z,
                           mxArray * p,
                           mxArray * k,
                           mxArray * fs,
                           mxArray * fp,
                           mxArray * fp1);

_mexLocalFunctionTable _local_function_table_bilinear
  = { 0, (mexFunctionTableEntry *)NULL };

/*
 * The function "mlfNBilinear" contains the nargout interface for the
 * "bilinear" M-function from file
 * "c:\matlab65\toolbox\signal\signal\bilinear.m" (lines 1-102). This interface
 * is only produced if the M-function uses the special variable "nargout". The
 * nargout interface allows the number of requested outputs to be specified via
 * the nargout argument, as opposed to the normal interface which dynamically
 * calculates the number of outputs based on the number of non-NULL inputs it
 * receives. This function processes any input arguments and passes them to the
 * implementation version of the function, appearing above.
 */
mxArray * mlfNBilinear(int nargout,
                       mxArray * * pd,
                       mxArray * * kd,
                       mxArray * * dd,
                       mxArray * z,
                       mxArray * p,
                       mxArray * k,
                       mxArray * fs,
                       mxArray * fp,
                       mxArray * fp1) {
    mxArray * zd = NULL;
    mxArray * pd__ = NULL;
    mxArray * kd__ = NULL;
    mxArray * dd__ = NULL;
    mlfEnterNewContext(3, 6, pd, kd, dd, z, p, k, fs, fp, fp1);
    zd = Mbilinear(&pd__, &kd__, &dd__, nargout, z, p, k, fs, fp, fp1);
    mlfRestorePreviousContext(3, 6, pd, kd, dd, z, p, k, fs, fp, fp1);
    if (pd != NULL) {
        mclCopyOutputArg(pd, pd__);
    } else {
        mxDestroyArray(pd__);
    }
    if (kd != NULL) {
        mclCopyOutputArg(kd, kd__);
    } else {
        mxDestroyArray(kd__);
    }
    if (dd != NULL) {
        mclCopyOutputArg(dd, dd__);
    } else {
        mxDestroyArray(dd__);
    }
    return mlfReturnValue(zd);
}

/*
 * The function "mlfBilinear" contains the normal interface for the "bilinear"
 * M-function from file "c:\matlab65\toolbox\signal\signal\bilinear.m" (lines
 * 1-102). This function processes any input arguments and passes them to the
 * implementation version of the function, appearing above.
 */
mxArray * mlfBilinear(mxArray * * pd,
                      mxArray * * kd,
                      mxArray * * dd,
                      mxArray * z,
                      mxArray * p,
                      mxArray * k,
                      mxArray * fs,
                      mxArray * fp,
                      mxArray * fp1) {
    int nargout = 1;
    mxArray * zd = NULL;
    mxArray * pd__ = NULL;
    mxArray * kd__ = NULL;
    mxArray * dd__ = NULL;
    mlfEnterNewContext(3, 6, pd, kd, dd, z, p, k, fs, fp, fp1);
    if (pd != NULL) {
        ++nargout;
    }
    if (kd != NULL) {
        ++nargout;
    }
    if (dd != NULL) {
        ++nargout;
    }
    zd = Mbilinear(&pd__, &kd__, &dd__, nargout, z, p, k, fs, fp, fp1);
    mlfRestorePreviousContext(3, 6, pd, kd, dd, z, p, k, fs, fp, fp1);
    if (pd != NULL) {
        mclCopyOutputArg(pd, pd__);
    } else {
        mxDestroyArray(pd__);
    }
    if (kd != NULL) {
        mclCopyOutputArg(kd, kd__);
    } else {
        mxDestroyArray(kd__);
    }
    if (dd != NULL) {
        mclCopyOutputArg(dd, dd__);
    } else {
        mxDestroyArray(dd__);
    }
    return mlfReturnValue(zd);
}

/*
 * The function "mlfVBilinear" contains the void interface for the "bilinear"
 * M-function from file "c:\matlab65\toolbox\signal\signal\bilinear.m" (lines
 * 1-102). The void interface is only produced if the M-function uses the
 * special variable "nargout", and has at least one output. The void interface
 * function specifies zero output arguments to the implementation version of
 * the function, and in the event that the implementation version still returns
 * an output (which, in MATLAB, would be assigned to the "ans" variable), it
 * deallocates the output. This function processes any input arguments and
 * passes them to the implementation version of the function, appearing above.
 */
void mlfVBilinear(mxArray * z,
                  mxArray * p,
                  mxArray * k,
                  mxArray * fs,
                  mxArray * fp,
                  mxArray * fp1) {
    mxArray * zd = NULL;
    mxArray * pd = NULL;
    mxArray * kd = NULL;
    mxArray * dd = NULL;
    mlfEnterNewContext(0, 6, z, p, k, fs, fp, fp1);
    zd = Mbilinear(&pd, &kd, &dd, 0, z, p, k, fs, fp, fp1);
    mlfRestorePreviousContext(0, 6, z, p, k, fs, fp, fp1);
    mxDestroyArray(zd);
    mxDestroyArray(pd);
    mxDestroyArray(kd);
    mxDestroyArray(dd);
}

/*
 * The function "mlxBilinear" contains the feval interface for the "bilinear"
 * M-function from file "c:\matlab65\toolbox\signal\signal\bilinear.m" (lines
 * 1-102). The feval function calls the implementation version of bilinear
 * through this function. This function processes any input arguments and
 * passes them to the implementation version of the function, appearing above.
 */
void mlxBilinear(int nlhs, mxArray * plhs[], int nrhs, mxArray * prhs[]) {
    mxArray * mprhs[6];
    mxArray * mplhs[4];
    int i;
    if (nlhs > 4) {
        mlfError(
          mxCreateString(
            "Run-time Error: File: bilinear Line: 1 Column:"
            " 1 The function \"bilinear\" was called with m"
            "ore than the declared number of outputs (4)."),
          NULL);
    }
    if (nrhs > 6) {
        mlfError(
          mxCreateString(
            "Run-time Error: File: bilinear Line: 1 Column:"
            " 1 The function \"bilinear\" was called with m"
            "ore than the declared number of inputs (6)."),
          NULL);
    }
    for (i = 0; i < 4; ++i) {
        mplhs[i] = NULL;
    }
    for (i = 0; i < 6 && i < nrhs; ++i) {
        mprhs[i] = prhs[i];
    }
    for (; i < 6; ++i) {
        mprhs[i] = NULL;
    }
    mlfEnterNewContext(
      0, 6, mprhs[0], mprhs[1], mprhs[2], mprhs[3], mprhs[4], mprhs[5]);
    mplhs[0]
      = Mbilinear(
          &mplhs[1],
          &mplhs[2],
          &mplhs[3],
          nlhs,
          mprhs[0],
          mprhs[1],
          mprhs[2],
          mprhs[3],
          mprhs[4],
          mprhs[5]);
    mlfRestorePreviousContext(
      0, 6, mprhs[0], mprhs[1], mprhs[2], mprhs[3], mprhs[4], mprhs[5]);
    plhs[0] = mplhs[0];
    for (i = 1; i < 4 && i < nlhs; ++i) {
        plhs[i] = mplhs[i];
    }
    for (; i < 4; ++i) {
        mxDestroyArray(mplhs[i]);
    }
}

/*
 * The function "Mbilinear" is the implementation version of the "bilinear"
 * M-function from file "c:\matlab65\toolbox\signal\signal\bilinear.m" (lines
 * 1-102). It contains the actual compiled code for that M-function. It is a
 * static function and must only be called from one of the interface functions,
 * appearing below.
 */
/*
 * function [zd, pd, kd, dd] = bilinear(z, p, k, fs, fp, fp1)
 */
static mxArray * Mbilinear(mxArray * * pd,
                           mxArray * * kd,
                           mxArray * * dd,
                           int nargout_,
                           mxArray * z,
                           mxArray * p,
                           mxArray * k,
                           mxArray * fs,
                           mxArray * fp,
                           mxArray * fp1) {
    mexLocalFunctionTable save_local_function_table_
      = mclSetCurrentLocalFunctionTable(&_local_function_table_bilinear);
    int nargin_ = mclNargin(6, z, p, k, fs, fp, fp1, NULL);
    mxArray * zd = NULL;
    mxArray * cd = NULL;
    mxArray * bd = NULL;
    mxArray * ad = NULL;
    mxArray * t2 = NULL;
    mxArray * t1 = NULL;
    mxArray * r = NULL;
    mxArray * t = NULL;
    mxArray * den = NULL;
    mxArray * num = NULL;
    mxArray * d = NULL;
    mxArray * c = NULL;
    mxArray * b = NULL;
    mxArray * a = NULL;
    mxArray * ans = NULL;
    mxArray * nd = NULL;
    mxArray * md = NULL;
    mxArray * nn = NULL;
    mxArray * mn = NULL;
    mclCopyArray(&z);
    mclCopyArray(&p);
    mclCopyArray(&k);
    mclCopyArray(&fs);
    mclCopyArray(&fp);
    mclCopyArray(&fp1);
    /*
     * %BILINEAR Bilinear transformation with optional frequency prewarping.
     * %   [Zd,Pd,Kd] = BILINEAR(Z,P,K,Fs) converts the s-domain transfer
     * %   function specified by Z, P, and K to a z-transform discrete
     * %   equivalent obtained from the bilinear transformation:
     * %
     * %      H(z) = H(s) |
     * %                  | s = 2*Fs*(z-1)/(z+1)
     * %
     * %   where column vectors Z and P specify the zeros and poles, scalar
     * %   K specifies the gain, and Fs is the sample frequency in Hz.
     * %
     * %   [NUMd,DENd] = BILINEAR(NUM,DEN,Fs), where NUM and DEN are 
     * %   row vectors containing numerator and denominator transfer
     * %   function coefficients, NUM(s)/DEN(s), in descending powers of
     * %   s, transforms to z-transform coefficients NUMd(z)/DENd(z).
     * %
     * %   [Ad,Bd,Cd,Dd] = BILINEAR(A,B,C,D,Fs) is a state-space version.
     * %
     * %   Each of the above three forms of BILINEAR accepts an optional
     * %   additional input argument that specifies prewarping. 
     * %
     * %   For example, [Zd,Pd,Kd] = BILINEAR(Z,P,K,Fs,Fp) applies prewarping 
     * %   before the bilinear transformation so that the frequency responses
     * %   before and after mapping match exactly at frequency point Fp
     * %   (match point Fp is specified in Hz).
     * %
     * %   See also IMPINVAR.
     * 
     * %   Author(s): J.N. Little, 4-28-87 
     * %   	   J.N. Little, 5-5-87, revised
     * %   Copyright 1988-2002 The MathWorks, Inc.
     * %   $Revision: 1.7 $  $Date: 2002/03/28 17:27:00 $
     * 
     * %   Gene Franklin, Stanford Univ., motivated the state-space
     * %   approach to the bilinear transformation.
     * 
     * [mn,nn] = size(z);
     */
    mlfSize(mlfVarargout(&mn, &nn, NULL), mclVa(z, "z"), NULL);
    /*
     * [md,nd] = size(p);
     */
    mlfSize(mlfVarargout(&md, &nd, NULL), mclVa(p, "p"), NULL);
    /*
     * 
     * if (nd == 1 & nn < 2) & nargout ~= 4	% In zero-pole-gain form
     */
    {
        mxArray * a_ = mclInitialize(mclEq(mclVv(nd, "nd"), _mxarray4_));
        if (mlfTobool(a_)) {
            mlfAssign(&a_, mclAnd(a_, mclLt(mclVv(nn, "nn"), _mxarray0_)));
        } else {
            mlfAssign(&a_, mlfScalar(0));
        }
        if (mlfTobool(a_)
            && mlfTobool(mclAnd(a_, mclBoolToArray(nargout_ != 4)))) {
            mxDestroyArray(a_);
            /*
             * if mn > md
             */
            if (mclGtBool(mclVv(mn, "mn"), mclVv(md, "md"))) {
                /*
                 * error('Numerator cannot be higher order than denominator.')
                 */
                mlfError(_mxarray1_, NULL);
            /*
             * end
             */
            }
            /*
             * if nargin == 5		% Prewarp
             */
            if (nargin_ == 5) {
                /*
                 * fp = 2*pi*fp;
                 */
                mlfAssign(&fp, mclMtimes(_mxarray3_, mclVa(fp, "fp")));
                /*
                 * fs = fp/tan(fp/fs/2);
                 */
                mlfAssign(
                  &fs,
                  mclMrdivide(
                    mclVa(fp, "fp"),
                    mlfTan(
                      mclMrdivide(
                        mclMrdivide(mclVa(fp, "fp"), mclVa(fs, "fs")),
                        _mxarray0_))));
            /*
             * else
             */
            } else {
                /*
                 * fs = 2*fs;
                 */
                mlfAssign(&fs, mclMtimes(_mxarray0_, mclVa(fs, "fs")));
            /*
             * end
             */
            }
            /*
             * z = z(finite(z));	 % Strip infinities from zeros
             */
            mlfAssign(
              &z, mclArrayRef1(mclVa(z, "z"), mlfFinite(mclVa(z, "z"))));
            /*
             * pd = (1+p/fs)./(1-p/fs); % Do bilinear transformation
             */
            mlfAssign(
              pd,
              mclRdivide(
                mclPlus(
                  _mxarray4_, mclMrdivide(mclVa(p, "p"), mclVa(fs, "fs"))),
                mclMinus(
                  _mxarray4_, mclMrdivide(mclVa(p, "p"), mclVa(fs, "fs")))));
            /*
             * zd = (1+z/fs)./(1-z/fs);
             */
            mlfAssign(
              &zd,
              mclRdivide(
                mclPlus(
                  _mxarray4_, mclMrdivide(mclVa(z, "z"), mclVa(fs, "fs"))),
                mclMinus(
                  _mxarray4_, mclMrdivide(mclVa(z, "z"), mclVa(fs, "fs")))));
            /*
             * % real(kd) or just kd?
             * kd = (k*prod(fs-z)./prod(fs-p));
             */
            mlfAssign(
              kd,
              mclRdivide(
                mclMtimes(
                  mclVa(k, "k"),
                  mlfProd(mclMinus(mclVa(fs, "fs"), mclVa(z, "z")), NULL)),
                mlfProd(mclMinus(mclVa(fs, "fs"), mclVa(p, "p")), NULL)));
            /*
             * zd = [zd;-ones(length(pd)-length(zd),1)];  % Add extra zeros at -1
             */
            mlfAssign(
              &zd,
              mlfVertcat(
                mclVv(zd, "zd"),
                mclUminus(
                  mlfOnes(
                    mlfScalar(
                      mclLengthInt(mclVv(*pd, "pd"))
                      - mclLengthInt(mclVv(zd, "zd"))),
                    _mxarray4_,
                    NULL)),
                NULL));
        /*
         * 
         * elseif (md == 1 & mn == 1) | nargout == 4 %
         */
        } else {
            mxDestroyArray(a_);
            {
                mxArray * a_0
                  = mclInitialize(mclEq(mclVv(md, "md"), _mxarray4_));
                if (mlfTobool(a_0)) {
                    mlfAssign(
                      &a_0, mclAnd(a_0, mclEq(mclVv(mn, "mn"), _mxarray4_)));
                } else {
                    mlfAssign(&a_0, mlfScalar(0));
                }
                if (mlfTobool(a_0)
                    || mlfTobool(mclOr(a_0, mclBoolToArray(nargout_ == 4)))) {
                    mxDestroyArray(a_0);
                    /*
                     * if nargout == 4		% State-space case
                     */
                    if (nargout_ == 4) {
                        /*
                         * a = z; b = p; c = k; d = fs; fs = fp;
                         */
                        mlfAssign(&a, mclVa(z, "z"));
                        mlfAssign(&b, mclVa(p, "p"));
                        mlfAssign(&c, mclVa(k, "k"));
                        mlfAssign(&d, mclVa(fs, "fs"));
                        mlfAssign(&fs, mclVa(fp, "fp"));
                        /*
                         * error(abcdchk(a,b,c,d));
                         */
                        mlfError(
                          mlfAbcdchk(
                            NULL,
                            NULL,
                            NULL,
                            NULL,
                            mclVv(a, "a"),
                            mclVv(b, "b"),
                            mclVv(c, "c"),
                            mclVv(d, "d")),
                          NULL);
                        /*
                         * if nargin == 6			% Prewarp
                         */
                        if (nargin_ == 6) {
                            /*
                             * fp = fp1;		% Decode arguments
                             */
                            mlfAssign(&fp, mclVa(fp1, "fp1"));
                            /*
                             * fp = 2*pi*fp;
                             */
                            mlfAssign(
                              &fp, mclMtimes(_mxarray3_, mclVa(fp, "fp")));
                            /*
                             * fs = fp/tan(fp/fs/2)/2;
                             */
                            mlfAssign(
                              &fs,
                              mclMrdivide(
                                mclMrdivide(
                                  mclVa(fp, "fp"),
                                  mlfTan(
                                    mclMrdivide(
                                      mclMrdivide(
                                        mclVa(fp, "fp"), mclVa(fs, "fs")),
                                      _mxarray0_))),
                                _mxarray0_));
                        /*
                         * end
                         */
                        }
                    /*
                     * else			% Transfer function case
                     */
                    } else {
                        /*
                         * if nn > nd
                         */
                        if (mclGtBool(mclVv(nn, "nn"), mclVv(nd, "nd"))) {
                            /*
                             * error('Numerator cannot be higher order than denominator.')
                             */
                            mlfError(_mxarray1_, NULL);
                        /*
                         * end
                         */
                        }
                        /*
                         * num = z; den = p;		% Decode arguments
                         */
                        mlfAssign(&num, mclVa(z, "z"));
                        mlfAssign(&den, mclVa(p, "p"));
                        /*
                         * if nargin == 4			% Prewarp
                         */
                        if (nargin_ == 4) {
                            /*
                             * fp = fs; fs = k;	% Decode arguments
                             */
                            mlfAssign(&fp, mclVa(fs, "fs"));
                            mlfAssign(&fs, mclVa(k, "k"));
                            /*
                             * fp = 2*pi*fp;
                             */
                            mlfAssign(
                              &fp, mclMtimes(_mxarray3_, mclVa(fp, "fp")));
                            /*
                             * fs = fp/tan(fp/fs/2)/2;
                             */
                            mlfAssign(
                              &fs,
                              mclMrdivide(
                                mclMrdivide(
                                  mclVa(fp, "fp"),
                                  mlfTan(
                                    mclMrdivide(
                                      mclMrdivide(
                                        mclVa(fp, "fp"), mclVa(fs, "fs")),
                                      _mxarray0_))),
                                _mxarray0_));
                        /*
                         * else
                         */
                        } else {
                            /*
                             * fs = k;			% Decode arguments
                             */
                            mlfAssign(&fs, mclVa(k, "k"));
                        /*
                         * end
                         */
                        }
                        /*
                         * % Put num(s)/den(s) in state-space canonical form.  
                         * [a,b,c,d] = tf2ss(num,den);
                         */
                        mlfAssign(
                          &a,
                          mlfTf2ss(
                            &b, &c, &d, mclVv(num, "num"), mclVv(den, "den")));
                    /*
                     * end
                     */
                    }
                    /*
                     * % Now do state-space version of bilinear transformation:
                     * t = 1/fs;
                     */
                    mlfAssign(&t, mclMrdivide(_mxarray4_, mclVa(fs, "fs")));
                    /*
                     * r = sqrt(t);
                     */
                    mlfAssign(&r, mlfSqrt(mclVv(t, "t")));
                    /*
                     * t1 = eye(size(a)) + a*t/2;
                     */
                    mlfAssign(
                      &t1,
                      mclPlus(
                        mlfEye(
                          mlfSize(mclValueVarargout(), mclVv(a, "a"), NULL),
                          NULL),
                        mclMrdivide(
                          mclMtimes(mclVv(a, "a"), mclVv(t, "t")),
                          _mxarray0_)));
                    /*
                     * t2 = eye(size(a)) - a*t/2;
                     */
                    mlfAssign(
                      &t2,
                      mclMinus(
                        mlfEye(
                          mlfSize(mclValueVarargout(), mclVv(a, "a"), NULL),
                          NULL),
                        mclMrdivide(
                          mclMtimes(mclVv(a, "a"), mclVv(t, "t")),
                          _mxarray0_)));
                    /*
                     * ad = t2\t1;
                     */
                    mlfAssign(
                      &ad, mlfMldivide(mclVv(t2, "t2"), mclVv(t1, "t1")));
                    /*
                     * bd = t/r*(t2\b);
                     */
                    mlfAssign(
                      &bd,
                      mclMtimes(
                        mclMrdivide(mclVv(t, "t"), mclVv(r, "r")),
                        mlfMldivide(mclVv(t2, "t2"), mclVv(b, "b"))));
                    /*
                     * cd = r*c/t2;
                     */
                    mlfAssign(
                      &cd,
                      mclMrdivide(
                        mclMtimes(mclVv(r, "r"), mclVv(c, "c")),
                        mclVv(t2, "t2")));
                    /*
                     * dd = c/t2*b*t/2 + d;
                     */
                    mlfAssign(
                      dd,
                      mclPlus(
                        mclMrdivide(
                          mclMtimes(
                            mclMtimes(
                              mclMrdivide(mclVv(c, "c"), mclVv(t2, "t2")),
                              mclVv(b, "b")),
                            mclVv(t, "t")),
                          _mxarray0_),
                        mclVv(d, "d")));
                    /*
                     * if nargout == 4
                     */
                    if (nargout_ == 4) {
                        /*
                         * zd = ad; pd = bd; kd = cd;
                         */
                        mlfAssign(&zd, mclVv(ad, "ad"));
                        mlfAssign(pd, mclVv(bd, "bd"));
                        mlfAssign(kd, mclVv(cd, "cd"));
                    /*
                     * else
                     */
                    } else {
                        /*
                         * % Convert back to transfer function form:
                         * p = poly(ad);
                         */
                        mlfAssign(&p, mlfPoly(mclVv(ad, "ad")));
                        /*
                         * zd = poly(ad-bd*cd)+(dd-1)*p;
                         */
                        mlfAssign(
                          &zd,
                          mclPlus(
                            mlfPoly(
                              mclMinus(
                                mclVv(ad, "ad"),
                                mclMtimes(mclVv(bd, "bd"), mclVv(cd, "cd")))),
                            mclMtimes(
                              mclMinus(mclVv(*dd, "dd"), _mxarray4_),
                              mclVa(p, "p"))));
                        /*
                         * pd = p;
                         */
                        mlfAssign(pd, mclVa(p, "p"));
                    /*
                     * end
                     */
                    }
                /*
                 * else
                 */
                } else {
                    mxDestroyArray(a_0);
                    /*
                     * error('First two arguments must have the same orientation.')
                     */
                    mlfError(_mxarray5_, NULL);
                }
            }
        }
    /*
     * end
     */
    }
    mclValidateOutput(zd, 1, nargout_, "zd", "bilinear");
    mclValidateOutput(*pd, 2, nargout_, "pd", "bilinear");
    mclValidateOutput(*kd, 3, nargout_, "kd", "bilinear");
    mclValidateOutput(*dd, 4, nargout_, "dd", "bilinear");
    mxDestroyArray(mn);
    mxDestroyArray(nn);
    mxDestroyArray(md);
    mxDestroyArray(nd);
    mxDestroyArray(ans);
    mxDestroyArray(a);
    mxDestroyArray(b);
    mxDestroyArray(c);
    mxDestroyArray(d);
    mxDestroyArray(num);
    mxDestroyArray(den);
    mxDestroyArray(t);
    mxDestroyArray(r);
    mxDestroyArray(t1);
    mxDestroyArray(t2);
    mxDestroyArray(ad);
    mxDestroyArray(bd);
    mxDestroyArray(cd);
    mxDestroyArray(fp1);
    mxDestroyArray(fp);
    mxDestroyArray(fs);
    mxDestroyArray(k);
    mxDestroyArray(p);
    mxDestroyArray(z);
    mclSetCurrentLocalFunctionTable(save_local_function_table_);
    return zd;
}
