#include "mex.h"
#include "math.h"
#include "matrix.h"
#include "sparse.h"
#include "myheader.h"
#include "allomorphy.h"

void mexFunction( int nlhs, mxArray *plhs[],
                  int nrhs, const mxArray *prhs[] )
{
  TENSOR2 u,v,psi,rho;
  TENSOR3 pmp;
  int ell;

  makeT2ad(u,           IN(1));
  makeT2ad(v,           IN(2));
  ell =  (int) *mxGetPr(IN(3));
  makeT3ad(pmp,         IN(4));

  /* allocate auxilliary variables */
  /* tmpar = allocTENSOR2(&psi,d2(u),d2(v)); /* must free later */

  OU(1) = allocTENSOR2(&psi,d2(u)+1,d2(v)+1);

  jointstrprob(psi,u,v,ell-1,pmp);
}