/*
 * MATLAB Compiler: 3.0
 * Date: Mon Sep 20 09:40:13 2004
 * 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" "-B" "sgl" "-m" "-W"
 * "main" "-L" "C" "-t" "-T" "link:exe" "-h" "libmmfile.mlib" "-W" "mainhg"
 * "libmwsglm.mlib" "recarga_multiple" 
 */
#include "relcond.h"
#include "libmatlbm.h"
static mxArray * _mxarray0_;
static mxArray * _mxarray1_;
static mxArray * _mxarray2_;
static mxArray * _mxarray3_;
static mxArray * _mxarray4_;
static mxArray * _mxarray5_;

void InitializeModule_relcond(void) {
    _mxarray0_ = mclInitializeDouble(1.0);
    _mxarray1_ = mclInitializeDouble(0.0);
    _mxarray2_ = mclInitializeDouble(.01);
    _mxarray3_ = mclInitializeDouble(.5);
    _mxarray4_ = mclInitializeDouble(-1.0);
    _mxarray5_ = mclInitializeDouble(2.0);
}

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

static mxArray * Mrelcond(int nargout_,
                          mxArray * th,
                          mxArray * a,
                          mxArray * nn,
                          mxArray * thr,
                          mxArray * ths);

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

/*
 * The function "mlfRelcond" contains the normal interface for the "relcond"
 * M-function from file "u:\final_recarga6\recarga\relcond.m" (lines 1-20).
 * This function processes any input arguments and passes them to the
 * implementation version of the function, appearing above.
 */
mxArray * mlfRelcond(mxArray * th,
                     mxArray * a,
                     mxArray * nn,
                     mxArray * thr,
                     mxArray * ths) {
    int nargout = 1;
    mxArray * krelative = NULL;
    mlfEnterNewContext(0, 5, th, a, nn, thr, ths);
    krelative = Mrelcond(nargout, th, a, nn, thr, ths);
    mlfRestorePreviousContext(0, 5, th, a, nn, thr, ths);
    return mlfReturnValue(krelative);
}

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

/*
 * The function "Mrelcond" is the implementation version of the "relcond"
 * M-function from file "u:\final_recarga6\recarga\relcond.m" (lines 1-20). 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 krelative=relcond(th,a,nn,thr,ths);
 */
static mxArray * Mrelcond(int nargout_,
                          mxArray * th,
                          mxArray * a,
                          mxArray * nn,
                          mxArray * thr,
                          mxArray * ths) {
    mexLocalFunctionTable save_local_function_table_
      = mclSetCurrentLocalFunctionTable(&_local_function_table_relcond);
    mxArray * krelative = NULL;
    mxArray * mm = NULL;
    mxArray * j = NULL;
    mxArray * k = NULL;
    mxArray * thn = NULL;
    mclCopyArray(&th);
    mclCopyArray(&a);
    mclCopyArray(&nn);
    mclCopyArray(&thr);
    mclCopyArray(&ths);
    /*
     * %
     * % RELCOND.M - Calculates the relative hydraulic conductivity as a
     * %    function of soil water content. Reference: van Genuchten, M. T.
     * %    1980. A closed-form equation for predicting the hydraulic
     * %    conductivity of unsaturated soils. Soil Sci. Soc. Amer.
     * %    J. 44, 892-898.
     * % INPUT: th (volumetric soil water content), a (van Genuchten paramater alpha),
     * %		nn (vG param. n), thr (residual water content), ths (sat. water content)
     * % OUTPUT: relative hydraulic conductivity [0,1]
     * %
     * thn=(th-thr)/(ths-thr);					% normalized water content
     */
    mlfAssign(
      &thn,
      mclMrdivide(
        mclMinus(mclVa(th, "th"), mclVa(thr, "thr")),
        mclMinus(mclVa(ths, "ths"), mclVa(thr, "thr"))));
    /*
     * if thn>1;thn=1;end; % check for oversaturation
     */
    if (mclGtBool(mclVv(thn, "thn"), _mxarray0_)) {
        mlfAssign(&thn, _mxarray0_);
    }
    /*
     * if thn<0;thn=0;end; % check for undervalues
     */
    if (mclLtBool(mclVv(thn, "thn"), _mxarray1_)) {
        mlfAssign(&thn, _mxarray1_);
    }
    /*
     * k=find(thn<=0); thn(k)=zeros(size(k))+0.01; %check for undervalues
     */
    mlfAssign(&k, mlfFind(NULL, NULL, mclLe(mclVv(thn, "thn"), _mxarray1_)));
    mclArrayAssign1(
      &thn,
      mclPlus(
        mlfZeros(mlfSize(mclValueVarargout(), mclVv(k, "k"), NULL), NULL),
        _mxarray2_),
      mclVv(k, "k"));
    /*
     * j=find(thn>=1); thn(j)=ones(size(j)); %check for oversaturation
     */
    mlfAssign(&j, mlfFind(NULL, NULL, mclGe(mclVv(thn, "thn"), _mxarray0_)));
    mclArrayAssign1(
      &thn,
      mlfOnes(mlfSize(mclValueVarargout(), mclVv(j, "j"), NULL), NULL),
      mclVv(j, "j"));
    /*
     * mm=1-(1/nn);
     */
    mlfAssign(
      &mm, mclMinus(_mxarray0_, mclMrdivide(_mxarray0_, mclVa(nn, "nn"))));
    /*
     * krelative=thn.^(0.5).*(1-thn.*(thn.^(-1/mm)-1).^mm).^2;
     */
    mlfAssign(
      &krelative,
      mclTimes(
        mlfPower(mclVv(thn, "thn"), _mxarray3_),
        mlfPower(
          mclMinus(
            _mxarray0_,
            mclTimes(
              mclVv(thn, "thn"),
              mlfPower(
                mclMinus(
                  mlfPower(
                    mclVv(thn, "thn"),
                    mclMrdivide(_mxarray4_, mclVv(mm, "mm"))),
                  _mxarray0_),
                mclVv(mm, "mm")))),
          _mxarray5_)));
    mclValidateOutput(krelative, 1, nargout_, "krelative", "relcond");
    mxDestroyArray(thn);
    mxDestroyArray(k);
    mxDestroyArray(j);
    mxDestroyArray(mm);
    mxDestroyArray(ths);
    mxDestroyArray(thr);
    mxDestroyArray(nn);
    mxDestroyArray(a);
    mxDestroyArray(th);
    mclSetCurrentLocalFunctionTable(save_local_function_table_);
    return krelative;
    /*
     * %if thn==0; krelative=0; end;
     */
}
