/* optm3.c */
#include <stdio.h>
#undef abs
#define abs(x) ((x)>=0.0)?(x):(-x)

void optm3(double *xf, double *yf, double *zf,
                                          /* initial and returned optimum */
           double *v,                     /* returned value */
           int    *cnt, int maxcnt,       /* number of function evaluations */
           double vtol,                   /* get v to this tolerance */
           double *step, double xyztol,   /* start and minimum xyz step size */
           double xmin, double xmax,      /* range for x */
           double ymin, double ymax,      /* range for y */
           double zmin, double zmax,      /* range for z */
           double (*f)(double x, double y, double z) )
{                                         /* function to be evaluated */
  double x, y, z;
  double dx, xl, xh, vxl, vxh;
  double dy, yl, yh, vyl, vyh;
  double dz, zl, zh, vzl, vzh;
  double v1, v2, v3;
  double dxv, dyv, dzv;
  int best;
  int debug = 1;

  while(1) /* coming in with  cnt==0  gets return each iteration */
  {
    best = 0;
    x = *xf;
    y = *yf;
    z = *zf;
    dx = (xmax-xmin)* *step; /* step size changes */
    dy = (ymax-ymin)* *step;
    dz = (zmax-zmin)* *step;
    v1 =f(x, y, z);
    v1 = abs(v1);
    v2 =v1;
    vxl=f(x-dx, y, z);
    vxl = abs(vxl);
    if(vxl < v2) { v2 = vxl; best = 1; }
    vxh=f(x+dx, y, z);
    vxh = abs(vxh);
    if(vxh < v2) { v2 = vxh; best = 2; }
    dxv = 0.0; /* no improvement moving x */
    if(vxl > v1 && vxh < v1) dxv = dx;
    if(vxl < v1 && vxh > v1) dxv = -dx;
    if(vxh < v1 && vxh < vxl) dxv = dx;
    if(vxl < v1 && vxl < vxh) dxv = -dx;
    vyl=f(x, y-dy, z);
    vyl = abs(vyl);
    if(vyl < v2) { v2 = vyl; best = 3; }
    vyh=f(x, y+dy, z);
    vyh = abs(vyh);
    if(vyh < v2) { v2 = vyh; best = 4; }
    dyv = 0.0; /* no improvement moving y */
    if(vyl > v1 && vyh < v1) dyv = dy;
    if(vyl < v1 && vyh > v1) dyv = -dy;
    if(vyh < v1 && vyh < vyl) dyv = dy;
    if(vyl < v1 && vyl < vyh) dyv = -dy;
    vzl=f(x, y, z-dz);
    vzl = abs(vzl);
    if(vzl < v2) { v2 = vzl; best = 5; }
    vzh=f(x, y, z+dz);
    vzh = abs(vzh);
    if(vzh < v2) { v2 = vzh; best = 6; }
    dzv = 0.0; /* no improvement moving z */
    if(vzl > v1 && vzh < v1) dzv = dz;
    if(vzl < v1 && vzh > v1) dzv = -dz;
    if(vzh < v1 && vzh < vzl) dzv = dz;
    if(vzl < v1 && vzl < vzh) dzv = -dz;
    v3 =f(x+dxv, y+dyv, z+dzv);
    v3 = abs(v3);
    *cnt += 7;
    if(v3 < v2) { v2 = v3; best = 7; }
    *v = v2;
    switch(best)
    {
      case 0: *step = *step*0.7071; break;
      case 1: *xf = x-dx; *yf = y; *zf = z; break;
      case 2: *xf = x+dx; *yf = y; *zf = z; break;
      case 3: *xf = x; *yf = y-dy; *zf = z; break;
      case 4: *xf = x; *yf = y+dy; *zf = z; break;
      case 5: *xf = x; *yf = y; *zf = z-dz; break;
      case 6: *xf = x; *yf = y; *zf = z+dz; break;
      case 7: *xf = x+dxv; *yf = y+dyv; *zf = z+dzv; break;
    } /* end this iteration, test for finish */
    if(debug) printf("best=%d, v=%g, step=%g, cnt=%d \n",
                     best, v2, *step, *cnt);
    if(*step < xyztol) return;
    if(vtol < v2) return;
    if(*cnt > maxcnt) return;
  } /* end while */
} /* end optm3.c */