// w=harten0(t,x) such that: w - g(w) = 0 where g(w) = sin(pi*(x-w*t))
struct harten0 {
  harten0 (Float t1) : t0(t1), pi(acos(Float(-1))), verbose(false) {}
  Float operator() (const point& x) const {
    Float x0 = fabs(x[0]);
    // find a good Newton's method starting point w0 vs w_star / g'(w_star) = 0
    Float wmin = 0, wmax = 1, w0 = 0.5;
    if (t0 > 0) {
      Float x_star  = (x0 < 0.5) ? x0+0.5 : x0-0.5;
      Float w_star  = x_star/t0;
      Float w_star2 = (x_star+1)/t0;
      Float g_star  = sin(pi*(x0-w_star*t0));
      if (w_star > wmax) {
        w0 = 0.5;
      } else {
        if (g_star < w_star) {
          w0 = (w_star + wmin)/2;
        } else {
          if (w_star2 > wmax) {
            w0 = (w_star + wmax)/2;
          } else {
            w0 = (w_star + w_star2)/2;
          }
        }
      }
    }
    // use a list of starting points for Newton
    Float wini_list[] = {w0, wmax, wmin};
    // loop
    Float r = 0;
    size_t i = 0;
    int status = 0;
    const Float tol = 1e2*numeric_limits<Float>::epsilon();
    Float w = 0;
    for (Float wini : wini_list) {
      if (verbose) derr << "[harten] wmin="<<wmin<<", wmax="<<wmax<<", wini="<<wini<<endl;
      if (verbose) derr << "[harten] i r dw w, for t0=" << t0 << ", x0="<<x0<< endl;
      w = wini;
      i = 0;
      status = 0;
      do {
        r = f(w,t0,x0);
        Float dw = -r/df_dw(w,t0,x0);
        if (verbose) derr << "[harten] " << i << " " << fabs(r)
                          << " " << dw << " " << w << endl;
        if (fabs(r) <= tol && fabs(dw) <= tol) { break; }
        if (++i >= max_iter) { status = 1; break; }
        if (w+dw > 0 && w+dw < wmax) {
          w += dw; 
        } else if (w+dw <= 0) {
          w = (w + wmin)/2;
          check_macro (1+w != w, "Newton iteration: machine precision problem.");
        } else {
          w = (w + wmax)/2;
        }
      } while (true);
      if (status == 0) break;
      if (verbose) derr << "[harten] failed: restart with another starting point..." << endl;
    }
    check_macro (status == 0, "t = " << t0 << ", x = " << x0 << " : Newton iteration " << i
	<< ": precision " << tol << " not reached: " << fabs(r));
    return (x[0] >= 0) ? w : -w;
  }
  void set_verbose() { verbose = true; }
  void set_noverbose() { verbose = false; }
protected:
  Float f     (Float w, Float t, Float x) const { return w -      sin(pi*(x-w*t)); }
  Float df_dw (Float w, Float t, Float x) const { return 1 + pi*t*cos(pi*(x-w*t)); }
  static const size_t max_iter = 100;
  Float t0, pi;
  bool verbose;
};
