/*
 * MATLAB Compiler: 3.0
 * Date: Mon Sep 20 09:40:11 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 "thetah.h"
#include "libmatlbm.h"
static mxArray * _mxarray0_;

void InitializeModule_thetah(void) {
    _mxarray0_ = mclInitializeDouble(1.0);
}

void TerminateModule_thetah(void) {
    mxDestroyArray(_mxarray0_);
}

static mxArray * Mthetah(int nargout_,
                         mxArray * h,
                         mxArray * a,
                         mxArray * nn,
                         mxArray * nsat,
                         mxArray * nres,
                         mxArray * ss,
                         mxArray * hh0);

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

/*
 * The function "mlfThetah" contains the normal interface for the "thetah"
 * M-function from file "u:\final_recarga6\recarga\thetah.m" (lines 1-21). This
 * function processes any input arguments and passes them to the implementation
 * version of the function, appearing above.
 */
mxArray * mlfThetah(mxArray * h,
                    mxArray * a,
                    mxArray * nn,
                    mxArray * nsat,
                    mxArray * nres,
                    mxArray * ss,
                    mxArray * hh0) {
    int nargout = 1;
    mxArray * theta = NULL;
    mlfEnterNewContext(0, 7, h, a, nn, nsat, nres, ss, hh0);
    theta = Mthetah(nargout, h, a, nn, nsat, nres, ss, hh0);
    mlfRestorePreviousContext(0, 7, h, a, nn, nsat, nres, ss, hh0);
    return mlfReturnValue(theta);
}

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

/*
 * The function "Mthetah" is the implementation version of the "thetah"
 * M-function from file "u:\final_recarga6\recarga\thetah.m" (lines 1-21). 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 theta=THETAH(h,a,nn,nsat,nres,ss,hh0);
 */
static mxArray * Mthetah(int nargout_,
                         mxArray * h,
                         mxArray * a,
                         mxArray * nn,
                         mxArray * nsat,
                         mxArray * nres,
                         mxArray * ss,
                         mxArray * hh0) {
    mexLocalFunctionTable save_local_function_table_
      = mclSetCurrentLocalFunctionTable(&_local_function_table_thetah);
    mxArray * theta = NULL;
    mxArray * b0 = NULL;
    mxArray * b = NULL;
    mxArray * mm = NULL;
    mxArray * x = NULL;
    mclCopyArray(&h);
    mclCopyArray(&a);
    mclCopyArray(&nn);
    mclCopyArray(&nsat);
    mclCopyArray(&nres);
    mclCopyArray(&ss);
    mclCopyArray(&hh0);
    /*
     * %
     * % thetah.M - Calculates the volumetric moisture content as a function
     * %    of pressure head. 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: h (pressure head), a (van Genuchten parameter alpha),
     * %		nn (vG param. n), nsat (sat. water content), nres (residual)
     * % OUTPUT: soil water content
     * %
     * %p=abs(h);					% use absolute value of h
     * %p0=abs(hh0);
     * %p=p-1/a/2;    % correct for air-entry suction=1/a
     * %if p<0; p=0; end;	% check for positive pressure heads
     * % thus, will return sat. moisture content for those saturated nodes
     * x=h>hh0;    % ones where hi>hh0
     */
    mlfAssign(&x, mclGt(mclVa(h, "h"), mclVa(hh0, "hh0")));
    /*
     * %p(k)=zeros(size(k));    % puts h=0 where h>0 saturated
     * mm=1-(1/nn);
     */
    mlfAssign(
      &mm, mclMinus(_mxarray0_, mclMrdivide(_mxarray0_, mclVa(nn, "nn"))));
    /*
     * b=(1+(-a.*h).^nn).^mm;
     */
    mlfAssign(
      &b,
      mlfPower(
        mclPlus(
          _mxarray0_,
          mlfPower(
            mclTimes(mclUminus(mclVa(a, "a")), mclVa(h, "h")),
            mclVa(nn, "nn"))),
        mclVv(mm, "mm")));
    /*
     * b0=(1+(-a.*hh0).^nn).^mm;
     */
    mlfAssign(
      &b0,
      mlfPower(
        mclPlus(
          _mxarray0_,
          mlfPower(
            mclTimes(mclUminus(mclVa(a, "a")), mclVa(hh0, "hh0")),
            mclVa(nn, "nn"))),
        mclVv(mm, "mm")));
    /*
     * theta=(1-x).*(nres+(nsat-nres)./b)+x.*(nres+(nsat-nres)./b0+ss.*(h-hh0));
     */
    mlfAssign(
      &theta,
      mclPlus(
        mclTimes(
          mclMinus(_mxarray0_, mclVv(x, "x")),
          mclPlus(
            mclVa(nres, "nres"),
            mclRdivide(
              mclMinus(mclVa(nsat, "nsat"), mclVa(nres, "nres")),
              mclVv(b, "b")))),
        mclTimes(
          mclVv(x, "x"),
          mclPlus(
            mclPlus(
              mclVa(nres, "nres"),
              mclRdivide(
                mclMinus(mclVa(nsat, "nsat"), mclVa(nres, "nres")),
                mclVv(b0, "b0"))),
            mclTimes(
              mclVa(ss, "ss"), mclMinus(mclVa(h, "h"), mclVa(hh0, "hh0")))))));
    mclValidateOutput(theta, 1, nargout_, "theta", "thetah");
    mxDestroyArray(x);
    mxDestroyArray(mm);
    mxDestroyArray(b);
    mxDestroyArray(b0);
    mxDestroyArray(hh0);
    mxDestroyArray(ss);
    mxDestroyArray(nres);
    mxDestroyArray(nsat);
    mxDestroyArray(nn);
    mxDestroyArray(a);
    mxDestroyArray(h);
    mclSetCurrentLocalFunctionTable(save_local_function_table_);
    return theta;
}
