/* 
 *   COMPLEX Construct complex data from real and imaginary parts.
 *   COMPLEX(A,B) returns the complex result A + Bi, where A and B
 *   are identically sized real arrays or scalars of the same data
 *   type.
 *
 *   COMPLEX(A) returns the complex result A + 0i, where A must
 *   be real.
 *
 *   Note that the expression A+i*B or A+j*B will give identical
 *   results if A and B are double-precision and i or j has not been
 *   assigned.  Use COMPLEX if ambiguity in the variables "i" or "j"
 *   might arise, or if A and B are not double-precision.
 *
 *  Author(s): R. Firtion, D. Orofino
 *  Copyright 1984-2000 The MathWorks, Inc. 
 *  $Revision: 1.3 $  $Date: 2000/06/01 00:40:15 $
 */
#include "mex.h"
#include "matrix.h"

enum {INPUT_RE_ARGC=0, INPUT_IM_ARGC, MAX_NUM_INPUTS};
enum {OUTPUT_ARGC=0};

#define INPUT_RE  prhs[INPUT_RE_ARGC]
#define INPUT_IM  prhs[INPUT_IM_ARGC]
#define OUTPUT    plhs[OUTPUT_ARGC]

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
    if (nrhs < 1) {
        mexErrMsgTxt("Not enough input arguments.");
    } else if (nrhs > MAX_NUM_INPUTS) {
        mexErrMsgTxt("Too many input arguments.");
    }

    {
        const int  Re_NumDims  = mxGetNumberOfDimensions(INPUT_RE);
        const int  Re_NumEle   = mxGetNumberOfElements(INPUT_RE);
        const int *pRe_Dims    = mxGetDimensions(INPUT_RE);
        char      *A           = (char *)mxGetPr(INPUT_RE);
        const int  BytesPerEle = mxGetElementSize(INPUT_RE);

        /* Get imag input, if present.
         * Also, guarantee that Im_NumEle != Re_NumEle if rhs==1
         */
        const int  Im_NumDims  = (nrhs>1) ? mxGetNumberOfDimensions(INPUT_IM) : -1;
        int        Im_NumEle   = (nrhs>1) ? mxGetNumberOfElements(INPUT_IM)   : -1;
        const int *pIm_Dims    = (nrhs>1) ? mxGetDimensions(INPUT_IM)         :  NULL;
        char      *B           = (nrhs>1) ? (char *)mxGetPr(INPUT_IM)         :  NULL;

        char	 *Cr, *Ci;
        int       i;

        /* Check data type is valid: */
        if(!mxIsNumeric(INPUT_RE) || mxIsSparse(INPUT_RE) || mxIsComplex(INPUT_RE)) {
            mexErrMsgTxt("Real input A must be numeric, real, and full.");
        }

        if (nrhs>1) {
            /* Only check B is two RHS arguments are passed: */

            if(!mxIsNumeric(INPUT_IM) || mxIsSparse(INPUT_IM) || mxIsComplex(INPUT_IM)) {
                mexErrMsgTxt("Imaginary input B must be numeric, real, and full.");
            }

            /* Check data types are the same: */
            if(mxGetClassID(INPUT_RE) != mxGetClassID(INPUT_IM)){
                mexErrMsgTxt("Input arrays must be the same data type.");
            }

            /* Check for conforming array sizes: */
            if ((Re_NumEle != Im_NumEle) && (Re_NumEle != 1) && (Im_NumEle != 1)) {
                mexErrMsgTxt("Input arrays must have equal number of dimensions.");
            }
        }

        if (Re_NumEle != Im_NumEle) {
            /* Scalar expansion: */

            /* Force imag part to expand if imag part isn't passed: */
            const int     ExpandRe  = (Re_NumEle == 1) && (nrhs > 1);
            const int     Re_Incr   = (ExpandRe) ? 0           : BytesPerEle;
            const int     Im_Incr   = (ExpandRe) ? BytesPerEle : 0;
            const int     out_nEle  = (ExpandRe) ? Im_NumEle   : Re_NumEle;
            const int     out_nDims = (ExpandRe) ? Im_NumDims  : Re_NumDims;
            const int    *out_pDims = (ExpandRe) ? pIm_Dims    : pRe_Dims;
            char         *Bzero     = NULL;

            OUTPUT = mxCreateNumericArray(out_nDims, out_pDims, mxGetClassID(INPUT_RE), mxCOMPLEX);
            Cr = (char *) mxGetPr(OUTPUT);
            Ci = (char *) mxGetPi(OUTPUT);

            if (nrhs == 1) {
                /* Allocate a zero for imag part: */
                B = Bzero = (char *)mxCalloc(BytesPerEle, sizeof(char));
                if (B == NULL) {
                    mexErrMsgTxt("Memory allocation failure.");
                }
            }
        
            for(i=0; i < out_nEle; i++) {
                memcpy(Cr, A, BytesPerEle);
                memcpy(Ci, B, BytesPerEle);
                A  += Re_Incr;
                B  += Im_Incr;
                Cr += BytesPerEle;
                Ci += BytesPerEle;
            }

            if (nrhs == 1) {
                mxFree(Bzero);
            }

        } else {
            /* No scalar expansion: */

            /* Check individual sizes of each dimension */
            for(i=0; i<Re_NumDims; i++){
                if(pRe_Dims[i] != pIm_Dims[i]) {
                    mexErrMsgTxt("Input arrays must have identical dimensions.");
                }
            }

            OUTPUT = mxCreateNumericArray(Re_NumDims, pRe_Dims, mxGetClassID(INPUT_RE), mxCOMPLEX);
            Cr = (char *) mxGetPr(OUTPUT);
            Ci = (char *) mxGetPi(OUTPUT);

            memcpy(Cr, A, Re_NumEle * BytesPerEle);
            memcpy(Ci, B, Im_NumEle * BytesPerEle);
        }
    }
}

/* [EOF] complex.c */