Problems with MPFIT and user-defined derivatives

I am trying to use the optimization library MPFIT to fit a Gaussian function to my data. Actually this code is part of the example code that comes with the MPFIT library. The original code automatically calculates internally the derivatives of the function numerically and it works perfectly. The MPFIT library also allows the user to provide the function derivatives. This is where the problem starts. Here is the function used to calculate the residuals and the function's first order partial derivatives.

int gaussfunc(int m, int n, double *p, double *dy, double **derivs, void *vars)
{
  int i,j;
  struct vars_struct *v = (struct vars_struct *) vars;
  double *x, *y, *ey;
  double a = p[1];
  double b = p[2];
  double c = p[3];
  double d = p[0];
  x = v->x;
  y = v->y;
  ey = v->ey;

  for (i=0; i<m; i++)
  {
    dy[i] = (y[i] - (a*exp(-(x[i]-b)*(x[i]-b)/(2*c*c))+d))/ey[i];
  }

  // the code below this point is the code I added to calculate the derivatives.
  if(derivs)
  {
      for(j = 0; j < n; j++)
      {
          if (derivs[j])
          {
              for (i = 0; i < m; i++)
              {
                  double da =     exp(-(x[i]-b)*(x[i]-b)/(2*c*c));
                  double db = a * exp(-(x[i]-b)*(x[i]-b)/(2*c*c)) * (x[i]-b)/(c*c);
                  double dc = a * exp(-(x[i]-b)*(x[i]-b)/(2*c*c)) * (x[i]-b)*(x[i]-b)/(c*c*c);
                  double dd = 1;
                  double foo;
                  if     (j == 0) foo = dd;
                  else if(j == 1) foo = da;
                  else if(j == 2) foo = db;
                  else if(j == 3) foo = dc;
                  derivs[j][i] = foo;
              }
          }
      }
  }

  return 0;
}

The code above the line 'if (derivs)' is the original one but refactored and the code below that is my code for computing the derivatives. I believe my maths are correct, and they are verified by https://math.stackexchange.com/questions/716545/calculating-the-first-order-partial-derivatives-of-the-gaussian-function/716553

Has anyone encountered the same problem while using MPFIT with user-defined derivatives?

Thank you.

Because the calculation of the residuals is (DATA-MODEL)/SIGMA, the derivatives should be:

[-d(MODEL)/d(PARAM)]/sigma.

So, this line:

derivs[j][i] = foo;

becomes:

derivs[j][i] = -foo/ey[i];

Problem solved! Thanks!