--- mpfit.c.orig	2014-06-05 10:02:14.819618639 -0500
+++ mpfit.c	2014-06-05 10:02:21.131645649 -0500
@@ -13,7 +13,7 @@
  */
 
 /* Main mpfit library routines (double precision) 
-   $Id: mpfit.c,v 1.20 2010/11/13 08:15:35 craigm Exp $
+   $Id: mpfit.c,v 1.23 2013/04/23 04:24:16 craigm Exp $
  */
 
 #include <stdio.h>
@@ -29,7 +29,8 @@
 	      double *wa, void *priv, int *nfev,
 	      double *step, double *dstep, int *dside,
 	      int *qulimited, double *ulimit,
-	      int *ddebug, double *ddrtol, double *ddatol);
+	      int *ddebug, double *ddrtol, double *ddatol,
+	      double *wa2, double **dvecptr);
 static void mp_qrfac(int m, int n, double *a, int lda, 
 	      int pivot, int *ipvt, int lipvt,
 	      double *rdiag, double *acnorm, double *wa);
@@ -274,7 +275,7 @@
 {
   mp_config conf;
   int i, j, info, iflag, nfree, npegged, iter;
-  int qanylim = 0, qanypegged = 0;
+  int qanylim = 0;
 
   int ij,jj,l;
   double actred,delta,dirder,fnorm,fnorm1,gnorm, orignorm;
@@ -297,6 +298,7 @@
   double *fvec = 0, *qtf = 0;
   double *x = 0, *xnew = 0, *fjac = 0, *diag = 0;
   double *wa1 = 0, *wa2 = 0, *wa3 = 0, *wa4 = 0;
+  double **dvecptr = 0;
   int *ipvt = 0;
 
   int ldfjac;
@@ -323,6 +325,7 @@
     if (config->nprint >= 0) conf.nprint = config->nprint;
     if (config->epsfcn > 0) conf.epsfcn = config->epsfcn;
     if (config->maxiter > 0) conf.maxiter = config->maxiter;
+    if (config->maxiter == MP_NO_ITER) conf.maxiter = 0;
     if (config->douserscale != 0) conf.douserscale = config->douserscale;
     if (config->covtol > 0) conf.covtol = config->covtol;
     if (config->nofinitecheck > 0) conf.nofinitecheck = config->nofinitecheck;
@@ -334,6 +337,7 @@
   nfree = 0;
   npegged = 0;
 
+  /* Basic error checking */
   if (funct == 0) {
     return MP_ERR_FUNC;
   }
@@ -442,6 +446,7 @@
   mp_malloc(wa3, double, npar);
   mp_malloc(wa4, double, m);
   mp_malloc(ipvt, int, npar);
+  mp_malloc(dvecptr, double *, npar);
 
   /* Evaluate user function with initial parameter values */
   iflag = mp_call(funct, m, npar, xall, fvec, 0, private_data);
@@ -483,30 +488,34 @@
   iflag = mp_fdjac2(funct, m, nfree, ifree, npar, xnew, fvec, fjac, ldfjac,
 		    conf.epsfcn, wa4, private_data, &nfev,
 		    step, dstep, mpside, qulim, ulim,
-		    ddebug, ddrtol, ddatol);
+		    ddebug, ddrtol, ddatol, wa2, dvecptr);
   if (iflag < 0) {
     goto CLEANUP;
   }
 
   /* Determine if any of the parameters are pegged at the limits */
-  qanypegged = 0;
   if (qanylim) {
     for (j=0; j<nfree; j++) {
       int lpegged = (qllim[j] && (x[j] == llim[j]));
       int upegged = (qulim[j] && (x[j] == ulim[j]));
       sum = 0;
-      
+
+      /* If the parameter is pegged at a limit, compute the gradient
+	 direction */
       if (lpegged || upegged) {
-	qanypegged = 1;
 	ij = j*ldfjac;
 	for (i=0; i<m; i++, ij++) {
 	  sum += fvec[i] * fjac[ij];
 	}
       }
+      /* If pegged at lower limit and gradient is toward negative then
+	 reset gradient to zero */
       if (lpegged && (sum > 0)) {
 	ij = j*ldfjac;
 	for (i=0; i<m; i++, ij++) fjac[ij] = 0;
       }
+      /* If pegged at upper limit and gradient is toward positive then
+	 reset gradient to zero */
       if (upegged && (sum < 0)) {
 	ij = j*ldfjac;
 	for (i=0; i<m; i++, ij++) fjac[ij] = 0;
@@ -624,7 +633,10 @@
    */
   if (gnorm <= conf.gtol) info = MP_OK_DIR;
   if (info != 0) goto L300;
-  if (conf.maxiter == 0) goto L300;
+  if (conf.maxiter == 0) {
+    info = MP_MAXITER;
+    goto L300;
+  }
 
   /*
    *	 rescale if necessary.
@@ -965,7 +977,7 @@
   if (qulim) free(qulim);
   if (llim)  free(llim);
   if (ulim)  free(ulim);
-
+  if (dvecptr) free(dvecptr);
 
   return info;
 }
@@ -980,7 +992,8 @@
 	      double *wa, void *priv, int *nfev,
 	      double *step, double *dstep, int *dside,
 	      int *qulimited, double *ulimit,
-	      int *ddebug, double *ddrtol, double *ddatol)
+	      int *ddebug, double *ddrtol, double *ddatol,
+	      double *wa2, double **dvec)
 {
 /*
 *     **********
@@ -1063,17 +1076,15 @@
   int iflag = 0;
   double eps,h,temp;
   static double zero = 0.0;
-  double **dvec = 0;
   int has_analytical_deriv = 0, has_numerical_deriv = 0;
   int has_debug_deriv = 0;
   
   temp = mp_dmax1(epsfcn,MP_MACHEP0);
   eps = sqrt(temp);
   ij = 0;
-  ldfjac = 0; /* Prevents compiler warning */
+  ldfjac = 0;   /* Prevent compiler warning */
+  if (ldfjac){} /* Prevent compiler warning */
 
-  dvec = (double **) malloc(sizeof(double **)*npar);
-  if (dvec == 0) return MP_ERR_MEMORY;
   for (j=0; j<npar; j++) dvec[j] = 0;
 
   /* Initialize the Jacobian derivative matrix */
@@ -1164,12 +1175,12 @@
 		   (fjold == 0)?(0):((fjold-fjac[ij])/fjold));
 	  }
 	}
-      }
+      } /* end debugging */
 
-    } else {
+    } else {  /* dside > 2 */
       /* COMPUTE THE TWO-SIDED DERIVATIVE */
-      for (i=0; i<m; i++, ij++) {
-	fjac[ij] = wa[i];    /* Store temp data: fjac[i+m*j] */
+      for (i=0; i<m; i++) {
+	wa2[i] = wa[i];
       }
 
       /* Evaluate at x - h */
@@ -1180,15 +1191,16 @@
       x[ifree[j]] = temp;
 
       /* Now compute derivative as (f(x+h) - f(x-h))/(2h) */
-      ij -= m;
       if (! debug ) {
+	/* Non-debug path for speed */
 	for (i=0; i<m; i++, ij++) {
 	  fjac[ij] = (fjac[ij] - wa[i])/(2*h); /* fjac[i+m*j] */
 	}
       } else {
+	/* Debug path for correctness */
 	for (i=0; i<m; i++, ij++) {
 	  double fjold = fjac[ij];
-	  fjac[ij] = (fjac[ij] - wa[i])/(2*h); /* fjac[i+m*j] */
+	  fjac[ij] = (wa2[i] - wa[i])/(2*h); /* fjac[i+m*j] */
 	  if ((da == 0 && dr == 0 && (fjold != 0 || fjac[ij] != 0)) ||
 	      ((da != 0 || dr != 0) && (fabs(fjold-fjac[ij]) > da + fabs(fjold)*dr))) {
 	    printf("   %10d %10.4g %10.4g %10.4g %10.4g %10.4g\n", 
@@ -1196,17 +1208,16 @@
 		   (fjold == 0)?(0):((fjold-fjac[ij])/fjold));
 	  }
 	}
-      }	
+      } /* end debugging */
       
-    }
-  }
+    } /* if (dside > 2) */
+  } /* if (has_numerical_derivative) */
 
   if (has_debug_deriv) {
     printf("FJAC DEBUG END\n");
   }
 
  DONE:
-  if (dvec) free(dvec);
   if (iflag < 0) return iflag;
   return 0; 
   /*
@@ -1305,8 +1316,10 @@
   static double one = 1.0;
   static double p05 = 0.05;
 
-  lda = 0;   /* Prevent compiler warning */
-  lipvt = 0; /* Prevent compiler warning */
+  lda = 0;      /* Prevent compiler warning */
+  lipvt = 0;    /* Prevent compiler warning */
+  if (lda) {}   /* Prevent compiler warning */
+  if (lipvt) {} /* Prevent compiler warning */
 
   /*
    *     compute the initial column norms and initialize several arrays.