mpfit.pro

      endif
      if n_elements(xall) GT 0 AND n_elements(xall) NE n_elements(parinfo) $
        then begin
          errmsg = 'ERROR: number of elements in PARINFO and P must agree'
          goto, TERMINATE
      endif
  endif

  ;; If the parameters were not specified at the command line, then
  ;; extract them from PARINFO
  if n_elements(xall) EQ 0 then begin
      mpfit_parinfo, parinfo, tagnames, 'VALUE', xall, status=status
      if status EQ 0 then begin
          errmsg = 'ERROR: either P or PARINFO[*].VALUE must be supplied.'
          goto, TERMINATE
      endif

      sz = size(xall)
      ;; Convert to double if parameters are not float or double
      if sz[sz[0]+1] NE 4 AND sz[sz[0]+1] NE 5 then $
        xall = double(xall)
  endif
  xall = xall[*]   ;; Make sure the array is 1-D
  npar = n_elements(xall)
  zero = xall[0] * 0.
  one  = zero    + 1.
  fnorm  = -one
  fnorm1 = -one

  ;; TIED parameters?
  mpfit_parinfo, parinfo, tagnames, 'TIED', ptied, default='', n=npar
  ptied = strtrim(ptied, 2)
  wh = where(ptied NE '', qanytied) 
  qanytied = qanytied GT 0
  mpconfig = create_struct(mpconfig, 'QANYTIED', qanytied, 'PTIED', ptied)

  ;; FIXED parameters ?
  mpfit_parinfo, parinfo, tagnames, 'FIXED', pfixed, default=0, n=npar
  pfixed = pfixed EQ 1
  pfixed = pfixed OR (ptied NE '');; Tied parameters are also effectively fixed
  
  ;; Finite differencing step, absolute and relative, and sidedness of deriv.
  mpfit_parinfo, parinfo, tagnames, 'STEP',     step, default=zero, n=npar
  mpfit_parinfo, parinfo, tagnames, 'RELSTEP', dstep, default=zero, n=npar
  mpfit_parinfo, parinfo, tagnames, 'MPSIDE',  dside, default=0,    n=npar
  ;; Debugging parameters
  mpfit_parinfo, parinfo, tagnames, 'MPDERIV_DEBUG',  ddebug, default=0, n=npar
  mpfit_parinfo, parinfo, tagnames, 'MPDERIV_RELTOL', ddrtol, default=zero, n=npar
  mpfit_parinfo, parinfo, tagnames, 'MPDERIV_ABSTOL', ddatol, default=zero, n=npar

  ;; Maximum and minimum steps allowed to be taken in one iteration
  mpfit_parinfo, parinfo, tagnames, 'MPMAXSTEP', maxstep, default=zero, n=npar
  mpfit_parinfo, parinfo, tagnames, 'MPMINSTEP', minstep, default=zero, n=npar
  qmin = minstep *  0  ;; Remove minstep for now!!
  qmax = maxstep NE 0
  wh = where(qmin AND qmax AND maxstep LT minstep, ct)
  if ct GT 0 then begin
      errmsg = 'ERROR: MPMINSTEP is greater than MPMAXSTEP'
      goto, TERMINATE
  endif

  ;; Finish up the free parameters
  ifree = where(pfixed NE 1, nfree)
  if nfree EQ 0 then begin
      errmsg = 'ERROR: no free parameters'
      goto, TERMINATE
  endif

  ;; An external Jacobian must be checked against the number of
  ;; parameters
  if isext then begin
      if (nj/nv) NE nfree then begin
          errmsg = string(nv, nfree, nfree, $
           format=('("ERROR: EXTERNAL_FJAC must be a ",I0," x ",I0,' + $
                   '" array, where ",I0," is the number of free parameters")'))
          goto, TERMINATE
      endif
  endif

  ;; Compose only VARYING parameters
  xnew = xall      ;; xnew is the set of parameters to be returned
  x = xnew[ifree]  ;; x is the set of free parameters
  ; Same for min/max step diagnostics
  qmin = qmin[ifree]  & minstep = minstep[ifree]
  qmax = qmax[ifree]  & maxstep = maxstep[ifree]
  wh = where(qmin OR qmax, ct)
  qminmax = ct GT 0


  ;; LIMITED parameters ?
  mpfit_parinfo, parinfo, tagnames, 'LIMITED', limited, status=st1
  mpfit_parinfo, parinfo, tagnames, 'LIMITS',  limits,  status=st2
  if st1 EQ 1 AND st2 EQ 1 then begin

      ;; Error checking on limits in parinfo
      wh = where((limited[0,*] AND xall LT limits[0,*]) OR $
                 (limited[1,*] AND xall GT limits[1,*]), ct)
      if ct GT 0 then begin
          errmsg = 'ERROR: parameters are not within PARINFO limits'
          goto, TERMINATE
      endif
      wh = where(limited[0,*] AND limited[1,*] AND $
                 limits[0,*] GE limits[1,*] AND $
                 pfixed EQ 0, ct)
      if ct GT 0 then begin
          errmsg = 'ERROR: PARINFO parameter limits are not consistent'
          goto, TERMINATE
      endif
      

      ;; Transfer structure values to local variables
      qulim = limited[1, ifree]
      ulim  = limits [1, ifree]
      qllim = limited[0, ifree]
      llim  = limits [0, ifree]

      wh = where(qulim OR qllim, ct)
      if ct GT 0 then qanylim = 1 else qanylim = 0

  endif else begin

      ;; Fill in local variables with dummy values
      qulim = lonarr(nfree)
      ulim  = x * 0.
      qllim = qulim
      llim  = x * 0.
      qanylim = 0

  endelse

  ;; Initialize the number of parameters pegged at a hard limit value
  wh = where((qulim AND (x EQ ulim)) OR (qllim AND (x EQ llim)), npegged)

  n = n_elements(x)
  if n_elements(maxiter) EQ 0 then maxiter = 200L

  ;; Check input parameters for errors
  if (n LE 0) OR (ftol LE 0) OR (xtol LE 0) OR (gtol LE 0) $
    OR (maxiter LT 0) OR (factor LE 0) then begin
      errmsg = 'ERROR: input keywords are inconsistent'
      goto, TERMINATE
  endif

  if keyword_set(rescale) then begin
      errmsg = 'ERROR: DIAG parameter scales are inconsistent'
      if n_elements(diag) LT n then goto, TERMINATE
      wh = where(diag LE 0, ct)
      if ct GT 0 then goto, TERMINATE
      errmsg = ''
  endif

  if n_elements(state) NE 0 AND NOT keyword_set(extinit) then begin
      szst = size(state)
      if szst[szst[0]+1] NE 8  then begin
          errmsg = 'EXTERNAL_STATE keyword was not preserved'
          status = 0
          goto, TERMINATE
      endif
      if nfree NE n_elements(state.ifree) then begin
          BAD_IFREE:
          errmsg = 'Number of free parameters must not change from one '+$
            'external iteration to the next'
          status = 0
          goto, TERMINATE
      endif
      wh = where(ifree NE state.ifree, ct)
      if ct GT 0 then goto, BAD_IFREE

      tnames = tag_names(state)
      for i = 0L, n_elements(tnames)-1 do begin
          dummy = execute(tnames[i]+' = state.'+tnames[i])
      endfor
      wa4 = reform(efvec, n_elements(efvec))

      goto, RESUME_FIT
  endif

  common mpfit_error, mperr

  if NOT isext then begin
      mperr = 0
      catch_msg = 'calling '+fcn
      fvec = mpfit_call(fcn, xnew, _EXTRA=fcnargs)
      iflag = mperr
      if iflag LT 0 then begin
          errmsg = 'ERROR: first call to "'+fcn+'" failed'
          goto, TERMINATE
      endif
  endif else begin
      fvec = reform(efvec, n_elements(efvec))
  endelse

  catch_msg = 'calling MPFIT_SETMACHAR'
  sz = size(fvec[0])
  isdouble = (sz[sz[0]+1] EQ 5)
  
  mpfit_setmachar, double=isdouble

  common mpfit_profile, profvals
;  prof_start = systime(1)

  MACHEP0 = machvals.machep
  DWARF   = machvals.minnum

  szx = size(x)
  ;; The parameters and the squared deviations should have the same
  ;; type.  Otherwise the MACHAR-based evaluation will fail.
  catch_msg = 'checking parameter data'
  tp = szx[szx[0]+1]
  if tp NE 4 AND tp NE 5 then begin
      if NOT keyword_set(quiet) then begin
          message, 'WARNING: input parameters must be at least FLOAT', /info
          message, '         (converting parameters to FLOAT)', /info
      endif
      x = float(x)
      xnew = float(x)
      szx = size(x)
  endif
  if isdouble AND tp NE 5 then begin
      if NOT keyword_set(quiet) then begin
          message, 'WARNING: data is DOUBLE but parameters are FLOAT', /info
          message, '         (converting parameters to DOUBLE)', /info
      endif
      x = double(x)
      xnew = double(xnew)
  endif

  m = n_elements(fvec)
  if (m LT n) then begin
      errmsg = 'ERROR: number of parameters must not exceed data'
      goto, TERMINATE
  endif

  fnorm = mpfit_enorm(fvec)

  ;; Initialize Levelberg-Marquardt parameter and iteration counter

  par = zero
  iter = 1L
  qtf = x * 0.

  ;; Beginning of the outer loop
  
  OUTER_LOOP:

  ;; If requested, call fcn to enable printing of iterates
  xnew[ifree] = x
  if qanytied then mpfit_tie, xnew, ptied
  dof = (n_elements(fvec) - nfree) > 1L

  if nprint GT 0 AND iterproc NE '' then begin
      catch_msg = 'calling '+iterproc
      iflag = 0L
      if (iter-1) MOD nprint EQ 0 then begin
          mperr = 0
          xnew0 = xnew

          call_procedure, iterproc, fcn, xnew, iter, fnorm^2, $
            FUNCTARGS=fcnargs, parinfo=parinfo, quiet=quiet, $
            dof=dof, _EXTRA=iterargs
          iflag = mperr

          ;; Check for user termination
          if iflag LT 0 then begin  
              errmsg = 'WARNING: premature termination by "'+iterproc+'"'
              goto, TERMINATE
          endif

          ;; If parameters were changed (grrr..) then re-tie
          if max(abs(xnew0-xnew)) GT 0 then begin
              if qanytied then mpfit_tie, xnew, ptied
              x = xnew[ifree]
          endif

      endif
  endif

  ;; Calculate the jacobian matrix
  iflag = 2
  if NOT isext then begin
      catch_msg = 'calling MPFIT_FDJAC2'
      ;; NOTE!  If you change this call then change the one during
      ;; clean-up as well!
      fjac = mpfit_fdjac2(fcn, x, fvec, step, qulim, ulim, dside, $
                          iflag=iflag, epsfcn=epsfcn, $
                          autoderiv=autoderiv, dstep=dstep, $
                          FUNCTARGS=fcnargs, ifree=ifree, xall=xnew, $
                          deriv_debug=ddebug, deriv_reltol=ddrtol, deriv_abstol=ddatol)
      if iflag LT 0 then begin
          errmsg = 'WARNING: premature termination by FDJAC2'
          goto, TERMINATE
      endif
  endif else begin
      fjac = reform(efjac,n_elements(fvec),npar, /overwrite)
  endelse

  ;; Rescale the residuals and gradient, for use with "alternative"
  ;; statistics such as the Cash statistic.
  catch_msg = 'prescaling residuals and gradient'
  if n_elements(scalfcn) GT 0 then begin
      call_procedure, strtrim(scalfcn[0],2), fvec, fjac
  endif

  ;; Determine if any of the parameters are pegged at the limits
  npegged = 0L
  if qanylim then begin
      catch_msg = 'zeroing derivatives of pegged parameters'
      whlpeg = where(qllim AND (x EQ llim), nlpeg)
      whupeg = where(qulim AND (x EQ ulim), nupeg)
      npegged = nlpeg + nupeg
      
      ;; See if any "pegged" values should keep their derivatives
      if (nlpeg GT 0) then begin
          ;; Total derivative of sum wrt lower pegged parameters
          ;;   Note: total(fvec*fjac[*,i]) is d(CHI^2)/dX[i]
          for i = 0L, nlpeg-1 do begin
              sum = total(fvec * fjac[*,whlpeg[i]])
              if sum GT 0 then fjac[*,whlpeg[i]] = 0
          endfor
      endif
      if (nupeg GT 0) then begin
          ;; Total derivative of sum wrt upper pegged parameters
          for i = 0L, nupeg-1 do begin
              sum = total(fvec * fjac[*,whupeg[i]])
              if sum LT 0 then fjac[*,whupeg[i]] = 0
          endfor
      endif
  endif

  ;; Save a copy of the Jacobian if the user requests it...
  if keyword_set(calc_fjac) then output_fjac = fjac

  ;; =====================
  ;; Compute the QR factorization of the jacobian
  catch_msg = 'calling MPFIT_QRFAC'
  ;;  IN:      Jacobian        
  ;; OUT:      Hh Vects  Permutation  RDIAG  ACNORM
  mpfit_qrfac, fjac,     ipvt,        wa1,   wa2, /pivot
  ;; Jacobian - jacobian matrix computed by mpfit_fdjac2
  ;; Hh vects - house holder vectors from QR factorization & R matrix
  ;; Permutation - permutation vector for pivoting
  ;; RDIAG - diagonal elements of R matrix
  ;; ACNORM - norms of input Jacobian matrix before factoring

  ;; =====================
  ;; On the first iteration if "diag" is unspecified, scale
  ;; according to the norms of the columns of the initial jacobian
  catch_msg = 'rescaling diagonal elements'
  if (iter EQ 1) then begin
      ;; Input: WA2 = root sum of squares of original Jacobian matrix
      ;;        DIAG = user-requested diagonal (not documented!)
      ;;        FACTOR = user-requested norm factor (not documented!)
      ;; Output: DIAG = Diagonal scaling values
      ;;         XNORM = sum of squared scaled residuals
      ;;         DELTA = rescaled XNORM

      if NOT keyword_set(rescale) OR (n_elements(diag) LT n) then begin
          diag = wa2 ;; Calculated from original Jacobian
          wh = where (diag EQ 0, ct) ;; Handle zero values
          if ct GT 0 then diag[wh] = one
      endif
      
      ;; On the first iteration, calculate the norm of the scaled x
      ;; and initialize the step bound delta 
      wa3 = diag * x           ;; WA3 is temp variable
      xnorm = mpfit_enorm(wa3)
      delta = factor*xnorm
      if delta EQ zero then delta = zero + factor
  endif

  ;; Form (q transpose)*fvec and store the first n components in qtf
  catch_msg = 'forming (q transpose)*fvec'
  wa4 = fvec
  for j=0L, n-1 do begin
      lj = ipvt[j]
      temp3 = fjac[j,lj]
      if temp3 NE 0 then begin
          fj = fjac[j:*,lj]
          wj = wa4[j:*]
          ;; *** optimization wa4(j:*)
          wa4[j] = wj - fj * total(fj*wj) / temp3  
      endif
      fjac[j,lj] = wa1[j]
      qtf[j] = wa4[j]
  endfor
  ;; From this point on, only the square matrix, consisting of the
  ;; triangle of R, is needed.
  fjac = fjac[0:n-1, 0:n-1]
  fjac = reform(fjac, n, n, /overwrite)
  fjac = fjac[*, ipvt]                    ;; Convert to permuted order
  fjac = reform(fjac, n, n, /overwrite)

  ;; Check for overflow.  This should be a cheap test here since FJAC
  ;; has been reduced to a (small) square matrix, and the test is
  ;; O(N^2).
  wh = where(finite(fjac) EQ 0, ct)
  if ct GT 0 then goto, FAIL_OVERFLOW

  ;; Compute the norm of the scaled gradient
  catch_msg = 'computing the scaled gradient'
  gnorm = zero
  if fnorm NE 0 then begin
      for j=0L, n-1 do begin
          l = ipvt[j]
          if wa2[l] NE 0 then begin
              sum = total(fjac[0:j,j]*qtf[0:j])/fnorm
              gnorm = max([gnorm,abs(sum/wa2[l])])
          endif
      endfor
  endif

  ;; Test for convergence of the gradient norm
  if gnorm LE gtol then info = 4
  if info NE 0 then goto, TERMINATE
  if maxiter EQ 0 then begin
     info = 5
     goto, TERMINATE
  endif

  ;; Rescale if necessary
  if NOT keyword_set(rescale) then $
    diag = diag > wa2

  ;; Beginning of the inner loop
  INNER_LOOP:
  
  ;; Determine the levenberg-marquardt parameter
  catch_msg = 'calculating LM parameter (MPFIT_LMPAR)'
  par = mpfit_lmpar(fjac, ipvt, diag, qtf, delta, wa1, wa2, par=par)

  ;; Store the direction p and x+p. Calculate the norm of p
  wa1 = -wa1

  if qanylim EQ 0 AND qminmax EQ 0 then begin
      ;; No parameter limits, so just move to new position WA2
      alpha = one
      wa2 = x + wa1

  endif else begin
      
      ;; Respect the limits.  If a step were to go out of bounds, then
      ;; we should take a step in the same direction but shorter distance.
      ;; The step should take us right to the limit in that case.
      alpha = one

      if qanylim EQ 1 then begin
          ;; Do not allow any steps out of bounds
          catch_msg = 'checking for a step out of bounds'
          if nlpeg GT 0 then wa1[whlpeg] = wa1[whlpeg] > 0
          if nupeg GT 0 then wa1[whupeg] = wa1[whupeg] < 0

          dwa1 = abs(wa1) GT MACHEP0
          whl = where(dwa1 AND qllim AND (x + wa1 LT llim), lct)
          if lct GT 0 then $
            alpha = min([alpha, (llim[whl]-x[whl])/wa1[whl]])
          whu = where(dwa1 AND qulim AND (x + wa1 GT ulim), uct)
          if uct GT 0 then $
            alpha = min([alpha, (ulim[whu]-x[whu])/wa1[whu]])
      endif

      ;; Obey any max step values.

      if qminmax EQ 1 then begin
          nwa1 = wa1 * alpha
          whmax = where(qmax AND maxstep GT 0, ct)
          if ct GT 0 then begin
              mrat = max(abs(nwa1[whmax])/abs(maxstep[whmax]))
              if mrat GT 1 then alpha = alpha / mrat
          endif
      endif          

      ;; Scale the resulting vector
      wa1 = wa1 * alpha
      wa2 = x + wa1

      ;; Adjust the final output values.  If the step put us exactly
      ;; on a boundary, make sure we peg it there.
      sgnu = (ulim GE 0)*2d - 1d
      sgnl = (llim GE 0)*2d - 1d

      ;; Handles case of 
      ;;      ... nonzero *LIM ...     ... zero *LIM ...
      ulim1 = ulim*(1-sgnu*MACHEP0) - (ulim EQ 0)*MACHEP0
      llim1 = llim*(1+sgnl*MACHEP0) + (llim EQ 0)*MACHEP0

      wh = where(qulim AND (wa2 GE ulim1), ct)
      if ct GT 0 then wa2[wh] = ulim[wh]

      wh = where(qllim AND (wa2 LE llim1), ct)
      if ct GT 0 then wa2[wh] = llim[wh]
  endelse

  wa3 = diag * wa1
  pnorm = mpfit_enorm(wa3)

  ;; On the first iteration, adjust the initial step bound
  if iter EQ 1 then delta = min([delta,pnorm])

  xnew[ifree] = wa2
  if isext then goto, SAVE_STATE

  ;; Evaluate the function at x+p and calculate its norm
  mperr = 0
  catch_msg = 'calling '+fcn
  wa4 = mpfit_call(fcn, xnew, _EXTRA=fcnargs)
  iflag = mperr
  if iflag LT 0 then begin
      errmsg = 'WARNING: premature termination by "'+fcn+'"'
      goto, TERMINATE
  endif
  RESUME_FIT:
  fnorm1 = mpfit_enorm(wa4)
  
  ;; Compute the scaled actual reduction
  catch_msg = 'computing convergence criteria'
  actred = -one
  if 0.1D * fnorm1 LT fnorm then actred = - (fnorm1/fnorm)^2 + 1.

  ;; Compute the scaled predicted reduction and the scaled directional
  ;; derivative
  for j = 0L, n-1 do begin
      wa3[j] = 0
      wa3[0:j] = wa3[0:j] + fjac[0:j,j]*wa1[ipvt[j]]
  endfor

  ;; Remember, alpha is the fraction of the full LM step actually
  ;; taken
  temp1 = mpfit_enorm(alpha*wa3)/fnorm
  temp2 = (sqrt(alpha*par)*pnorm)/fnorm
  half  = zero + 0.5
  prered = temp1*temp1 + (temp2*temp2)/half
  dirder = -(temp1*temp1 + temp2*temp2)

  ;; Compute the ratio of the actual to the predicted reduction.
  ratio = zero
  tenth = zero + 0.1
  if prered NE 0 then ratio = actred/prered

  ;; Update the step bound
  if ratio LE 0.25D then begin
      if actred GE 0 then temp = half $
      else temp = half*dirder/(dirder + half*actred)
      if ((0.1D*fnorm1) GE fnorm) OR (temp LT 0.1D) then temp = tenth
      delta = temp*min([delta,pnorm/tenth])
      par = par/temp
  endif else begin
      if (par EQ 0) OR (ratio GE 0.75) then begin
          delta = pnorm/half
          par = half*par
      endif
  endelse

  ;; Test for successful iteration
  if ratio GE 0.0001 then begin
      ;; Successful iteration.  Update x, fvec, and their norms
      x = wa2
      wa2 = diag * x

      fvec = wa4
      xnorm = mpfit_enorm(wa2)
      fnorm = fnorm1
      iter = iter + 1
  endif

  ;; Tests for convergence
  if (abs(actred) LE ftol) AND (prered LE ftol) $
    AND  (0.5D * ratio LE 1) then info = 1
  if delta LE xtol*xnorm then info = 2
  if (abs(actred) LE ftol) AND (prered LE ftol) $
    AND (0.5D * ratio LE 1) AND (info EQ 2) then info = 3
  if info NE 0 then goto, TERMINATE

  ;; Tests for termination and stringent tolerances
  if iter GE maxiter then info = 5
  if (abs(actred) LE MACHEP0) AND (prered LE MACHEP0) $
    AND (0.5*ratio LE 1) then info = 6
  if delta LE MACHEP0*xnorm then info = 7
  if gnorm LE MACHEP0 then info = 8
  if info NE 0 then goto, TERMINATE

  ;; End of inner loop. Repeat if iteration unsuccessful
  if ratio LT 0.0001 then begin
      goto, INNER_LOOP
  endif

  ;; Check for over/underflow
  wh = where(finite(wa1) EQ 0 OR finite(wa2) EQ 0 OR finite(x) EQ 0, ct)
  if ct GT 0 OR finite(ratio) EQ 0 then begin
      FAIL_OVERFLOW:
      errmsg = ('ERROR: parameter or function value(s) have become '+$
                'infinite; check model function for over- '+$
                'and underflow')
      info = -16
      goto, TERMINATE
  endif

  ;; End of outer loop.
  goto, OUTER_LOOP

TERMINATE:
  catch_msg = 'in the termination phase'
  ;; Termination, either normal or user imposed.
  if iflag LT 0 then info = iflag
  iflag = 0
  if n_elements(xnew) EQ 0 then goto, FINAL_RETURN
  if nfree EQ 0 then xnew = xall else xnew[ifree] = x
  if n_elements(qanytied) GT 0 then if qanytied then mpfit_tie, xnew, ptied
  dof = n_elements(fvec) - nfree


  ;; Call the ITERPROC at the end of the fit, if the fit status is
  ;; okay.  Don't call it if the fit failed for some reason.
  if info GT 0 then begin
      
      mperr = 0
      xnew0 = xnew
      
      call_procedure, iterproc, fcn, xnew, iter, fnorm^2, $
        FUNCTARGS=fcnargs, parinfo=parinfo, quiet=quiet, $
        dof=dof, _EXTRA=iterargs
      iflag = mperr

      if iflag LT 0 then begin  
          errmsg = 'WARNING: premature termination by "'+iterproc+'"'
      endif else begin
          ;; If parameters were changed (grrr..) then re-tie
          if max(abs(xnew0-xnew)) GT 0 then begin
              if qanytied then mpfit_tie, xnew, ptied
              x = xnew[ifree]
          endif
      endelse

  endif

  ;; Initialize the number of parameters pegged at a hard limit value
  npegged = 0L
  if n_elements(qanylim) GT 0 then if qanylim then begin
      wh = where((qulim AND (x EQ ulim)) OR $
                 (qllim AND (x EQ llim)), npegged)
  endif

  ;; Calculate final function value (FNORM) and residuals (FVEC)
  if isext EQ 0 AND nprint GT 0 AND info GT 0 then begin
      catch_msg = 'calling '+fcn
      fvec = mpfit_call(fcn, xnew, _EXTRA=fcnargs)
      catch_msg = 'in the termination phase'
      fnorm = mpfit_enorm(fvec)
  endif

  if n_elements(fnorm) GT 0 AND n_elements(fnorm1) GT 0 then begin
      fnorm = max([fnorm, fnorm1])
      fnorm = fnorm^2.
  endif

  covar = !values.d_nan
  ;; (very carefully) set the covariance matrix COVAR
  if info GT 0 AND NOT keyword_set(nocovar) $
    AND n_elements(n) GT 0 $
    AND n_elements(fjac) GT 0 AND n_elements(ipvt) GT 0 then begin
      sz = size(fjac)
      if n GT 0 AND sz[0] GT 1 AND sz[1] GE n AND sz[2] GE n $
        AND n_elements(ipvt) GE n then begin
          catch_msg = 'computing the covariance matrix'
          if n EQ 1 then $
            cv = mpfit_covar(reform([fjac[0,0]],1,1), ipvt[0]) $
          else $
            cv = mpfit_covar(fjac[0:n-1,0:n-1], ipvt[0:n-1])
          cv = reform(cv, n, n, /overwrite)
          nn = n_elements(xall)
          
          ;; Fill in actual covariance matrix, accounting for fixed
          ;; parameters.
          covar = replicate(zero, nn, nn)
          for i = 0L, n-1 do begin
              covar[ifree, ifree[i]] = cv[*,i]
          end
          
          ;; Compute errors in parameters
          catch_msg = 'computing parameter errors'
          i = lindgen(nn)
          perror = replicate(abs(covar[0])*0., nn)
          wh = where(covar[i,i] GE 0, ct)
          if ct GT 0 then $
            perror[wh] = sqrt(covar[wh, wh])
      endif
  endif

;  catch_msg = 'returning the result'
;  profvals.mpfit = profvals.mpfit + (systime(1) - prof_start)

  FINAL_RETURN:
  mpfit_fencepost_active = 0
  nfev = mpconfig.nfev
  if n_elements(xnew) EQ 0 then return, !values.d_nan
  return, xnew

  
  ;; ------------------------------------------------------------------
  ;; Alternate ending if the user supplies the function and gradients
  ;; externally
  ;; ------------------------------------------------------------------

  SAVE_STATE:

  catch_msg = 'saving MPFIT state'

  ;; Names of variables to save
  varlist = ['alpha', 'delta', 'diag', 'dwarf', 'factor', 'fnorm', $
             'fjac', 'gnorm', 'nfree', 'ifree', 'ipvt', 'iter', $
             'm', 'n', 'machvals', 'machep0', 'npegged', $
             'whlpeg', 'whupeg', 'nlpeg', 'nupeg', $
             'mpconfig', 'par', 'pnorm', 'qtf', $
             'wa1', 'wa2', 'wa3', 'xnorm', 'x', 'xnew']
  cmd = ''

  ;; Construct an expression that will save them
  for i = 0L, n_elements(varlist)-1 do begin
      ival = 0
      dummy = execute('ival = n_elements('+varlist[i]+')')
      if ival GT 0 then begin
          cmd = cmd + ',' + varlist[i]+':'+varlist[i]
      endif
  endfor
  cmd = 'state = create_struct({'+strmid(cmd,1)+'})'
  state = 0

  if execute(cmd) NE 1 then $
    message, 'ERROR: could not save MPFIT state'

  ;; Set STATUS keyword to prepare for next iteration, and reset init
  ;; so we do not init the next time
  info = 9
  extinit = 0

  return, xnew

end