;+
; NAME:
;   PHUNWRAP
;
; AUTHOR:
;   Craig B. Markwardt, NASA/GSFC Code 662, Greenbelt, MD 20770
;   craigm@lheamail.gsfc.nasa.gov
;   UPDATED VERSIONs can be found on my WEB PAGE: 
;      http://cow.physics.wisc.edu/~craigm/idl/idl.html
;
; PURPOSE:
;   Unwrap phase jumps to recover cycle counts
;
; MAJOR TOPICS:
;   Mathematics
;
; CALLING SEQUENCE:
;   CYCLES = PHUNWRAP(PHASE, TOLERANCE=, MAXVAL=)
;
; DESCRIPTION:
;
;   PHUNWRAP unwraps a sequence of phases to produce a new series of
;   cycle counts.  Phase jumps due to crossing over the PHASE=0
;   boundary are removed by adding an integral number of cycles.  The
;   algorithm is based on the MATLAB "unwrap" function.
;
;   NOTE: the unwrapping process can be ambiguous if there is a phase
;   jump of more than a half cycle in the series.  For example, if the
;   phase changes by ~0.5 cycles, it is not possible to distinguish
;   whether there wasa +0.5 cycle or -0.5 cycle jump.  The most
;   accurate unwrapping can be performed if the PHASE series is nearly
;   continuous and does not have rapid phase changes.
;
;   Users can select the tolerance used to determine the phase jump.
;   Users can also select the definition of "1 cycle" by changing
;   MAXVAL.  By default, MAXVAL is 2*!DPI, which correspondes to 1
;   cycle = 2*!DPI radians, but other values of 1 (cycle), or 360
;   (degrees) are possible.
;
; INPUTS:
;
;   PHASE - phase series to be unwrapped.  Values should range from 0
;           to MAXVAL.  The ordering of the series is important.
;
; RETURNS:
;
;   A new series, expressed in cycles, with cycle jumps larger than
;   TOLERANCE removed.
;
; OPTIONAL KEYWORDS:
;
;   TOLERANCE - phase jump tolerance.  If the phase from one sample to
;               the next changes by more than TOLERANCE, then a single
;               cycle jump is assumed to have occurred.
;               DEFAULT: 0.5*MAXVAL
;
;   MAXVAL - Maximum value for phase. Common values are: 2*!DPI
;            (radians; DEFAULT); 1 (cycle); 360 (degrees), but any
;            positive value may be used.
;
; EXAMPLE:
;
;  ;; Set up some fake data
;  x = dindgen(100)/10d
;  y = x/2
;  ph = y MOD 1.0    ;; Mock phases
;
;  cycles = phunwrap(ph, maxval=1)
;
; MODIFICATION HISTORY:
;   Written and documented, CM, July 2003
;   Handle the case of unsigned integer input, CM, Feb 2006
;
;  $Id: phunwrap.pro,v 1.3 2006/03/28 14:19:53 craigm Exp $
;
;-
; Copyright (C) 2003, 2006, Craig Markwardt
; This software is provided as is without any warranty whatsoever.
; Permission to use, copy, modify, and distribute modified or
; unmodified copies is granted, provided this copyright and disclaimer
; are included unchanged.
;-

function phunwrap, ph, tolerance=tol0, maxval=maxval0

  common phunwrap_common, idlver
  if n_elements(idlver) EQ 0 then begin
      idlver = !version.release
  endif


  if n_elements(maxval0) EQ 0 then maxval = 2d*!dpi else maxval = maxval0(0)
  if n_elements(tol0) EQ 0 then tol = 0.5*maxval else tol = tol0(0)*maxval

  if n_elements(ph) LT 2 then return, ph

  sz = size(ph)
  tp = sz(sz(0)+1)

  ;; First order difference 
  case tp of 
      12: dph = [0, long(ph)-long(ph(1:*))]
      13: dph = [0, long64(ph)-long64(ph(1:*))]
      15: dph = [0, long64(ph)-long64(ph(1:*))]
      else: dph = [0, ph - ph(1:*)]
  endcase
  
  p = maxval * (fix((dph GT tol) EQ 1) - fix((dph LT (-tol)) EQ 1))
  if idlver GT 5.25 then begin
      ;; Use built-in version if available
      r = total(p, /cumulative) 
  endif else begin
      ;; .. if not, then use the lame FOR loop
      r = p
      for i = 1L, n_elements(r)-1 do $
        r(i) = r(i) + r(i-1)
  endelse

  return, ph+r
end