yorick banner

Home

Manual

Packages

Global Index

Keywords

Quick Reference


#if 0
u = double(indgen(100));
v = double(indgen(105))(-,);
data = 20*exp(-4*log(2)*(((u-21.7)/4.1)^2 + ((v-60.33)/4.1)^2)) + 10*(random(numberof(u), numberof(v)) - 0.5);
gg_fit_2d_spike(data);
#endif

require, "yeti_gist_gui.i";

func gg_draw_box (x0, y0, x1, y1, color=, width=)
{
  plg, color=color, type=1, width=width, marks=0,
    [y0,y1,y1,y0,y0], [x0,x0,x1,x1,x0];
}

func gg_draw_target (xc, yc, dx, dy, color=, width=)
{
  x0 = (x1 = xc - 0.5*dx) - dx;
  x3 = (x2 = xc + 0.5*dx) + dx;
  y0 = (y1 = yc - 0.5*dy) - dy;
  y3 = (y2 = yc + 0.5*dy) + dy;
  //write, "target:",xc,yc," - ",plsys();
  pldj, color=color, type=1, width=width,
    [x0,x2,x1,x1,xc,xc,x1,x2], [yc,yc,y1,y2,y0,y2,y1,y1],  
    [x1,x3,x2,x2,xc,xc,x1,x2], [yc,yc,y1,y2,y1,y3,y2,y2];
}

func gg_draw_roi (xc, yc, dx, dy, size, color=, width=)
{
  w = (size - 1)/2;
  wx = dx*(w + 0.5);
  wy = dy*(w + 0.5);
  x1 = xc - wx;
  x2 = xc - 0.5*dx;
  x3 = xc + 0.5*dx;
  x4 = xc + wx;
  y1 = yc - wy;
  y2 = yc - 0.5*dy;
  y3 = yc + 0.5*dy;
  y4 = yc + wy;
  pldj, color=color, type=1, width=width,
    [x1,x4,x4,x1,x2,x3,x3,x2], [y1,y1,y4,y4,y2,y2,y3,y3],  
    [x4,x4,x1,x1,x3,x3,x2,x2], [y1,y4,y4,y1,y2,y3,y3,y2];
}

func gg_fit_2d_spike (z, x0, y0, x1, y1, background=, size=)
{
  local model, param;

  dims = dimsof(z);
  nx = dims(2);
  ny = dims(3);
  if (is_void(x0)) x0 = 1.0;
  if (is_void(y0)) y0 = 1.0;
  if (is_void(x1)) x1 = x0 + (nx - 1);
  if (is_void(y1)) y1 = y0 + (ny - 1);
  xmin = double(x0);
  xmax = double(x1);
  ymin = double(y0);
  ymax = double(y1);
  xstep = (xmax - xmin)/(nx - 1);
  ystep = (ymax - ymin)/(ny - 1);
  x = xmin + xstep*indgen(0:nx-1);
  y = ymin + ystep*indgen(0:ny-1);
  x0 -= 0.5*xstep;
  x1 += 0.5*xstep;
  y0 -= 0.5*ystep;
  y1 += 0.5*ystep;

  /* Create widgets. */
  bw = 77; // button width
  bh = 25; // button height
  pad = 2; // padding space
  bd = 2;  // border width

  //Button = gg_button(text=string(0), width=bw, height=bh, bd=bd);
  
  Top = gg_toplevel(bd=bd, relief="groove", padx=pad, pady=pad);
  Zoom = gg_button(text="zoom", width=bw, height=bh, bd=bd);
  Unzoom = gg_button(text="unzoom", width=bw, height=bh, bd=bd);
  Verbose = gg_button(text="verbose", width=bw, height=bh, bd=bd);
  Redraw = gg_button(text="redraw", width=bw, height=bh, bd=bd);
  Ok = gg_button(text="ok", width=bw, height=bh, bd=bd);
  Cancel = gg_button(text="cancel", width=bw, height=bh, bd=bd);
  Fit = gg_button(text="fit", width=bw, height=bh, bd=bd);

  Message = gg_label(text="this is the help...", width=7*bw, height=5*bh,
                     bd=1, relief=GG_SUNKEN, bg=GG_WHITE, fg=GG_BLUE,
                     font="helveticaB", fontsize=18, justify="CH");

  /* Widgets for "size" selection. */
  size_max = min(nx, ny);
  size_min = 3;
  if (is_void(size)) size = 11;
  if (numberof(size)!=1 || dimsof(size)(1)!=0 || structof((size += 0)) != long ||
      size%2!=1) error, "bad value for keyword SIZE";
  if (size < size_min) error, "value of SIZE keyword too small";
  SizeLabel = gg_label(text="size:", width=bw, height=bh,
                        bd=bd, relief=GG_FLAT, justify="RH");
  SizeValue = gg_label(text=swrite(format="%d", size), width=bw, height=bh,
                       bd=1, relief=GG_SUNKEN, justify="LH",
                       bg=GG_WHITE, fg=GG_BLACK);
  SizePlus = gg_button(text="+", width=12, height=12, bd=1);
  SizeMinus = gg_button(text="-", width=12, height=12, bd=1);

  /* Widgets for "background" selection. */
  Background = gg_label(text="background:", width=bw, height=bh,
                        bd=bd, relief="flat", justify="RH");
  BackgroundNone = gg_button(text="none",
                             relief=GG_RAISED, fg=GG_DEFAULT_FG,
                             width=bw, height=bh, bd=bd);
  BackgroundConstant = gg_button(text="constant",
                                 relief=GG_RAISED, fg=GG_DEFAULT_FG,
                                 width=bw, height=bh, bd=bd);
  BackgroundSlope = gg_button(text="slope",
                                 relief=GG_RAISED, fg=GG_DEFAULT_FG,
                                 width=bw, height=bh, bd=bd);
  if (is_void(background)) background = 0;
  if      (background == 0) selection = BackgroundNone;
  else if (background == 1) selection = BackgroundConstant;
  else if (background == 2) selection = BackgroundSlope;
  else error, "bad value for keyword BACKGROUND";
  h_set, selection, relief=GG_SUNKEN, fg=GG_RED;

  /* Arrange widgets in a grid. */
  gg_manage, Top, 1,1, Zoom;
  gg_manage, Top, 2,1, Unzoom;
  gg_manage, Top, 3,1, Verbose;
  gg_manage, Top, 6,1, Redraw;

  gg_manage, Top, 1,2, Background, anchor="e";
  gg_manage, Top, 2,2, BackgroundNone;
  gg_manage, Top, 3,2, BackgroundConstant;
  gg_manage, Top, 4,2, BackgroundSlope;

  gg_manage, Top, 1,3, SizeLabel, rowspan=2;
  gg_manage, Top, 2,3, SizeValue, rowspan=2;
  gg_manage, Top, 3,3, SizePlus, anchor="ws", padx=1, pady=1;
  gg_manage, Top, 3,4, SizeMinus, anchor="wn", padx=1, pady=1;
  
  gg_manage, Top, 6,3, Fit, rowspan=2;

  gg_manage, Top, 1,5, Message, colspan=6;

  gg_manage, Top, 1,6, Ok;
  gg_manage, Top, 6,6, Cancel;

  fma;
  pli, z, x0, y0, x1, y1;
  limits;
  gg_paint, Top;
  
  /* Event loop. */
  npaints = 0;
  do_zoom = 0;
  verbose = 0;
  point = xpoint = ypoint = 0;
  for (;;) {
    ms = mouse(-1, 0, "");
    if (is_void(ms)) continue;
    clicked = gg_find_clicked(Top, ms);
    if (is_void(clicked)) {
      if (do_zoom) {
        gg_zoom, ms;
      } else {
        lim = limits();
        if (long(ms(9))==1 &&
            (ms(1) - lim(1))*(ms(1) - lim(2)) < 0.0 &&
            (ms(2) - lim(3))*(ms(2) - lim(4)) < 0.0 &&
            (ms(1) - x0)*(ms(1) - x1) < 0.0 &&
            (ms(2) - y0)*(ms(2) - y1) < 0.0) {
          if (point) {
            /* Erase previous mark. */
            gg_draw_roi, x(xpoint), y(ypoint), xstep, ystep, size,
              color=GG_XOR;
          }
          xpoint = long(floor((ms(1) - xmin)/xstep + 1.5));
          ypoint = long(floor((ms(2) - ymin)/ystep + 1.5));
          gg_draw_roi, x(xpoint), y(ypoint), xstep, ystep, size,
            color=GG_XOR;
          point = 1;
          ++npaints;
        }
      }
    } else if (clicked == BackgroundNone) {
      gg_paint, h_set(BackgroundNone, relief=GG_SUNKEN, fg=GG_RED);
      gg_paint, h_set(BackgroundConstant, relief=GG_RAISED, fg=GG_DEFAULT_FG);
      gg_paint, h_set(BackgroundSlope, relief=GG_RAISED, fg=GG_DEFAULT_FG);
      background = 0;
      ++npaints;
    } else if (clicked == BackgroundConstant) {
      gg_paint, h_set(BackgroundNone, relief=GG_RAISED, fg=GG_DEFAULT_FG);
      gg_paint, h_set(BackgroundConstant, relief=GG_SUNKEN, fg=GG_RED);
      gg_paint, h_set(BackgroundSlope, relief=GG_RAISED, fg=GG_DEFAULT_FG);
      background = 1;
      ++npaints;
    } else if (clicked == BackgroundSlope) {
      gg_paint, h_set(BackgroundNone, relief=GG_RAISED, fg=GG_DEFAULT_FG);
      gg_paint, h_set(BackgroundConstant, relief=GG_RAISED, fg=GG_DEFAULT_FG);
      gg_paint, h_set(BackgroundSlope, relief=GG_SUNKEN, fg=GG_RED);
      background = 2;
      ++npaints;
    } else if (clicked == SizePlus) {
      if (size < size_max) {
        if (point) gg_draw_roi, x(xpoint), y(ypoint), xstep, ystep, size,
                     color=GG_XOR;
        size = min(size+2, size_max);
        gg_paint, h_set(SizeValue, text=swrite(format="%d", size));   
        gg_draw_roi, x(xpoint), y(ypoint), xstep, ystep, size,
          color=GG_XOR;
        ++npaints;
      }
    } else if (clicked == SizeMinus) {
      if (size > size_min) {
        if (point) gg_draw_roi, x(xpoint), y(ypoint), xstep, ystep, size,
                     color=GG_XOR;
        size = max(size-2, size_min);
        gg_paint, h_set(SizeValue, text=swrite(format="%d", size));
        gg_draw_roi, x(xpoint), y(ypoint), xstep, ystep, size,
          color=GG_XOR;
        ++npaints;
      }
    } else if (clicked == Fit) {
      if (! point) {
        gg_paint, h_set(Message, text="You must select a position first:\nunselect \"zoom\" and\nclick mouse near spike in image.");
      } else {
        /* Region of interest (ROI). */
        w = (size - 1)/2;
        xrng = max(1, xpoint-w):min(nx, xpoint+w);
        yrng = max(1, ypoint-w):min(ny, ypoint+w);
        zroi = z(xrng,yrng);
        xroi = x(xrng);
        yroi = y(yrng)(-,);
        param = [max(zroi), w, x(xpoint), y(ypoint)];
        if (background > 0) {
          grow, param, min(zroi);
          if (background == 2) grow, param, 0.0, 0.0;
        }
        param = fit_2d_spike(zroi, yroi, xroi, /*model, weight=, */
                             shape="spline", 
                             guess=param, bg=background,
                             absorption=0, verb=verbose);
        xpeak = param(3);
        ypeak = param(4);
        if (abs(x(xpoint) - xpeak) > 0.5*xstep ||
            abs(y(ypoint) - ypeak) > 0.5*ystep) {
          gg_draw_roi, x(xpoint), y(ypoint), xstep, ystep, size, 
            color=GG_XOR;
          xpoint = long(floor((xpeak - xmin)/xstep + 1.5));
          ypoint = long(floor((ypeak - ymin)/ystep + 1.5));
          gg_draw_roi, x(xpoint), y(ypoint), xstep, ystep, size,
            color=GG_XOR;
          msg = "*** WARNING *** box has moved: click fit again.\n";
        } else {
          msg = "";
        }
        fwhm = param(2);
        msg += swrite(format="peakmax = %g   fwhm = %g (%.1f%.1f)",
                     param(1), fwhm, fwhm/xstep, fwhm/ystep);
        msg += swrite(format="\nxpeak = %g   ypeak = %g", xpeak, ypeak);
        if (background > 0) {
          msg += swrite(format="\nbackground = %g", param(5));
          if (background == 2)
            msg += swrite(format="+ %g*(x - xpeak) + %g*(y - ypeak)",
                          param(6), param(7));
        }
        gg_paint, h_set(Message, text=msg);
      }
      ++npaints;
    } else if (clicked == Zoom) {
      do_zoom = ! do_zoom;
      if (do_zoom) h_set, clicked, relief=GG_SUNKEN, fg=GG_RED;
      else         h_set, clicked, relief=GG_RAISED, fg=GG_DEFAULT_FG;
      gg_paint, clicked;
      ++npaints;
    } else if (clicked == Verbose) {
      verbose = ! verbose;
      if (verbose) h_set, clicked, relief=GG_SUNKEN, fg=GG_RED;
      else         h_set, clicked, relief=GG_RAISED, fg=GG_DEFAULT_FG;
      gg_paint, clicked;
      ++npaints;
    } else if (clicked == Unzoom || clicked == Redraw) {
      fma;
      pli, z, x0, y0, x1, y1;
      if (clicked == Unzoom) limits;
      gg_paint, Top;
      if (point) gg_draw_roi, x(xpoint), y(ypoint), xstep, ystep, size,
                   color=GG_XOR;
      npaints = 0;
    } else if (clicked == Ok) {
      return param;
    } else if (clicked == Cancel) {
      return;
    }

    if (npaints > 30) {
      fma;
      pli, z, x0, y0, x1, y1;
      gg_paint, Top;
      if (point) gg_draw_roi, x(xpoint), y(ypoint), xstep, ystep, size,
                   color=GG_XOR;
      npaints = 0;
    }
  }
  
}

func __fit_spike_gauss (&p, &dp, u) {
  /* Gaussian shape: exp(-C1*u*u)
   * To have FWHM = 1, use:
   *   C1 = 4*log(2) ~ 2.77258872223978123766892848583
   *   C2 = 2*C1     ~ 5.54517744447956247533785697167
   */
  p = exp(-2.77258872223978123766892848583*u*u);
  dp = -5.54517744447956247533785697167*u*p;
}

/*
 * CUBIC B-SPLINE BASIS FUNCTION
 *   S(x) = (|x|/2 - 1)*x^2 + 2/3   if |x|<=1
 *        = 1/6*(2 - |x|)^3         if 1<=|x|<=2
 *        = 0                       otherwise
 *   S(x) is normalised (integral S(x) dx = 1) but we want a function
 *   F(x) = alpha*S(beta*x) so that:
 *     f(0) = 1
 *     f(0.5) = 0.5
 *     f(-u) = f(u)                              (i.e. f is symmetric)
 *     f(u) >= 0     everywhere
 *     f'(u)>0 for u<0, f'(0)=0, f'(u)<0 for u>0 (i.e. f is "bell-shaped")
 *   the 3 first properties are needed to properly define the FWHM (Full
 *   Half Width at Half Maximum).
 *
 *     f(0) = 1  <==>  alpha = 3/2
 *     f(x) = core(|x|)  if if |x|<=1
 *          = wing(|x|)  if 1<=|x|<=2
 *   with:
 *     core(x) = (beta*x*3/4 - 3/2)*(beta*x)^2 + 1
 *     wing(x) = (3/2)*(2 - beta*x)^3
 *   FWHM can be in the "core":
 *     core(1/2) = 1/2 <=> beta = (4 + 8*cos((arccos(-1/8) - 2*pi)/3))/3
 *                              ~ 1.444703448928752590331004298512891328576
 *   in this case, beta*1/2 <= 1, so its is possible.  FWHM can be in the
 *   "wing":
 *     wing(1/2) = 1/2 <=> beta = 4 - 2^(4/3)
 *                              ~ 1.48015790021025367046557878544354329886
 *   in this case, beta*1/2 < 1 so it is not possible.
 *
 *   Finally BETA = (4 + 8*cos((arccos(-1/8) - 2*pi)/3))/3
 *                ~ 1.444703448928752590331004298512891328576
 */
  
func __fit_spike_spline (&p, &dp, u) {
  /* To have FWHM = 1, use:
   * BETA = (4 + 8*cos((arccos(-1/8) - 2*pi)/3))/3
   *      ~ 1.444703448928752590331004298512891328576
   */
  BETA = 1.444703448928752590331004298512891328576;
  s = BETA*sign(u);
  u *= s;
  p = dp = array(double, dimsof(u));
  if (is_array((i = where(u <= 1.0)))) {
    ui = u(i);
    p(i) = (0.75*ui - 1.5)*ui*ui + 1.0;
    dp(i)= (2.25*ui - 3.0)*ui;
  }
  if (is_array((i = where((u > 1.0)*(u < 2.0))))) {
    ui = 2.0 - u(i);
    p(i) = 0.25*ui*ui*ui;
    dp(i) = -0.75*ui*ui;
  }
  dp *= s;
}

func __fit_lp (&lmin, &lmax, &lgain)
{
  if (is_void(lmin))  lmin = 1e-3;
  else if (lmin <= 0.0 || lmin > 1.0) error, "bad value for keyword LMIN";
  if (is_void(lmax))  lmax  = 1e30;
  else if (lmax <= 1.0) error, "bad value for keyword LMAX";
  if (is_void(lgain)) lgain = 10.0;
  else if (lgain <= 1.0) error, "bad value for keyword LGAIN";
  return double(lmin);
}

func fit_2d_spike (z, y, x, &model, weight=, shape=, 
                  guess=, bg=, absorption=, verb=,
                  maxeval=, maxiter=, conv=,
                  lmin=, lmax=, lgain=)
/* DOCUMENT

     Parameters:
       A(1) = peak intensity of spike
       A(2) = FWHM of spike
       A(3) = abscissa of spike
       A(4) = ordinate of spike
       A(5) = background level
       A(6) = background X-slope
       A(7) = background Y-slope
     Model:
       MODEL(X,Y) = A1*P((X-A3)/A2)*P((Y-A4)/A2) + A5 + A6*(X-A3) + A7*(Y-A4)
     where P(U) is the shape function.

     SHAPE is the 1-D shape of the spike, it can be "spline" to use a
     cubic B-spline function, "gauss" to use a gaussian or any user
     defined routine which has the following prototype:
       func SHAPE(&p, &dp, u) {gg_fit_2d_spike(data)
         p = f(u);
         dp = f'(u);
       }
  
     where f(u) is a "bell-shaped" function and f'(u) is
     its derivative; f(u) should have
     with the following properties:
       f(0) = 1
       f(0.5) = 0.5
       f(-u) = f(u)                              (i.e. f is symmetric)
       f(u) >= 0     everywhere
       f'(u)>0 for u<0, f'(0)=0, f'(u)<0 for u>0 (i.e. f is "bell-shaped")
     the 3 first properties are needed to properly define the FWHM (Full
     Half Width at Half Maximum).
       
     Background:
       BG = nil/0  no background (BG = 0 everywhere)
       BG = 1      uniform background
       BG = 2      BG(x,y) = a() + a()*(x-x0) + a()*(y-y0)
     Spike:
     Support: */
{
  local g, h, chi2_best, a_best, delta_chi2, delta_a, p1, dp1, p2, dp2;

  if (is_void(conv)) conv = 1e-5;
  
  /* Initializes seetings for Levenberg-Marquardt parameter (smallness of
     trust region). */
  lp = __fit_lp(lmin, lmax, lgain);

  if (! is_func(shape)) {
    if (is_void(shape) || shape == "spline") shape = __fit_spike_spline;
    else if (shape == "gauss") shape = __fit_spike_gauss;
    else error, "bad value for keyword SHAPE";
  }

  if (! bg) nparams = 4;
  else if (bg == 1) nparams = 5;
  else if (bg == 2) nparams = 7;
  else error, "bad value for keyword BG";

  /* ONE is used to make conformable arrays the same dimension list as Z. */
  one = array(1.0, dimsof(z));

  /* Check WEIGHT, compute W=sqrt(WEIGHT). */
  if (is_void(weight)) {
    weight = 0n;
  } else {
    s = structof(weight);
    if (s==double || s==float || s==int || s==long || s==short) {
      /* Numerical weight. */
      if (min(weight) < 0.0) error, "WEIGHT must be non-negative";
      if (max(weight) <= 0.0) error, "WEIGHT must not be zero everywhere";
      if (dimsof(weight)(1) != 0) {
        /* Non-scalar weight. */
        w = one*sqrt(weight);gg_fit_2d_spike(data)
        weight = 1n;
      } else {
        /* Scalar weight: same as with no weight. */
        weight = 0n;
      }
    } else error, "bad data type for WEIGHT";
  }

  if (is_void(guess)) {
    /* Starting solution. PEAK is the list of indices in the half topmost
       part of the spike (bottommost part for an absorption spike) and is
       used to estimate the spike position and FWHM. */
    a = array(double, nparams);
    zmin = min(z);
    zmax = max(z);
    if (absorption) {
      if (bg) {
        a(5) = zmax;
        a(1) = zmin - zmax;
        peak = where(z <= 0.5*(zmax + zmin));
      } else {
        a(1) = zmin;
        peak = where(z <= 0.5*zmin);
      }
    } else {
      if (bg) {
        a(5) = zmin;
        a(1) = zmax - zmin;
        peak = where(z >= 0.5*(zmax + zmin));
      } else {
        a(1) = zmax;
        peak = where(z >= 0.5*zmax);
      }
    }
    x_peak = (one*x)(peak);
    y_peak = (one*y)(peak);
    a(2) = 0.5*(abs(x_peak(ptp)) + abs(y_peak(ptp))); /* FWHM */
    a(3) = avg(x_peak);
    a(4) = avg(y_peak);
    peak = x_peak = y_peak = [];
  } else {
    s = structof(guess);
    if ((s==double || s==float || s==long || s==int || s==short || s==char) &&
        dimsof(guess)(1)==1 && numberof(guess) == nparams) {
      a = double(guess);
    } else {
      error, "bad GUESS parameters";
    }
  }


  /* Minimum allowed FWHM (taken to be a small fraction of the sampling
     to avoid division by zero). */
  fwhm_min = 1e-5*sqrt(((max(x)-min(x))*(max(y)-min(y)))/double(numberof(z)));
  //write, format="FWHM_min = %g\n", fwhm_min;
  
  /* Model:
   *   m = a1*p((x-a3)/a2)*p((y-a4)/a2) + a5 + a6*(x-a3) + a7*(y-a4)
   *
   * Jacobian:
   *   jac1 = p((x-a3)/a2)*p((y-a4)/a2)
   *   jac2 = -(a1/a2)*((x-a3)/a2*p'((x-a3)/a2)*p((y-a4)/a2) +
   *                    (y-a4)/a2*p((x-a3)/a2)*p'((y-a4)/a2))
   *   jac3 = -(a1/a2)*p((y-a4)/a2)*p'((x-a3)/a2) - a6
   *   jac4 = -(a1/a2)*p((x-a3)/a2)*p'((y-a4)/a2) - a7
   *   jac5 = 1
   *   jac6 = x-a3
   *   jac7 = y-a4
   *
   * Local quadratic approximation (l = likelihood, a = parameters,
   * g = gradient vector with respect to parameters, H = hessian matrix,
   * J = jacobian, i.e. paratial derivatives of model with repsect to
   * parameters):
   *   l(a) = sum_x weight(x) * (model(a,x) - data(x))^2
   *   g(a) = 2 * sum_x  weight(x) * (model(a,x) - data(x)) * J(a,x)
   *   H(a,b) ~ 2 * sum_x  weight(x) * J(a,x) * J(b,x)
   *   ==> l(a+da) ~ l(a) + da'.g(a) + 1/2 da'.H.da
   *   ==> da_best = -LUsolve(H,g)
   *               = -LUsolve(1/2 H, 1/2 g)
   * trust region control (lp = Levenberg-Marquardt parameter
   *                          = smallness of trust region):
   *   da(lp) = -LUsolve(1/2 H(lp), 1/2 g)
   * with:
   *   H(lp)(a,b) = H(a,b)          where a!=b
   *              = (1 + lp)*H(a,a) where a=b
   */
  
  warn = string(0);
  iter = -1;
  eval = 0; 
  for ( ; ; ) {
    /* Apply constraints. */
    if (a(2) < fwhm_min) a(2) = fwhm_min;
    
    /* Compute model M, (weighted) residuals R and Chi2. */
    s = 1.0/a(2); /* 1/FWHM */
    jac6 = x - a(3);
    jac7 = y - a(4);
    shape, p1, dp1, (u1 = s*jac6);
    shape, p2, dp2, (u2 = s*jac7);
    if (nparams == 4)      m = a(1)*p1*p2;
    else if (nparams == 5) m = a(1)*p1*p2 + a(5);
    else                   m = a(1)*p1*p2 + a(5) + a(6)*jac6 + a(7)*jac7;
    r = m - z;
    if (weight) r *= w;
    chi2 = sum(r*r);

    ++eval;
    if (! is_void(maxeval) && eval >= maxeval) {
      warn = swrite(format="too many evaluations (%d)", eval);
    }

    if (eval <= 1 || chi2 < chi2_best) {
      /* Solution improved: register current model and current solution as
         the best one found so far. */
      ++iter;
      eq_nocopy, model, m;
      if (iter) {
        delta_chi2 = abs(chi2_best - chi2);
        delta_a    = abs(a - a_best);
      }
      chi2_best = chi2;
      a_best = a;
      if (verb && iter%verb == 0) {
        if (eval==1) {
          write, format="%s\n%s\n",
            "ITER  EVAL  CHI2         FWHM       X0            Y0",
            "----  ----  -----------  ---------  ------------  ------------";
        }
        write, format="%4d %5d  %11.5e  %-9.3g  %-+12.5g  %-+12.5g\n",
          iter, eval, chi2, a(2), a(3), a(4);
      }
      
      /* Check for termination. */
      if (chi2 <= 0.0 || (iter && delta_chi2 <= conv*abs(chi2))) {
        warn = string(0);
        break;
      }
      if (! is_void(maxiter) && iter >= maxiter) {
        warn = swrite(format="too many iterations (%d)", iter);
        break;
      }
      if (warn) break; /* too many evaluations */
      
      /* Solution improved: augment trust region. */
      lp = max(lp/lgain, lmin);
      
      /* Compute Jacobian:
       *   jac1 = p1*p2
       *   jac2 = -(a1/a2)*(u1*dp1*p2 + u2*p1*dp2)
       *   jac3 = -(a1/a2)*p2*dp1 - a6
       *   jac4 = -(a1/a2)*p1*dp2 - a7
       *   jac5 = 1
       *   jac6 = x-a3
       *   jac7 = y-a4
       */
      jac = array(pointer, nparams);
      q = -a(1)/a(2);
      if (weight) {
        jac(1) = &(w*p1*p2);
        jac(2) = &(w*q*(u1*dp1*p2 + u2*p1*dp2));
        if (nparams >= 5) jac(5) = &w;
        if (nparams == 7) {
          jac(3) = &(w*(q*p2*dp1 - a(6)));
          jac(4) = &(w*(q*p1*dp2 - a(7)));
          jac(6) = &(w*jac6);
          jac(7) = &(w*jac7);
        } else {
          jac(3) = &(w*q*p2*dp1);
          jac(4) = &(w*q*p1*dp2);
        }
      } else {
        jac(1) = &(p1*p2);
        jac(2) = &(q*(u1*dp1*p2 + u2*p1*dp2));
        if (nparams >= 5) jac(5) = &one;
        if (nparams == 7) {
          jac(3) = &(q*p2*dp1 - a(6));
          jac(4) = &(q*p1*dp2 - a(7));
          jac(6) = &jac6;
          jac(7) = &jac7;
        } else {
          jac(3) = &(q*p2*dp1);
          jac(4) = &(q*p1*dp2);
        }
      }

      /* Compute (half of) Hessian and (half of) gradient to obtain local
         quadratic approximation. */
      h = array(double, nparams, nparams);
      g = array(double, nparams);
      for (i=1 ; i<=nparams ; ++i) {
        if (weight) {
          jac_i = w*(*jac(i));
        } else {
          local jac_i;
          eq_nocopy, jac_i, *jac(i);
        }
        for (j=1 ; j<i ; ++j) h(i,j) = (h(j,i) = sum(jac_i * *jac(j)));
        h(i,i) = sum(jac_i * jac_i);
        g(i) = sum(jac_i * r);
      }

      /* Account for bound constraints on the FWHM. */
      if (a(2) <= fwhm_min && g(2) > 0.0) g(2) = 0.0;
      
    } else {
      /* Solution not improved, reduce trust region. */
      if (lp >= lmax)
        warn = "Levenberg-Marquardt parameter overflow (convergence)";
      if (warn) break; /* parameter overflow or too many evaluations */
      lp = min(lp*lgain, lmax);
    }

    /* Solve the regularized local quadratic problem. */
    tmp = h;
    tmp(::nparams+1) *= 1.0 + lp;
    a = a_best - LUsolve(tmp, g);
  }
  if (warn && verb!=0) write, format="WARNING - %s\n", warn;
  return a;
}