File mttroot/mtt/lib/rep/sfun_rep/mex_ae.c.tmpl artifact ca0d3e6d3b part of check-in 9f2a0a7535

/* -*-c-*-	Put emacs into c-mode
 * <mtt_model_name>_sfun_ae.c:
 * Matlab mex algebraic equations for <mtt_model_name>
 */

#if <have_algebraic_equations> /* have_algebraic_equations */

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <mex.h>
#include "sfun_debug.h"
#include "useful-functions.hh"
#include "<mtt_model_name>_def.h"
#include "<mtt_model_name>_sympar.h"
#include "<mtt_model_name>_cr.h"

/* utility procedures */

double *
array_of_double (size_t n)
{
  void *p = calloc (n, sizeof (double));
  if (! p) {
    fprintf (stderr, "*** Error: failed to allocate memory\n");
  }
  return (double *) p;
}

/* system equations */

static void
<mtt_model_name>_ae (double *mttyz,
		     const double *mttx,
		     const double *mttu,
		     const double mttt,
		     const double *mttpar)
{
#include "<mtt_model_name>_ae.c"
  PRINT_LEAVE;
}

/* generic mex function */

#ifdef __cplusplus
extern "C" {
#endif

void
mexFunction (int nlhs, mxArray *plhs[],
	     int nrhs, const mxArray *prhs[])
{
  double *mttyz;		/* residuals */

  double *mttx;			/* states */
  double *mttu;			/* known + unknown inputs */
  double mttt;			/* time */
  double *mttpar;		/* parameters */

  unsigned int i;
  double *p;

  PRINT_ENTER;

  mttyz		= array_of_double (MTTNYZ);

  mttx		= array_of_double (MTTNX);
  mttu		= array_of_double (MTTNU + MTTNYZ);
  mttpar	= array_of_double (MTTNPAR);

  /* get trial values */
  p = mxGetPr (prhs[0]);
  for (i = 0; i < MTTNYZ; i++) {
    mttu[MTTNU + i] = p[i];
  }

  /* get states */
  p = mxGetPr (prhs[1]);
  for (i = 0; i < MTTNX; i++) {
    mttx[i] = p[i];
  }

  /* get known inputs */
  p = mxGetPr (prhs[2]);
  for (i = 0; i < MTTNU; i++) {
    mttu[i] = p[i];
  }

  /* get time */
  p = mxGetPr (prhs[3]);
  mttt = *p;

  /* get parameters */
  p = mxGetPr (prhs[4]);
  for (i = 0; i < MTTNPAR; i++) {
    mttpar[i] = p[i];
  }
  
  /* evaluate residuals */
  <mtt_model_name>_ae (mttyz, mttx, mttu, mttt, mttpar);
  
  /* return residuals */
  plhs[0] = mxCreateDoubleMatrix (MTTNYZ, 1, mxREAL);
  p = mxGetPr (plhs[0]);
  for (i = 0; i < MTTNYZ; i++) {
    p[i] = mttyz[i];
  }

  /* release memory */
  free (mttx);
  free (mttu);
  free (mttpar);
  free (mttyz);

  PRINT_LEAVE;
}

#ifdef __cplusplus
}
#endif

#endif /* have_algebraic_equations */