undef("unwrap_phase")
function unwrap_phase(pha[*]:numeric)    ; pha <==> wrapped phase
;
; This should replicate the Matlab function: unwrap
;
;---This function is based on Matlab code at: 
;---https://www.ljmu.ac.uk/~/media/files/ljmu/about-us/faculties-and-schools/tae/geri/onedimensionalphaseunwrapping_finalpdf

local N, K, increments, i, ii, nii, iStrt, iLast, diff, pha2, pi, pi2
begin
  N    = dimsizes(pha)
  K    = totype( 0, typeof(pha))

  increments = new( N, typeof(pha), getVarFillValue(pha))
  increments = totype( 0, typeof(pha))

  if (any(ismissing(pha))) then
      ii    = ind(.not.ismissing(pha))
      nii   = dimsizes(ii)
      iStrt = ii(0)+1
      iLast = ii(nii-1)
  else
      iStrt = 1             ; NCL is 0-based
      iLast = N-1
  end if

  pi  = get_pi(pha)         ; NCL 6.4.0
  pi2 = 2*pi

  do i=iStrt,iLast
     diff = pha(i)-pha(i-1)
     if (diff.gt.pi) then
         K = K - pi2
     elseif (diff.lt.(-pi)) then
         K = K + pi2
     end if
     increments(i) = K
  end do

  pha2 = pha + increments

;;pha2(iStrt-1) = pha2(iStrt)     ; Matlab does this ... near as I can tell

  pha2@long_name    = "unwrapped phase"
  pha2@units        = "radians"
  pha2@NCL_function = "unwrap_phase"
  return(pha2)
end