;+
;  NAME: newton_raphson_3d.pro
;
;  PURPOSE:
;    A quick-and-dirty three-dimensional version of the iterative
;    Newton-Raphson method for finding a null point in a vector field.  For a
;    vector field Bvec(xvec), the main step is:
;      xvec_{n+1} = xvec_{n} - [ dB_i/dx_j | xvec_{n} ] ^ {-1} * Bvec(xvec)
;
;  CALLING SEQUENCE:
;    result = newton_raphson_3d(funcname,xstart,status=status)
;
;  INPUTS:
;    funcname = the name of a function that evaluates Bvec given a field point
;               xvec (syntax = Bvec = funcname(xvec,dmatrix=dmatrix), where
;               xvec and Bvec are 3-element arrays of type FLOAT or DOUBLE),
;               and dmatrix is a keyword that, on output, contains the matrix
;               of dB_i/dx_j values at xvec 
;    xstart = a 3-element array containing the starting point for the iteration
;
;  KEYWORDS
;    status = a status parameter that indicates problems with the minimization
;             process: 
;               0 = no problems encountered
;               1 = inversion attempted on a singular array (which indicates
;                   that the inversion is invalid) (returned from IDL's invert
;                   function)
;               2 = warning that a small pivot element was used during the
;                   matrix inversion step, and that significant accuracy was
;                   probably lost (returned from IDL's invert function)
;               3 = iteration step limit was reached before an acceptable
;                   minimum was found
;    itmax = iteration limit (default=10)
;    loud = if set, prints iteration number and value of x and B
;
;  OUTPUTS:
;    result = a 3-element array containing a root of Bvec, assuming everything
;             went well 
;
;  NOTES:
;    -The routine has not undergone extensive testing (so be careful!).
;
;  MODIFICATION HISTORY:
;    M.DeRosa - 22 Nov 2010 - created
;
;-

function newton_raphson_3d,funcname,xstart,status=status,itmax=itmax,loud=loud

;  usage message
if n_params() ne 2 then begin
  print,'  result = newton_raphson_3d(funcname,xstart,status=status)'
  return,-1
endif

;  idiot checking
if n_elements(xstart) ne 3 then begin
  print,'  ERROR in newton_raphson_3d: xstart must have 3 elements'
  return,-1
endif
if strlen(funcname) eq 0 then begin
  print,'  ERROR in newton_raphson_3d: funcname must not be null'
endif
if not keyword_set(itmax) then itmax=10

;  starting point
p0=xstart
bp0=call_function(funcname,p0,dmatrix=dbidxj)

;  iterate
iter=0
status=0
flag=0
if keyword_set(loud) then print,iter,[p0,bp0]
repeat begin

  ;  invert matrix
  dbidxj_inv=invert(dbidxj,status_invert)
  status=status>status_invert

  ;  compute next point
  p1=p0-dbidxj_inv#bp0

  ;  test new value
  bp1=call_function(funcname,p1,dmatrix=dbidxj)
  if sqrt(total(bp1*bp1)) lt 1e-5 then flag=1  ;  hit null
  if iter gt itmax then begin
    flag=2
    status=3
  endif

  ;  reset values of p0
  p0=p1
  bp0=bp1
  iter=iter+1
  if keyword_set(loud) then print,iter,[p0,bp0]

endrep until (flag gt 0)

;  return
return,p0

end