;; Fu Jie Huang , 8/2002
;; a quick test of the  one-dimensional root finding 
;; example from <GSL reference manual> pp.441
;; find roots of function "f(x)=x^2-5" within bound [0, 5]
;; result: x=2.236

(libload "gsl/const.lsh")
(libload "gsl/aux_globalvar")
(libload "gsl/aux_structure")
(libload "gsl/roots")

(de quadratic ( x params)
    ((-double-) x)
    ((-gptr-) params)
    (let ((ret (double 0)))
      #{ {
      double * p = (double *) $params;
      double a = *p++;
      double b = *p++;
      double c = *p;
      $ret = (a * $x + b) * $x +c;
      } #} 
      ret))

(dhc-make () quadratic)

(setq pp (double-matrix 3))
(pp () '(1.0 0.0 -5.0))

(setq f  (alloc_gsl_function))
(setval_function_gsl_function f (gptr quadratic))
(setval_params_gsl_function f (idx-ptr pp))

(setq ss (gsl_root_fsolver_alloc (var_gsl_root_fsolver_brent)))

(gsl_root_fsolver_set  ss  f   0.0    5.0)

(setq iter 0)
(do-while (and (= status GSL_CONTINUE) (< iter 100))
  (incr iter)
  (setq status (gsl_root_fsolver_iterate ss))
  (setq r      (gsl_root_fsolver_root ss))
  (setq x_lo   (gsl_root_fsolver_x_lower ss))
  (setq x_hi   (gsl_root_fsolver_x_upper ss))
  (setq status (gsl_root_test_interval    x_lo x_hi  0   0.001))
  (when (= status GSL_SUCCESS) (printf "converged: %g\n" r))
)



