Index: ThirdParty/levmar/matlab/CMakeLists.txt
===================================================================
--- ThirdParty/levmar/matlab/CMakeLists.txt	(revision 8ce1a95e94bdc842511db41cb292238a19558df0)
+++ ThirdParty/levmar/matlab/CMakeLists.txt	(revision 8ce1a95e94bdc842511db41cb292238a19558df0)
@@ -0,0 +1,58 @@
+# CMake file for levmar's MEX-file; see http://www.cmake.org
+# Requires FindMatlab.cmake included with cmake
+
+PROJECT(LEVMARMEX)
+#CMAKE_MINIMUM_REQUIRED(VERSION 1.4)
+
+INCLUDE("C:/Program Files/CMake 2.4/share/cmake-2.4/Modules/FindMatlab.cmake")
+
+# f2c is sometimes equivalent to libF77 & libI77; in that case, set HAVE_F2C to 0
+SET(HAVE_F2C 1 CACHE BOOL "Do we have f2c or F77/I77?" )
+
+# the directory where the lapack/blas/f2c libraries reside
+SET(LAPACKBLAS_DIR /usr/lib CACHE PATH "Path to lapack/blas libraries")
+
+# the directory where levmar.h resides
+SET(LM_H_DIR .. CACHE PATH "Path to levmar.h")
+# the directory where the levmar library resides
+SET(LEVMAR_DIR .. CACHE PATH "Path to levmar library")
+
+# actual names for the lapack/blas/f2c libraries
+SET(LAPACK_LIB lapack CACHE STRING "The name of the lapack library")
+SET(BLAS_LIB blas CACHE STRING "The name of the blas library")
+IF(HAVE_F2C)
+  SET(F2C_LIB f2c CACHE STRING "The name of the f2c library")
+ELSE(HAVE_F2C)
+  SET(F77_LIB libF77 CACHE STRING "The name of the F77 library")
+  SET(I77_LIB libI77 CACHE STRING "The name of the I77 library")
+ENDIF(HAVE_F2C)
+
+########################## NO CHANGES BEYOND THIS POINT ##########################
+
+INCLUDE_DIRECTORIES(${LM_H_DIR})
+LINK_DIRECTORIES(${LAPACKBLAS_DIR} ${LEVMAR_DIR})
+
+SET(SRC levmar.c)
+
+# naming conventions for the generated file's suffix
+IF(WIN32)
+  SET(SUFFIX ".mexw32")
+ELSE(WIN32)
+  SET(SUFFIX ".mexglx")
+ENDIF(WIN32)
+
+SET(OUTNAME "levmar${SUFFIX}")
+
+ADD_LIBRARY(${OUTNAME} MODULE ${SRC})
+
+IF(HAVE_F2C)
+	ADD_CUSTOM_COMMAND(OUTPUT ${OUTNAME}
+                   DEPENDS ${SRC}
+                   COMMAND mex
+                   ARGS -O -I${LM_H_DIR} ${SRC} -I${MATLAB_INCLUDE_DIR} -L${LAPACKBLAS_DIR} -L${LEVMAR_DIR} -L${MATLAB_MEX_LIBRARY} -llevmar -l${LAPACK_LIB} -l${BLAS_LIB} -l${F2C_LIB} -output ${MATLAB_LIBRARIES} ${OUTNAME})
+ELSE(HAVE_F2C)
+	ADD_CUSTOM_COMMAND(OUTPUT ${OUTNAME}
+                   DEPENDS ${SRC}
+                   COMMAND mex
+                   ARGS -O -I${LM_H_DIR} ${SRC} -I${MATLAB_INCLUDE_DIR} -L${LAPACKBLAS_DIR} -L${LEVMAR_DIR} -L${MATLAB_MEX_LIBRARY} -llevmar -l${LAPACK_LIB} -l${BLAS_LIB} -l${F77_LIB} -l${I77_LIB} ${MATLAB_LIBRARIES} -output ${OUTNAME})
+ENDIF(HAVE_F2C)
Index: ThirdParty/levmar/matlab/Makefile
===================================================================
--- ThirdParty/levmar/matlab/Makefile	(revision 8ce1a95e94bdc842511db41cb292238a19558df0)
+++ ThirdParty/levmar/matlab/Makefile	(revision 8ce1a95e94bdc842511db41cb292238a19558df0)
@@ -0,0 +1,30 @@
+#
+# Unix/Linux Makefile for MATLAB interface to levmar
+#
+
+MEX=mex
+MEXCFLAGS=-I.. -O #-g
+# WHEN USING LAPACK, CHANGE THE NEXT TWO LINES TO WHERE YOUR COMPILED LAPACK/BLAS & F2C LIBS ARE!
+LAPACKBLASLIBS_PATH=/usr/lib
+F2CLIBS_PATH=/usr/local/lib
+
+
+# I had to specify the absolute path to the libs, otherwise mex linked against their dynamic versions...
+INTFACESRCS=levmar.c
+LAPACKLIBS=$(LAPACKBLASLIBS_PATH)/liblapack.a $(LAPACKBLASLIBS_PATH)/libblas.a $(F2CLIBS_PATH)/libf2c.a
+                                 # On systems with a FORTRAN (not f2c'ed) version of LAPACK, libf2c.a is
+                                 # not necessary; on others, libf2c.a comes in two parts: libF77.a and libI77.a
+
+LIBS=$(LAPACKLIBS)
+
+dummy: $(INTFACESRCS)
+	$(MEX) $(MEXCFLAGS) $(INTFACESRCS) ../liblevmar.a $(LIBS)
+
+clean:
+	@rm -f levmar.mexglx
+
+depend:
+	makedepend -f Makefile $(INTFACESRCS)
+
+# DO NOT DELETE THIS LINE -- make depend depends on it.
+
Index: ThirdParty/levmar/matlab/Makefile.w32
===================================================================
--- ThirdParty/levmar/matlab/Makefile.w32	(revision 8ce1a95e94bdc842511db41cb292238a19558df0)
+++ ThirdParty/levmar/matlab/Makefile.w32	(revision 8ce1a95e94bdc842511db41cb292238a19558df0)
@@ -0,0 +1,26 @@
+#
+# Windows Makefile for MATLAB interface to levmar
+#
+
+MEX=mex
+MEXCFLAGS=-I.. -O #-g
+# WHEN USING LAPACK, CHANGE THE NEXT TWO LINES TO WHERE YOUR COMPILED LAPACK/BLAS & F2C LIBS ARE!
+LAPACKBLASLIBS_PATH=C:\src\lib
+F2CLIBS_PATH=$(LAPACKBLASLIBS_PATH) # define appropriately if not identical to LAPACKBLASLIBS_PATH
+
+
+INTFACESRCS=levmar.c
+LAPACKLIBS=$(LAPACKBLASLIBS_PATH)/clapack.lib $(LAPACKBLASLIBS_PATH)/blas.lib $(F2CLIBS_PATH)/libF77.lib $(F2CLIBS_PATH)/libI77.lib
+LIBS=$(LAPACKLIBS)
+
+dummy: $(INTFACESRCS)
+	$(MEX) $(MEXCFLAGS) $(INTFACESRCS) ../levmar.lib $(LIBS)
+
+clean:
+	-del levmar.mexw32
+
+depend:
+	makedepend -f Makefile $(INTFACESRCS)
+
+# DO NOT DELETE THIS LINE -- make depend depends on it.
+
Index: ThirdParty/levmar/matlab/README.txt
===================================================================
--- ThirdParty/levmar/matlab/README.txt	(revision 8ce1a95e94bdc842511db41cb292238a19558df0)
+++ ThirdParty/levmar/matlab/README.txt	(revision 8ce1a95e94bdc842511db41cb292238a19558df0)
@@ -0,0 +1,35 @@
+This directory contains a matlab MEX interface to levmar. This interface
+has been tested with Matlab v. 6.5 R13 under linux and v. 7.4 R2007 under Windows.
+Users have also reported success with GNU Octave.
+
+FILES
+The following files are included:
+levmar.c: C MEX-file for levmar
+Makefile: UNIX makefile for compiling levmar.c using mex
+Makefile.w32: Windows makefile for compiling levmar.c using mex
+levmar.m: Documentation for the MEX interface
+lmdemo.m: Demonstration of using the MEX interface; run as matlab < lmdemo.m
+
+*.m: Matlab functions implementing various objective functions and their Jacobians.
+     For instance, meyer.m implements the objective function for Meyer's (reformulated)
+     problem and jacmeyer.m implements its Jacobian.
+
+
+
+COMPILING
+Use the provided Makefile or Makefile.w32, depending on your platform.
+Alternatively, levmar.c can be compiled from matlab's prompt with a
+command like
+
+mex -DHAVE_LAPACK -I.. -O -L<levmar library dir> -L<blas/lapack libraries dir> levmar.c -llevmar -lclapack -lblas -lf2c
+          
+Make sure that you substitute the angle brackets with the correct paths to
+the levmar and the blas/lapack directories. Also, on certain systems,
+-lf2c should be changed to -llibF77 -llibI77
+If your mex compiler has not been configured, the following command should be run first:
+
+mex -setup 
+
+
+TESTING
+After compiling, execute lmdemo.m with matlab < lmdemo.m 
Index: ThirdParty/levmar/matlab/bt3.m
===================================================================
--- ThirdParty/levmar/matlab/bt3.m	(revision 8ce1a95e94bdc842511db41cb292238a19558df0)
+++ ThirdParty/levmar/matlab/bt3.m	(revision 8ce1a95e94bdc842511db41cb292238a19558df0)
@@ -0,0 +1,11 @@
+function x = bt3(p, adata)
+  n=5;
+
+  t1=p(1)-p(2);
+  t2=p(2)+p(3)-2.0;
+  t3=p(4)-1.0;
+  t4=p(5)-1.0;
+
+  for i=1:n
+    x(i)=t1*t1 + t2*t2 + t3*t3 + t4*t4;
+  end
Index: ThirdParty/levmar/matlab/expfit.m
===================================================================
--- ThirdParty/levmar/matlab/expfit.m	(revision 8ce1a95e94bdc842511db41cb292238a19558df0)
+++ ThirdParty/levmar/matlab/expfit.m	(revision 8ce1a95e94bdc842511db41cb292238a19558df0)
@@ -0,0 +1,8 @@
+function x = expfit(p, data)
+  n=data;
+
+% data1, data2 are actually unused
+
+  for i=1:n
+    x(i)=p(1)*exp(-p(2)*i) + p(3);
+  end
Index: ThirdParty/levmar/matlab/hs01.m
===================================================================
--- ThirdParty/levmar/matlab/hs01.m	(revision 8ce1a95e94bdc842511db41cb292238a19558df0)
+++ ThirdParty/levmar/matlab/hs01.m	(revision 8ce1a95e94bdc842511db41cb292238a19558df0)
@@ -0,0 +1,6 @@
+function x = hs01(p)
+  n=2;
+
+  t=p(1)*p(1);
+  x(1)=10.0*(p(2)-t);
+  x(2)=1.0-p(1);
Index: ThirdParty/levmar/matlab/jacbt3.m
===================================================================
--- ThirdParty/levmar/matlab/jacbt3.m	(revision 8ce1a95e94bdc842511db41cb292238a19558df0)
+++ ThirdParty/levmar/matlab/jacbt3.m	(revision 8ce1a95e94bdc842511db41cb292238a19558df0)
@@ -0,0 +1,13 @@
+function jac = jacbt3(p, adata)
+  n=5;
+  m=5;
+
+  t1=p(1)-p(2);
+  t2=p(2)+p(3)-2.0;
+  t3=p(4)-1.0;
+  t4=p(5)-1.0;
+
+  for i=1:n
+    jac(i, 1:m)=[2.0*t1, 2.0*(t2-t1), 2.0*t2, 2.0*t3, 2.0*t4];
+  end
+
Index: ThirdParty/levmar/matlab/jacexpfit.m
===================================================================
--- ThirdParty/levmar/matlab/jacexpfit.m	(revision 8ce1a95e94bdc842511db41cb292238a19558df0)
+++ ThirdParty/levmar/matlab/jacexpfit.m	(revision 8ce1a95e94bdc842511db41cb292238a19558df0)
@@ -0,0 +1,7 @@
+function jac = jacexpfit(p, data)
+  n=data;
+  m=max(size(p));
+
+  for i=1:n
+    jac(i, 1:m)=[exp(-p(2)*i), -p(1)*i*exp(-p(2)*i), 1.0];
+  end
Index: ThirdParty/levmar/matlab/jachs01.m
===================================================================
--- ThirdParty/levmar/matlab/jachs01.m	(revision 8ce1a95e94bdc842511db41cb292238a19558df0)
+++ ThirdParty/levmar/matlab/jachs01.m	(revision 8ce1a95e94bdc842511db41cb292238a19558df0)
@@ -0,0 +1,5 @@
+function jac = jachs01(p)
+  m=2;
+
+  jac(1, 1:m)=[-20.0*p(1), 10.0];
+  jac(2, 1:m)=[-1.0, 0.0];
Index: ThirdParty/levmar/matlab/jacmeyer.m
===================================================================
--- ThirdParty/levmar/matlab/jacmeyer.m	(revision 8ce1a95e94bdc842511db41cb292238a19558df0)
+++ ThirdParty/levmar/matlab/jacmeyer.m	(revision 8ce1a95e94bdc842511db41cb292238a19558df0)
@@ -0,0 +1,10 @@
+function jac = jacmeyer(p, data1, data2)
+  n=16;
+  m=3;
+
+  for i=1:n
+    ui=0.45+0.05*i;
+    tmp=exp(10.0*p(2)/(ui+p(3)) - 13.0);
+
+    jac(i, 1:m)=[tmp, 10.0*p(1)*tmp/(ui+p(3)), -10.0*p(1)*p(2)*tmp/((ui+p(3))*(ui+p(3)))];
+  end
Index: ThirdParty/levmar/matlab/jacmodhs38.m
===================================================================
--- ThirdParty/levmar/matlab/jacmodhs38.m	(revision 8ce1a95e94bdc842511db41cb292238a19558df0)
+++ ThirdParty/levmar/matlab/jacmodhs38.m	(revision 8ce1a95e94bdc842511db41cb292238a19558df0)
@@ -0,0 +1,9 @@
+function jac = jacmodhs38(p)
+  m=4;
+
+  jac(1, 1:m)=[-20.0*p(1), 10.0, 0.0, 0.0];
+  jac(2, 1:m)=[-1.0, 0.0, 0.0, 0.0];
+  jac(3, 1:m)=[0.0, 0.0, -2.0*sqrt(90.0)*p(3), sqrt(90.0)];
+  jac(4, 1:m)=[0.0, 0.0, -1.0, 0.0];
+  jac(5, 1:m)=[0.0, sqrt(10.1), 0.0, 0.0];
+  jac(6, 1:m)=[0.0, 0.0, 0.0, sqrt(10.1)];
Index: ThirdParty/levmar/matlab/jacmodhs52.m
===================================================================
--- ThirdParty/levmar/matlab/jacmodhs52.m	(revision 8ce1a95e94bdc842511db41cb292238a19558df0)
+++ ThirdParty/levmar/matlab/jacmodhs52.m	(revision 8ce1a95e94bdc842511db41cb292238a19558df0)
@@ -0,0 +1,7 @@
+function jac = jacmodhs52(p)
+  m=5;
+
+  jac(1, 1:m)=[4.0, -1.0, 0.0, 0.0, 0.0];
+  jac(2, 1:m)=[0.0, 1.0, 1.0, 0.0, 0.0];
+  jac(3, 1:m)=[0.0, 0.0, 0.0, 1.0, 0.0];
+  jac(4, 1:m)=[0.0, 0.0, 0.0, 0.0, 1.0];
Index: ThirdParty/levmar/matlab/jacmodhs76.m
===================================================================
--- ThirdParty/levmar/matlab/jacmodhs76.m	(revision 8ce1a95e94bdc842511db41cb292238a19558df0)
+++ ThirdParty/levmar/matlab/jacmodhs76.m	(revision 8ce1a95e94bdc842511db41cb292238a19558df0)
@@ -0,0 +1,7 @@
+function jac = jacmodhs76(p)
+  m=4;
+
+  jac(1, 1:m)=[1.0, 0.0, 0.0, 0.0];
+  jac(2, 1:m)=[0.0, sqrt(0.5), 0.0, 0.0];
+  jac(3, 1:m)=[0.0, 0.0, 1.0, 0.0];
+  jac(4, 1:m)=[0.0, 0.0, 0.0, sqrt(0.5)];
Index: ThirdParty/levmar/matlab/jacosborne.m
===================================================================
--- ThirdParty/levmar/matlab/jacosborne.m	(revision 8ce1a95e94bdc842511db41cb292238a19558df0)
+++ ThirdParty/levmar/matlab/jacosborne.m	(revision 8ce1a95e94bdc842511db41cb292238a19558df0)
@@ -0,0 +1,11 @@
+function jac = jacosborne(p)
+  n=33;
+  m=5;
+
+  for i=1:n
+    t=10*(i-1);
+    tmp1=exp(-p(4)*t);
+    tmp2=exp(-p(5)*t);
+
+    jac(i, 1:m)=[1.0, tmp1, tmp2, -p(2)*t*tmp1, -p(3)*t*tmp2];
+  end
Index: ThirdParty/levmar/matlab/levmar.c
===================================================================
--- ThirdParty/levmar/matlab/levmar.c	(revision 8ce1a95e94bdc842511db41cb292238a19558df0)
+++ ThirdParty/levmar/matlab/levmar.c	(revision 8ce1a95e94bdc842511db41cb292238a19558df0)
@@ -0,0 +1,692 @@
+/* ////////////////////////////////////////////////////////////////////////////////
+// 
+//  Matlab MEX file for the Levenberg - Marquardt minimization algorithm
+//  Copyright (C) 2007  Manolis Lourakis (lourakis at ics forth gr)
+//  Institute of Computer Science, Foundation for Research & Technology - Hellas
+//  Heraklion, Crete, Greece.
+//
+//  This program is free software; you can redistribute it and/or modify
+//  it under the terms of the GNU General Public License as published by
+//  the Free Software Foundation; either version 2 of the License, or
+//  (at your option) any later version.
+//
+//  This program is distributed in the hope that it will be useful,
+//  but WITHOUT ANY WARRANTY; without even the implied warranty of
+//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+//  GNU General Public License for more details.
+//
+//////////////////////////////////////////////////////////////////////////////// */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <stdarg.h>
+#include <math.h>
+#include <string.h>
+#include <ctype.h>
+
+#include <levmar.h>
+
+#include <mex.h>
+
+/**
+#define DEBUG
+**/
+
+#ifndef HAVE_LAPACK
+#ifdef _MSC_VER
+#pragma message("LAPACK not available, certain functionalities cannot be compiled!")
+#else
+#warning LAPACK not available, certain functionalities cannot be compiled
+#endif /* _MSC_VER */
+#endif /* HAVE_LAPACK */
+
+#define __MAX__(A, B)     ((A)>=(B)? (A) : (B))
+
+#define MIN_UNCONSTRAINED     0
+#define MIN_CONSTRAINED_BC    1
+#define MIN_CONSTRAINED_LEC   2
+#define MIN_CONSTRAINED_BLEC  3
+#define MIN_CONSTRAINED_BLEIC 4
+#define MIN_CONSTRAINED_BLIC  5
+#define MIN_CONSTRAINED_LEIC  6
+#define MIN_CONSTRAINED_LIC   7
+
+struct mexdata {
+  /* matlab names of the fitting function & its Jacobian */
+  char *fname, *jacname;
+
+  /* binary flags specifying if input p0 is a row or column vector */
+  int isrow_p0;
+
+  /* rhs args to be passed to matlab. rhs[0] is reserved for
+   * passing the parameter vector. If present, problem-specific
+   * data are passed in rhs[1], rhs[2], etc
+   */
+  mxArray **rhs;
+  int nrhs; /* >= 1 */
+};
+
+/* display printf-style error messages in matlab */
+static void matlabFmtdErrMsgTxt(char *fmt, ...)
+{
+char  buf[256];
+va_list args;
+
+  va_start(args, fmt);
+  vsprintf(buf, fmt, args);
+  va_end(args);
+
+  mexErrMsgTxt(buf);
+}
+
+/* display printf-style warning messages in matlab */
+static void matlabFmtdWarnMsgTxt(char *fmt, ...)
+{
+char  buf[256];
+va_list args;
+
+  va_start(args, fmt);
+  vsprintf(buf, fmt, args);
+  va_end(args);
+
+  mexWarnMsgTxt(buf);
+}
+
+static void func(double *p, double *hx, int m, int n, void *adata)
+{
+mxArray *lhs[1];
+double *mp, *mx;
+register int i;
+struct mexdata *dat=(struct mexdata *)adata;
+    
+  /* prepare to call matlab */
+  mp=mxGetPr(dat->rhs[0]);
+  for(i=0; i<m; ++i)
+    mp[i]=p[i];
+    
+  /* invoke matlab */
+  mexCallMATLAB(1, lhs, dat->nrhs, dat->rhs, dat->fname);
+
+  /* copy back results & cleanup */
+  mx=mxGetPr(lhs[0]);
+  for(i=0; i<n; ++i)
+    hx[i]=mx[i];
+
+  /* delete the matrix created by matlab */
+  mxDestroyArray(lhs[0]);
+}
+
+static void jacfunc(double *p, double *j, int m, int n, void *adata)
+{
+mxArray *lhs[1];
+double *mp;
+double *mj;
+register int i, k;
+struct mexdata *dat=(struct mexdata *)adata;
+    
+  /* prepare to call matlab */
+  mp=mxGetPr(dat->rhs[0]);
+  for(i=0; i<m; ++i)
+    mp[i]=p[i];
+
+  /* invoke matlab */
+  mexCallMATLAB(1, lhs, dat->nrhs, dat->rhs, dat->jacname);
+    
+  /* copy back results & cleanup. Note that the nxm Jacobian 
+   * computed by matlab should be transposed so that
+   * levmar gets it in row major, as expected
+   */
+  mj=mxGetPr(lhs[0]);
+  for(i=0; i<n; ++i)
+    for(k=0; k<m; ++k)
+      j[i*m+k]=mj[i+k*n];
+
+  /* delete the matrix created by matlab */
+  mxDestroyArray(lhs[0]);
+}
+
+/* matlab matrices are in column-major, this routine converts them to row major for levmar */
+static double *getTranspose(mxArray *Am)
+{
+int m, n;
+double *At, *A;
+register int i, j;
+
+  m=mxGetM(Am);
+  n=mxGetN(Am);
+  A=mxGetPr(Am);
+  At=mxMalloc(m*n*sizeof(double));
+
+  for(i=0; i<m; i++)
+    for(j=0; j<n; j++)
+      At[i*n+j]=A[i+j*m];
+  
+  return At;
+}
+
+/* check the supplied matlab function and its Jacobian. Returns 1 on error, 0 otherwise */
+static int checkFuncAndJacobian(double *p, int  m, int n, int havejac, struct mexdata *dat)
+{
+mxArray *lhs[1];
+register int i;
+int ret=0;
+double *mp;
+
+  mexSetTrapFlag(1); /* handle errors in the MEX-file */
+
+  mp=mxGetPr(dat->rhs[0]);
+  for(i=0; i<m; ++i)
+    mp[i]=p[i];
+
+  /* attempt to call the supplied func */
+  i=mexCallMATLAB(1, lhs, dat->nrhs, dat->rhs, dat->fname);
+  if(i){
+    fprintf(stderr, "levmar: error calling '%s'.\n", dat->fname);
+    ret=1;
+  }
+  else if(!mxIsDouble(lhs[0]) || mxIsComplex(lhs[0]) || !(mxGetM(lhs[0])==1 || mxGetN(lhs[0])==1) ||
+      __MAX__(mxGetM(lhs[0]), mxGetN(lhs[0]))!=n){
+    fprintf(stderr, "levmar: '%s' should produce a real vector with %d elements (got %d).\n",
+                    dat->fname, n, __MAX__(mxGetM(lhs[0]), mxGetN(lhs[0])));
+    ret=1;
+  }
+  /* delete the matrix created by matlab */
+  mxDestroyArray(lhs[0]);
+
+  if(havejac){
+    /* attempt to call the supplied jac  */
+    i=mexCallMATLAB(1, lhs, dat->nrhs, dat->rhs, dat->jacname);
+    if(i){
+      fprintf(stderr, "levmar: error calling '%s'.\n", dat->jacname);
+      ret=1;
+    }
+    else if(!mxIsDouble(lhs[0]) || mxIsComplex(lhs[0]) || mxGetM(lhs[0])!=n || mxGetN(lhs[0])!=m){
+      fprintf(stderr, "levmar: '%s' should produce a real %dx%d matrix (got %dx%d).\n",
+                      dat->jacname, n, m, mxGetM(lhs[0]), mxGetN(lhs[0]));
+      ret=1;
+    }
+    else if(mxIsSparse(lhs[0])){
+      fprintf(stderr, "levmar: '%s' should produce a real dense matrix (got a sparse one).\n", dat->jacname);
+      ret=1;
+    }
+    /* delete the matrix created by matlab */
+    mxDestroyArray(lhs[0]);
+  }
+
+  mexSetTrapFlag(0); /* on error terminate the MEX-file and return control to the MATLAB prompt */
+
+  return ret;
+}
+
+
+/*
+[ret, p, info, covar]=levmar_der (f, j, p0, x, itmax, opts, 'unc'                              ...)
+[ret, p, info, covar]=levmar_bc  (f, j, p0, x, itmax, opts, 'bc',   lb, ub,                    ...)
+[ret, p, info, covar]=levmar_bc  (f, j, p0, x, itmax, opts, 'bc',   lb, ub, dscl,              ...)
+[ret, p, info, covar]=levmar_lec (f, j, p0, x, itmax, opts, 'lec',                A, b,        ...)
+[ret, p, info, covar]=levmar_blec(f, j, p0, x, itmax, opts, 'blec', lb, ub,       A, b, wghts, ...)
+
+[ret, p, info, covar]=levmar_bleic(f, j, p0, x, itmax, opts, 'bleic', lb, ub,       A, b, C, d, ...)
+[ret, p, info, covar]=levmar_blic (f, j, p0, x, itmax, opts, 'blic',  lb, ub,             C, d, ...)
+[ret, p, info, covar]=levmar_leic (f, j, p0, x, itmax, opts, 'leic',                A, b, C, d, ...)
+[ret, p, info, covar]=levmar_lic  (f, j, p0, x, itmax, opts, 'lic',                       C, d, ...)
+
+*/
+
+void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *Prhs[])
+{
+register int i;
+register double *pdbl;
+mxArray **prhs=(mxArray **)&Prhs[0], *At, *Ct;
+struct mexdata mdata;
+int len, status;
+double *p, *p0, *ret, *x;
+int m, n, havejac, Arows, Crows, itmax, nopts, mintype, nextra;
+double opts[LM_OPTS_SZ]={LM_INIT_MU, LM_STOP_THRESH, LM_STOP_THRESH, LM_STOP_THRESH, LM_DIFF_DELTA};
+double info[LM_INFO_SZ];
+double *lb=NULL, *ub=NULL, *dscl=NULL, *A=NULL, *b=NULL, *wghts=NULL, *C=NULL, *d=NULL, *covar=NULL;
+
+  /* parse input args; start by checking their number */
+  if((nrhs<5))
+    matlabFmtdErrMsgTxt("levmar: at least 5 input arguments required (got %d).", nrhs);
+  if(nlhs>4)
+    matlabFmtdErrMsgTxt("levmar: too many output arguments (max. 4, got %d).", nlhs);
+  else if(nlhs<2)
+    matlabFmtdErrMsgTxt("levmar: too few output arguments (min. 2, got %d).", nlhs);
+    
+  /* note that in order to accommodate optional args, prhs & nrhs are adjusted accordingly below */
+
+  /** func **/
+  /* first argument must be a string , i.e. a char row vector */
+  if(mxIsChar(prhs[0])!=1)
+    mexErrMsgTxt("levmar: first argument must be a string.");
+  if(mxGetM(prhs[0])!=1)
+    mexErrMsgTxt("levmar: first argument must be a string (i.e. char row vector).");
+  /* store supplied name */
+  len=mxGetN(prhs[0])+1;
+  mdata.fname=mxCalloc(len, sizeof(char));
+  status=mxGetString(prhs[0], mdata.fname, len);
+  if(status!=0)
+    mexErrMsgTxt("levmar: not enough space. String is truncated.");
+
+  /** jac (optional) **/
+  /* check whether second argument is a string */
+  if(mxIsChar(prhs[1])==1){
+    if(mxGetM(prhs[1])!=1)
+      mexErrMsgTxt("levmar: second argument must be a string (i.e. row vector).");
+    /* store supplied name */
+    len=mxGetN(prhs[1])+1;
+    mdata.jacname=mxCalloc(len, sizeof(char));
+    status=mxGetString(prhs[1], mdata.jacname, len);
+    if(status!=0)
+      mexErrMsgTxt("levmar: not enough space. String is truncated.");
+    havejac=1;
+
+    ++prhs;
+    --nrhs;
+  }
+  else{
+    mdata.jacname=NULL;
+    havejac=0;
+  }
+
+#ifdef DEBUG
+  fflush(stderr);
+  fprintf(stderr, "LEVMAR: %s analytic Jacobian\n", havejac? "with" : "no");
+#endif /* DEBUG */
+
+/* CHECK 
+if(!mxIsDouble(prhs[1]) || mxIsComplex(prhs[1]) || !(mxGetM(prhs[1])==1 && mxGetN(prhs[1])==1))
+*/
+
+  /** p0 **/
+  /* the second required argument must be a real row or column vector */
+  if(!mxIsDouble(prhs[1]) || mxIsComplex(prhs[1]) || !(mxGetM(prhs[1])==1 || mxGetN(prhs[1])==1))
+    mexErrMsgTxt("levmar: p0 must be a real vector.");
+  p0=mxGetPr(prhs[1]);
+  /* determine if we have a row or column vector and retrieve its 
+   * size, i.e. the number of parameters
+   */
+  if(mxGetM(prhs[1])==1){
+    m=mxGetN(prhs[1]);
+    mdata.isrow_p0=1;
+  }
+  else{
+    m=mxGetM(prhs[1]);
+    mdata.isrow_p0=0;
+  }
+  /* copy input parameter vector to avoid destroying it */
+  p=mxMalloc(m*sizeof(double));
+  for(i=0; i<m; ++i)
+    p[i]=p0[i];
+    
+  /** x **/
+  /* the third required argument must be a real row or column vector */
+  if(!mxIsDouble(prhs[2]) || mxIsComplex(prhs[2]) || !(mxGetM(prhs[2])==1 || mxGetN(prhs[2])==1))
+    mexErrMsgTxt("levmar: x must be a real vector.");
+  x=mxGetPr(prhs[2]);
+  n=__MAX__(mxGetM(prhs[2]), mxGetN(prhs[2]));
+
+  /** itmax **/
+  /* the fourth required argument must be a scalar */
+  if(!mxIsDouble(prhs[3]) || mxIsComplex(prhs[3]) || mxGetM(prhs[3])!=1 || mxGetN(prhs[3])!=1)
+    mexErrMsgTxt("levmar: itmax must be a scalar.");
+  itmax=(int)mxGetScalar(prhs[3]);
+    
+  /** opts **/
+  /* the fifth required argument must be a real row or column vector */
+  if(!mxIsDouble(prhs[4]) || mxIsComplex(prhs[4]) || (!(mxGetM(prhs[4])==1 || mxGetN(prhs[4])==1) &&
+                                                      !(mxGetM(prhs[4])==0 && mxGetN(prhs[4])==0)))
+    mexErrMsgTxt("levmar: opts must be a real vector.");
+  pdbl=mxGetPr(prhs[4]);
+  nopts=__MAX__(mxGetM(prhs[4]), mxGetN(prhs[4]));
+  if(nopts!=0){ /* if opts==[], nothing needs to be done and the defaults are used */
+    if(nopts>LM_OPTS_SZ)
+      matlabFmtdErrMsgTxt("levmar: opts must have at most %d elements, got %d.", LM_OPTS_SZ, nopts);
+    else if(nopts<((havejac)? LM_OPTS_SZ-1 : LM_OPTS_SZ))
+      matlabFmtdWarnMsgTxt("levmar: only the %d first elements of opts specified, remaining set to defaults.", nopts);
+    for(i=0; i<nopts; ++i)
+      opts[i]=pdbl[i];
+  }
+#ifdef DEBUG
+  else{
+    fflush(stderr);
+    fprintf(stderr, "LEVMAR: empty options vector, using defaults\n");
+  }
+#endif /* DEBUG */
+
+  /** mintype (optional) **/
+  /* check whether sixth argument is a string */
+  if(nrhs>=6 && mxIsChar(prhs[5])==1 && mxGetM(prhs[5])==1){
+    char *minhowto;
+
+    /* examine supplied name */
+    len=mxGetN(prhs[5])+1;
+    minhowto=mxCalloc(len, sizeof(char));
+    status=mxGetString(prhs[5], minhowto, len);
+    if(status!=0)
+      mexErrMsgTxt("levmar: not enough space. String is truncated.");
+
+    for(i=0; minhowto[i]; ++i)
+      minhowto[i]=tolower(minhowto[i]);
+    if(!strncmp(minhowto, "unc", 3)) mintype=MIN_UNCONSTRAINED;
+    else if(!strncmp(minhowto, "bc", 2)) mintype=MIN_CONSTRAINED_BC;
+    else if(!strncmp(minhowto, "lec", 3)) mintype=MIN_CONSTRAINED_LEC;
+    else if(!strncmp(minhowto, "blec", 4)) mintype=MIN_CONSTRAINED_BLEC;
+    else if(!strncmp(minhowto, "bleic", 5)) mintype=MIN_CONSTRAINED_BLEIC;
+    else if(!strncmp(minhowto, "blic", 4)) mintype=MIN_CONSTRAINED_BLIC;
+    else if(!strncmp(minhowto, "leic", 4)) mintype=MIN_CONSTRAINED_LEIC;
+    else if(!strncmp(minhowto, "lic", 3)) mintype=MIN_CONSTRAINED_BLIC;
+    else matlabFmtdErrMsgTxt("levmar: unknown minimization type '%s'.", minhowto);
+
+    mxFree(minhowto);
+
+    ++prhs;
+    --nrhs;
+  }
+  else
+    mintype=MIN_UNCONSTRAINED;
+
+  if(mintype==MIN_UNCONSTRAINED) goto extraargs;
+
+  /* arguments below this point are optional and their presence depends
+   * upon the minimization type determined above
+   */
+  /** lb, ub **/
+  if(nrhs>=7 && (mintype==MIN_CONSTRAINED_BC || mintype==MIN_CONSTRAINED_BLEC || mintype==MIN_CONSTRAINED_BLIC || mintype==MIN_CONSTRAINED_BLEIC)){
+    /* check if the next two arguments are real row or column vectors */
+    if(mxIsDouble(prhs[5]) && !mxIsComplex(prhs[5]) && (mxGetM(prhs[5])==1 || mxGetN(prhs[5])==1)){
+      if(mxIsDouble(prhs[6]) && !mxIsComplex(prhs[6]) && (mxGetM(prhs[6])==1 || mxGetN(prhs[6])==1)){
+        if((i=__MAX__(mxGetM(prhs[5]), mxGetN(prhs[5])))!=m)
+          matlabFmtdErrMsgTxt("levmar: lb must have %d elements, got %d.", m, i);
+        if((i=__MAX__(mxGetM(prhs[6]), mxGetN(prhs[6])))!=m)
+          matlabFmtdErrMsgTxt("levmar: ub must have %d elements, got %d.", m, i);
+
+        lb=mxGetPr(prhs[5]);
+        ub=mxGetPr(prhs[6]);
+
+        prhs+=2;
+        nrhs-=2;
+      }
+    }
+  }
+
+  /** dscl **/
+  if(nrhs>=7 && (mintype==MIN_CONSTRAINED_BC)){
+    /* check if the next argument is a real row or column vector */
+    if(mxIsDouble(prhs[5]) && !mxIsComplex(prhs[5]) && (mxGetM(prhs[5])==1 || mxGetN(prhs[5])==1)){
+      if((i=__MAX__(mxGetM(prhs[5]), mxGetN(prhs[5])))!=m)
+        matlabFmtdErrMsgTxt("levmar: dscl must have %d elements, got %d.", m, i);
+
+      dscl=mxGetPr(prhs[5]);
+
+      ++prhs;
+      --nrhs;
+    }
+  }
+
+  /** A, b **/
+  if(nrhs>=7 && (mintype==MIN_CONSTRAINED_LEC || mintype==MIN_CONSTRAINED_BLEC || mintype==MIN_CONSTRAINED_LEIC || mintype==MIN_CONSTRAINED_BLEIC)){
+    /* check if the next two arguments are a real matrix and a real row or column vector */
+    if(mxIsDouble(prhs[5]) && !mxIsComplex(prhs[5]) && mxGetM(prhs[5])>=1 && mxGetN(prhs[5])>=1){
+      if(mxIsDouble(prhs[6]) && !mxIsComplex(prhs[6]) && (mxGetM(prhs[6])==1 || mxGetN(prhs[6])==1)){
+        if((i=mxGetN(prhs[5]))!=m)
+          matlabFmtdErrMsgTxt("levmar: A must have %d columns, got %d.", m, i);
+        if((i=__MAX__(mxGetM(prhs[6]), mxGetN(prhs[6])))!=(Arows=mxGetM(prhs[5])))
+          matlabFmtdErrMsgTxt("levmar: b must have %d elements, got %d.", Arows, i);
+
+        At=prhs[5];
+        b=mxGetPr(prhs[6]);
+        A=getTranspose(At);
+
+        prhs+=2;
+        nrhs-=2;
+      }
+    }
+  }
+
+  /* wghts */
+  /* check if we have a weights vector */
+  if(nrhs>=6 && mintype==MIN_CONSTRAINED_BLEC){ /* only check if we have seen both box & linear constraints */
+    if(mxIsDouble(prhs[5]) && !mxIsComplex(prhs[5]) && (mxGetM(prhs[5])==1 || mxGetN(prhs[5])==1)){
+      if(__MAX__(mxGetM(prhs[5]), mxGetN(prhs[5]))==m){
+        wghts=mxGetPr(prhs[5]);
+
+        ++prhs;
+        --nrhs;
+      }
+    }
+  }
+
+  /** C, d **/
+  if(nrhs>=7 && (mintype==MIN_CONSTRAINED_BLEIC || mintype==MIN_CONSTRAINED_BLIC || mintype==MIN_CONSTRAINED_LEIC || mintype==MIN_CONSTRAINED_LIC)){
+    /* check if the next two arguments are a real matrix and a real row or column vector */
+    if(mxIsDouble(prhs[5]) && !mxIsComplex(prhs[5]) && mxGetM(prhs[5])>=1 && mxGetN(prhs[5])>=1){
+      if(mxIsDouble(prhs[6]) && !mxIsComplex(prhs[6]) && (mxGetM(prhs[6])==1 || mxGetN(prhs[6])==1)){
+        if((i=mxGetN(prhs[5]))!=m)
+          matlabFmtdErrMsgTxt("levmar: C must have %d columns, got %d.", m, i);
+        if((i=__MAX__(mxGetM(prhs[6]), mxGetN(prhs[6])))!=(Crows=mxGetM(prhs[5])))
+          matlabFmtdErrMsgTxt("levmar: d must have %d elements, got %d.", Crows, i);
+
+        Ct=prhs[5];
+        d=mxGetPr(prhs[6]);
+        C=getTranspose(Ct);
+
+        prhs+=2;
+        nrhs-=2;
+      }
+    }
+  }
+
+  /* arguments below this point are assumed to be extra arguments passed
+   * to every invocation of the fitting function and its Jacobian
+   */
+
+extraargs:
+  /* handle any extra args and allocate memory for
+   * passing the current parameter estimate to matlab
+   */
+  nextra=nrhs-5;
+  mdata.nrhs=nextra+1;
+  mdata.rhs=(mxArray **)mxMalloc(mdata.nrhs*sizeof(mxArray *));
+  for(i=0; i<nextra; ++i)
+    mdata.rhs[i+1]=(mxArray *)prhs[nrhs-nextra+i]; /* discard 'const' modifier */
+#ifdef DEBUG
+  fflush(stderr);
+  fprintf(stderr, "LEVMAR: %d extra args\n", nextra);
+#endif /* DEBUG */
+
+  if(mdata.isrow_p0){ /* row vector */
+    mdata.rhs[0]=mxCreateDoubleMatrix(1, m, mxREAL);
+    /*
+    mxSetM(mdata.rhs[0], 1);
+    mxSetN(mdata.rhs[0], m);
+    */
+  }
+  else{ /* column vector */
+    mdata.rhs[0]=mxCreateDoubleMatrix(m, 1, mxREAL);
+    /*
+    mxSetM(mdata.rhs[0], m);
+    mxSetN(mdata.rhs[0], 1);
+    */
+  }
+
+  /* ensure that the supplied function & Jacobian are as expected */
+  if(checkFuncAndJacobian(p, m, n, havejac, &mdata)){
+    status=LM_ERROR;
+    goto cleanup;
+  }
+
+  if(nlhs>3) /* covariance output required */
+    covar=mxMalloc(m*m*sizeof(double));
+
+  /* invoke levmar */
+  switch(mintype){
+    case MIN_UNCONSTRAINED: /* no constraints */
+      if(havejac)
+        status=dlevmar_der(func, jacfunc, p, x, m, n, itmax, opts, info, NULL, covar, (void *)&mdata);
+      else
+        status=dlevmar_dif(func,          p, x, m, n, itmax, opts, info, NULL, covar, (void *)&mdata);
+#ifdef DEBUG
+  fflush(stderr);
+  fprintf(stderr, "LEVMAR: calling dlevmar_der()/dlevmar_dif()\n");
+#endif /* DEBUG */
+    break;
+    case MIN_CONSTRAINED_BC: /* box constraints */
+      if(havejac)
+        status=dlevmar_bc_der(func, jacfunc, p, x, m, n, lb, ub, dscl, itmax, opts, info, NULL, covar, (void *)&mdata);
+      else
+        status=dlevmar_bc_dif(func,          p, x, m, n, lb, ub, dscl, itmax, opts, info, NULL, covar, (void *)&mdata);
+#ifdef DEBUG
+  fflush(stderr);
+  fprintf(stderr, "LEVMAR: calling dlevmar_bc_der()/dlevmar_bc_dif()\n");
+#endif /* DEBUG */
+    break;
+    case MIN_CONSTRAINED_LEC:  /* linear equation constraints */
+#ifdef HAVE_LAPACK
+      if(havejac)
+        status=dlevmar_lec_der(func, jacfunc, p, x, m, n, A, b, Arows, itmax, opts, info, NULL, covar, (void *)&mdata);
+      else
+        status=dlevmar_lec_dif(func,          p, x, m, n, A, b, Arows, itmax, opts, info, NULL, covar, (void *)&mdata);
+#else
+      mexErrMsgTxt("levmar: no linear constraints support, HAVE_LAPACK was not defined during MEX-file compilation.");
+#endif /* HAVE_LAPACK */
+
+#ifdef DEBUG
+  fflush(stderr);
+  fprintf(stderr, "LEVMAR: calling dlevmar_lec_der()/dlevmar_lec_dif()\n");
+#endif /* DEBUG */
+    break;
+    case MIN_CONSTRAINED_BLEC: /* box & linear equation constraints */
+#ifdef HAVE_LAPACK
+      if(havejac)
+        status=dlevmar_blec_der(func, jacfunc, p, x, m, n, lb, ub, A, b, Arows, wghts, itmax, opts, info, NULL, covar, (void *)&mdata);
+      else
+        status=dlevmar_blec_dif(func,          p, x, m, n, lb, ub, A, b, Arows, wghts, itmax, opts, info, NULL, covar, (void *)&mdata);
+#else
+      mexErrMsgTxt("levmar: no box & linear constraints support, HAVE_LAPACK was not defined during MEX-file compilation.");
+#endif /* HAVE_LAPACK */
+
+#ifdef DEBUG
+  fflush(stderr);
+  fprintf(stderr, "LEVMAR: calling dlevmar_blec_der()/dlevmar_blec_dif()\n");
+#endif /* DEBUG */
+    break;
+    case MIN_CONSTRAINED_BLEIC: /* box, linear equation & inequalities constraints */
+#ifdef HAVE_LAPACK
+      if(havejac)
+        status=dlevmar_bleic_der(func, jacfunc, p, x, m, n, lb, ub, A, b, Arows, C, d, Crows, itmax, opts, info, NULL, covar, (void *)&mdata);
+      else
+        status=dlevmar_bleic_dif(func, p, x, m, n, lb, ub, A, b, Arows, C, d, Crows, itmax, opts, info, NULL, covar, (void *)&mdata);
+#else
+      mexErrMsgTxt("levmar: no box, linear equation & inequality constraints support, HAVE_LAPACK was not defined during MEX-file compilation.");
+#endif /* HAVE_LAPACK */
+
+#ifdef DEBUG
+  fflush(stderr);
+  fprintf(stderr, "LEVMAR: calling dlevmar_bleic_der()/dlevmar_bleic_dif()\n");
+#endif /* DEBUG */
+    break;
+    case MIN_CONSTRAINED_BLIC: /* box, linear inequalities constraints */
+#ifdef HAVE_LAPACK
+      if(havejac)
+        status=dlevmar_bleic_der(func, jacfunc, p, x, m, n, lb, ub, NULL, NULL, 0, C, d, Crows, itmax, opts, info, NULL, covar, (void *)&mdata);
+      else
+        status=dlevmar_bleic_dif(func, p, x, m, n, lb, ub, NULL, NULL, 0, C, d, Crows, itmax, opts, info, NULL, covar, (void *)&mdata);
+#else
+      mexErrMsgTxt("levmar: no box & linear inequality constraints support, HAVE_LAPACK was not defined during MEX-file compilation.");
+#endif /* HAVE_LAPACK */
+
+#ifdef DEBUG
+  fflush(stderr);
+  fprintf(stderr, "LEVMAR: calling dlevmar_blic_der()/dlevmar_blic_dif()\n");
+#endif /* DEBUG */
+    break;
+    case MIN_CONSTRAINED_LEIC: /* linear equation & inequalities constraints */
+#ifdef HAVE_LAPACK
+      if(havejac)
+        status=dlevmar_bleic_der(func, jacfunc, p, x, m, n, NULL, NULL, A, b, Arows, C, d, Crows, itmax, opts, info, NULL, covar, (void *)&mdata);
+      else
+        status=dlevmar_bleic_dif(func, p, x, m, n, NULL, NULL, A, b, Arows, C, d, Crows, itmax, opts, info, NULL, covar, (void *)&mdata);
+#else
+      mexErrMsgTxt("levmar: no linear equation & inequality constraints support, HAVE_LAPACK was not defined during MEX-file compilation.");
+#endif /* HAVE_LAPACK */
+
+#ifdef DEBUG
+  fflush(stderr);
+  fprintf(stderr, "LEVMAR: calling dlevmar_leic_der()/dlevmar_leic_dif()\n");
+#endif /* DEBUG */
+    break;
+    case MIN_CONSTRAINED_LIC: /* linear inequalities constraints */
+#ifdef HAVE_LAPACK
+      if(havejac)
+        status=dlevmar_bleic_der(func, jacfunc, p, x, m, n, NULL, NULL, NULL, NULL, 0, C, d, Crows, itmax, opts, info, NULL, covar, (void *)&mdata);
+      else
+        status=dlevmar_bleic_dif(func, p, x, m, n, NULL, NULL, NULL, NULL, 0, C, d, Crows, itmax, opts, info, NULL, covar, (void *)&mdata);
+#else
+      mexErrMsgTxt("levmar: no linear equation & inequality constraints support, HAVE_LAPACK was not defined during MEX-file compilation.");
+#endif /* HAVE_LAPACK */
+
+#ifdef DEBUG
+  fflush(stderr);
+  fprintf(stderr, "LEVMAR: calling dlevmar_lic_der()/dlevmar_lic_dif()\n");
+#endif /* DEBUG */
+    break;
+    default:
+      mexErrMsgTxt("levmar: unexpected internal error.");
+  }
+
+#ifdef DEBUG
+  fflush(stderr);
+  printf("LEVMAR: minimization returned %d in %g iter, reason %g\n\tSolution: ", status, info[5], info[6]);
+  for(i=0; i<m; ++i)
+    printf("%.7g ", p[i]);
+  printf("\n\n\tMinimization info:\n\t");
+  for(i=0; i<LM_INFO_SZ; ++i)
+    printf("%g ", info[i]);
+  printf("\n");
+#endif /* DEBUG */
+
+  /* copy back return results */
+  /** ret **/
+  plhs[0]=mxCreateDoubleMatrix(1, 1, mxREAL);
+  ret=mxGetPr(plhs[0]);
+  ret[0]=(double)status;
+
+  /** popt **/
+  plhs[1]=(mdata.isrow_p0==1)? mxCreateDoubleMatrix(1, m, mxREAL) : mxCreateDoubleMatrix(m, 1, mxREAL);
+  pdbl=mxGetPr(plhs[1]);
+  for(i=0; i<m; ++i)
+    pdbl[i]=p[i];
+
+  /** info **/
+  if(nlhs>2){
+    plhs[2]=mxCreateDoubleMatrix(1, LM_INFO_SZ, mxREAL);
+    pdbl=mxGetPr(plhs[2]);
+    for(i=0; i<LM_INFO_SZ; ++i)
+      pdbl[i]=info[i];
+  }
+
+  /** covar **/
+  if(nlhs>3){
+    plhs[3]=mxCreateDoubleMatrix(m, m, mxREAL);
+    pdbl=mxGetPr(plhs[3]);
+    for(i=0; i<m*m; ++i) /* covariance matrices are symmetric, thus no need to transpose! */
+      pdbl[i]=covar[i];
+  }
+
+cleanup:
+  /* cleanup */
+  mxDestroyArray(mdata.rhs[0]);
+  if(A) mxFree(A);
+  if(C) mxFree(C);
+
+  mxFree(mdata.fname);
+  if(havejac) mxFree(mdata.jacname);
+  mxFree(p);
+  mxFree(mdata.rhs);
+  if(covar) mxFree(covar);
+
+  if(status==LM_ERROR)
+    mexWarnMsgTxt("levmar: optimization returned with an error!");
+}
Index: ThirdParty/levmar/matlab/levmar.m
===================================================================
--- ThirdParty/levmar/matlab/levmar.m	(revision 8ce1a95e94bdc842511db41cb292238a19558df0)
+++ ThirdParty/levmar/matlab/levmar.m	(revision 8ce1a95e94bdc842511db41cb292238a19558df0)
@@ -0,0 +1,83 @@
+% LEVMAR  matlab MEX interface to the levmar non-linear least squares minimization
+% library available from http://www.ics.forth.gr/~lourakis/levmar/
+% 
+% Usage: levmar can be used in any of the 8 following ways:
+% [ret, popt, info, covar]=levmar(fname, jacname, p0, x, itmax, opts, 'unc', ...)
+% [ret, popt, info, covar]=levmar(fname, jacname, p0, x, itmax, opts, 'bc', lb, ub, ...)
+% [ret, popt, info, covar]=levmar(fname, jacname, p0, x, itmax, opts, 'lec', A, b, ...)
+% [ret, popt, info, covar]=levmar(fname, jacname, p0, x, itmax, opts, 'blec', lb, ub, A, b, wghts, ...)
+%
+% [ret, popt, info, covar]=levmar(fname, jacname, p0, x, itmax, opts, 'bleic', lb, ub, A, b, C, d, ...)
+% [ret, popt, info, covar]=levmar(fname, jacname, p0, x, itmax, opts, 'blic', lb, ub, C, d, ...)
+% [ret, popt, info, covar]=levmar(fname, jacname, p0, x, itmax, opts, 'leic', A, b, C, d, ...)
+% [ret, popt, info, covar]=levmar(fname, jacname, p0, x, itmax, opts, 'lic', C, d, ...)
+%
+%  
+% The dots at the end denote any additional, problem specific data that are passed uninterpreted to
+% all invocations of fname and jacname, see below for details.
+%
+% In the following, the word "vector" is meant to imply either a row or a column vector.
+%
+% required input arguments:
+% - fname: String defining the name of a matlab function implementing the function to be minimized.
+%      fname will be called as fname(p, ...), where p denotes the parameter vector and the dots any
+%      additional data passed as extra arguments during the invocation of levmar (refer to Meyer's
+%      problem in lmdemo.m for an example).
+%
+% - p0: vector of doubles holding the initial parameters estimates.
+%
+% - x: vector of doubles holding the measurements vector.
+%
+% - itmax: maximum number of iterations.
+%
+% - opts: vector of doubles specifying the minimization parameters, as follows:
+%      opts(1) scale factor for the initial damping factor
+%      opts(2) stopping threshold for ||J^T e||_inf
+%      opts(3) stopping threshold for ||Dp||_2
+%      opts(4) stopping threshold for ||e||_2
+%      opts(5) step used in finite difference approximation to the Jacobian.
+%      If an empty vector (i.e. []) is specified, defaults are used.
+%  
+% optional input arguments:
+% - jacname: String defining the name of matlab function implementing the Jacobian of function fname.
+%      jacname will be called as jacname(p, ...) where p is again the parameter vector and the dots
+%      denote any additional data passed as extra arguments to the invocation of levmar. If omitted,
+%      the Jacobian is approximated with finite differences through repeated invocations of fname.
+%
+% - type: String defining the minimization type. It should be one of the following:
+%      'unc' specifies unconstrained minimization.
+%      'bc' specifies minimization subject to box constraints.
+%      'lec' specifies minimization subject to linear equation constraints.
+%      'blec' specifies minimization subject to box and linear equation constraints.
+%      'bleic' specifies minimization subject to box, linear equation and inequality constraints.
+%      'blic' specifies minimization subject to box and linear inequality constraints.
+%      'leic' specifies minimization subject to linear equation and inequality constraints.
+%      'lic' specifies minimization subject to linear inequality constraints.
+%      If omitted, a default of 'unc' is assumed. Depending on the minimization type, the MEX
+%      interface will invoke one of dlevmar_XXX, dlevmar_bc_XXX, dlevmar_lec_XXX, dlevmar_blec_XXX or dlevmar_bleic_XXX
+%
+% - lb, ub: vectors of doubles specifying lower and upper bounds for p, respectively
+%
+% - A, b: k x m matrix and k vector specifying linear equation constraints for p, i.e. A*p=b
+%      A should have full rank.
+%
+% - C, d: k x m matrix and k vector specifying linear inequality constraints for p, i.e. C*p>=d
+%      A should have full rank.
+%
+% - wghts: vector of doubles specifying the weights for the penalty terms corresponding to
+%      the box constraints, see lmblec_core.c for more details. If omitted and a 'blec' type
+%      minimization is to be carried out, default weights are used.
+%  
+%
+% output arguments
+% - ret: return value of levmar, corresponding to the number of iterations if successful, -1 otherwise.
+%
+% - popt: estimated minimizer, i.e. minimized parameters vector.
+%
+% - info: optional array of doubles, which upon return provides information regarding the minimization.
+%      See lm_core.c for more details.
+%
+% - covar: optional covariance matrix corresponding to the estimated minimizer.
+%
+ 
+error('levmar.m is used only for providing documentation to levmar; make sure that levmar.c has been compiled using mex');
Index: ThirdParty/levmar/matlab/lmdemo.m
===================================================================
--- ThirdParty/levmar/matlab/lmdemo.m	(revision 8ce1a95e94bdc842511db41cb292238a19558df0)
+++ ThirdParty/levmar/matlab/lmdemo.m	(revision 8ce1a95e94bdc842511db41cb292238a19558df0)
@@ -0,0 +1,140 @@
+% Demo program for levmar's MEX-file interface
+% Performs minimization of several test problems
+
+format long;
+
+% Unconstrained minimization
+
+% fitting the exponential model x_i=p(1)*exp(-p(2)*i)+p(3) of expfit.c to noisy measurements obtained with (5.0 0.1 1.0)
+p0=[1.0, 0.0, 0.0];
+x=[5.8728, 5.4948, 5.0081, 4.5929, 4.3574, 4.1198, 3.6843, 3.3642, 2.9742, 3.0237, 2.7002, 2.8781,...
+   2.5144, 2.4432, 2.2894, 2.0938, 1.9265, 2.1271, 1.8387, 1.7791, 1.6686, 1.6232, 1.571, 1.6057,...
+   1.3825, 1.5087, 1.3624, 1.4206, 1.2097, 1.3129, 1.131, 1.306, 1.2008, 1.3469, 1.1837, 1.2102,...
+   0.96518, 1.2129, 1.2003, 1.0743];
+
+options=[1E-03, 1E-15, 1E-15, 1E-20, 1E-06];
+% arg demonstrates additional data passing to expfit/jacexpfit
+arg=[40];
+
+[ret, popt, info]=levmar('expfit', 'jacexpfit', p0, x, 200, options, arg);
+disp('Exponential model fitting (see also ../expfit.c)');
+popt
+
+
+% Meyer's (reformulated) problem
+p0=[8.85, 4.0, 2.5];
+
+x=[];
+x(1:4)=[34.780, 28.610, 23.650, 19.630];
+x(5:8)=[16.370, 13.720, 11.540, 9.744];
+x(9:12)=[8.261, 7.030, 6.005, 5.147];
+x(13:16)=[4.427, 3.820, 3.307, 2.872];
+
+options=[1E-03, 1E-15, 1E-15, 1E-20, 1E-06];
+% arg1, arg2 demonstrate additional dummy data passing to meyer/jacmeyer
+arg1=[17];
+arg2=[27];
+
+%[ret, popt, info]=levmar('meyer', 'jacmeyer', p0, x, 200, options, arg1, arg2);
+
+%[ret, popt, info, covar]=levmar('meyer', 'jacmeyer', p0, x, 200, options, arg1, arg2);
+[ret, popt, info, covar]=levmar('meyer', p0, x, 200, options, 'unc', arg1, arg2);
+disp('Meyer''s (reformulated) problem');
+popt
+
+
+% Osborne's problem
+p0=[0.5, 1.5, -1.0, 1.0E-2, 2.0E-2];
+
+x=[8.44E-1, 9.08E-1, 9.32E-1, 9.36E-1, 9.25E-1, 9.08E-1, 8.81E-1,...
+8.5E-1, 8.18E-1, 7.84E-1, 7.51E-1, 7.18E-1, 6.85E-1, 6.58E-1,...
+6.28E-1, 6.03E-1, 5.8E-1, 5.58E-1, 5.38E-1, 5.22E-1, 5.06E-1,...
+4.9E-1, 4.78E-1, 4.67E-1, 4.57E-1, 4.48E-1, 4.38E-1, 4.31E-1,...
+4.24E-1, 4.2E-1, 4.14E-1, 4.11E-1, 4.06E-1];
+
+
+options=[1E-03, 1E-15, 1E-15, 1E-20, 1E-06];
+
+[ret, popt, info, covar]=levmar('osborne', 'jacosborne', p0, x, 200, options);
+%[ret, popt, info, covar]=levmar('osborne', p0, x, 200, options, 'unc');
+disp('Osborne''s problem');
+popt
+
+
+% Linear constraints
+
+% Boggs-Tolle problem 3
+p0=[2.0, 2.0, 2.0, 2.0, 2.0];
+x=[0.0, 0.0, 0.0, 0.0, 0.0];
+options=[1E-03, 1E-15, 1E-15, 1E-20];
+adata=[];
+
+A=[1.0, 3.0, 0.0, 0.0, 0.0;
+   0.0, 0.0, 1.0, 1.0, -2.0;
+   0.0, 1.0, 0.0, 0.0, -1.0];
+b=[0.0, 0.0, 0.0]';
+
+[ret, popt, info, covar]=levmar('bt3', 'jacbt3', p0, x, 200, options, 'lec', A, b, adata);
+disp('Boggs-Tolle problem 3');
+popt
+
+
+% Box constraints
+
+% Hock-Schittkowski problem 01
+p0=[-2.0, 1.0];
+x=[0.0, 0.0];
+lb=[-realmax, -1.5];
+ub=[realmax, realmax];
+options=[1E-03, 1E-15, 1E-15, 1E-20];
+
+[ret, popt, info, covar]=levmar('hs01', 'jachs01', p0, x, 200, options, 'bc', lb, ub);
+disp('Hock-Schittkowski problem 01');
+popt
+
+
+% Box and linear constraints
+
+% Hock-Schittkowski modified problem 52 (#1)
+p0=[2.0, 2.0, 2.0, 2.0, 2.0];
+x=[0.0, 0.0, 0.0, 0.0];
+lb=[-0.09, 0.0, -realmax, -0.2, 0.0];
+ub=[realmax, 0.3, 0.25, 0.3, 0.3];
+A=[1.0, 3.0, 0.0, 0.0, 0.0;
+   0.0, 0.0, 1.0, 1.0, -2.0;
+   0.0, 1.0, 0.0, 0.0, -1.0];
+b=[0.0, 0.0, 0.0]';
+options=[1E-03, 1E-15, 1E-15, 1E-20];
+
+[ret, popt, info, covar]=levmar('modhs52', 'jacmodhs52', p0, x, 200, options, 'blec', lb, ub, A, b);
+disp('Hock-Schittkowski modified problem 52 (#1)');
+popt
+
+% Schittkowski modified problem 235
+p0=[-2.0, 3.0, 1.0];
+x=[0.0, 0.0];
+lb=[-realmax, 0.1, 0.7];
+ub=[realmax, 2.9, realmax];
+A=[1.0, 0.0, 1.0;
+   0.0, 1.0, -4.0];
+b=[-1.0, 0.0]';
+options=[1E-03, 1E-15, 1E-15, 1E-20];
+
+[ret, popt, info, covar]=levmar('mods235', p0, x, 200, options, 'blec', lb, ub, A, b);
+disp('Hock-Schittkowski modified problem 235');
+popt
+
+% Box, linear equation & inequality constraints
+p0=[0.5, 0.5, 0.5, 0.5];
+x=[0.0, 0.0, 0.0, 0.0];
+lb=[0.0, 0.0, 0.0, 0.0];
+ub=[realmax, realmax, realmax, realmax];
+A=[0.0, 1.0, 4.0, 0.0];
+b=[1.5]';
+C=[-1.0, -2.0, -1.0, -1.0;
+   -3.0, -1.0, -2.0, 1.0];
+d=[-5.0, -0.4]';
+
+[ret, popt, info, covar]=levmar('modhs76', 'jacmodhs76', p0, x, 200, options, 'bleic', lb, ub, A, b, C, d);
+disp('Hock-Schittkowski modified problem 76');
+popt
Index: ThirdParty/levmar/matlab/meyer.m
===================================================================
--- ThirdParty/levmar/matlab/meyer.m	(revision 8ce1a95e94bdc842511db41cb292238a19558df0)
+++ ThirdParty/levmar/matlab/meyer.m	(revision 8ce1a95e94bdc842511db41cb292238a19558df0)
@@ -0,0 +1,9 @@
+function x = meyer(p, data1, data2)
+  n=16;
+
+% data1, data2 are actually unused
+
+  for i=1:n
+    ui=0.45+0.05*i;
+    x(i)=p(1)*exp(10.0*p(2)/(ui+p(3)) - 13.0);
+  end
Index: ThirdParty/levmar/matlab/modhs38.m
===================================================================
--- ThirdParty/levmar/matlab/modhs38.m	(revision 8ce1a95e94bdc842511db41cb292238a19558df0)
+++ ThirdParty/levmar/matlab/modhs38.m	(revision 8ce1a95e94bdc842511db41cb292238a19558df0)
@@ -0,0 +1,8 @@
+function x = modhs38(p)
+
+  x(1)=10.0*(p(2)-p(1)^2);
+  x(2)=(1.0-p(1));
+  x(3)=sqrt(90.0)*(p(4)-p(3)^2);
+  x(4)=(1.0-p(3));
+  x(5)=sqrt(10.1)*(p(2)-1.0);
+  x(6)=sqrt(10.1)*(p(4)-1.0);
Index: ThirdParty/levmar/matlab/modhs52.m
===================================================================
--- ThirdParty/levmar/matlab/modhs52.m	(revision 8ce1a95e94bdc842511db41cb292238a19558df0)
+++ ThirdParty/levmar/matlab/modhs52.m	(revision 8ce1a95e94bdc842511db41cb292238a19558df0)
@@ -0,0 +1,7 @@
+function x = modhs52(p)
+  n=4;
+
+  x(1)=4.0*p(1)-p(2);
+  x(2)=p(2)+p(3)-2.0;
+  x(3)=p(4)-1.0;
+  x(4)=p(5)-1.0;
Index: ThirdParty/levmar/matlab/modhs76.m
===================================================================
--- ThirdParty/levmar/matlab/modhs76.m	(revision 8ce1a95e94bdc842511db41cb292238a19558df0)
+++ ThirdParty/levmar/matlab/modhs76.m	(revision 8ce1a95e94bdc842511db41cb292238a19558df0)
@@ -0,0 +1,7 @@
+function x = modhs76(p)
+
+  x(1)=p(1);
+  x(2)=sqrt(0.5)*p(2);
+  x(3)=p(3);
+  x(4)=sqrt(0.5)*p(4);
+
Index: ThirdParty/levmar/matlab/mods235.m
===================================================================
--- ThirdParty/levmar/matlab/mods235.m	(revision 8ce1a95e94bdc842511db41cb292238a19558df0)
+++ ThirdParty/levmar/matlab/mods235.m	(revision 8ce1a95e94bdc842511db41cb292238a19558df0)
@@ -0,0 +1,5 @@
+function x = mods235(p)
+  n=2;
+
+  x(1)=0.1*(p(1)-1.0);
+  x(2)=p(2)-p(1)*p(1);
Index: ThirdParty/levmar/matlab/osborne.m
===================================================================
--- ThirdParty/levmar/matlab/osborne.m	(revision 8ce1a95e94bdc842511db41cb292238a19558df0)
+++ ThirdParty/levmar/matlab/osborne.m	(revision 8ce1a95e94bdc842511db41cb292238a19558df0)
@@ -0,0 +1,7 @@
+function x = osborne(p)
+  n=33;
+
+  for i=1:n
+    t=10*(i-1);
+    x(i)=p(1) + p(2)*exp(-p(4)*t) + p(3)*exp(-p(5)*t);
+  end
