function local_maxmin,vec,in_width, $
                  MAXIMA=maxima,COUNT=count,POINT_ORDER=point_order

   count = 0
   ; make sure width is positive
   width = abs(in_width)
   n = n_elements(vec)
   ; Compute the first derivative of the vector
   vecp = vec[1:n-1] - vec[0:n-2]
   np = n - 1
   ; Collapse the derivative to just +1, 0, or -1
   vecps = intarr(np)
   z = where(vecp gt 0.0,count)
   if count ne 0 then vecps[z]=1
   z = where(vecp lt 0.0,count)
   if count ne 0 then vecps[z]=-1

   ; Compute the derivative of just the sign vectors (vecps)
   vecpps = vecps[1:np-1] - vecps[0:np-2]

   ; Keep the appropriate extremum
   if maxima then begin
      z = where(vecpps lt 0,nidx)
   endif else begin
      z = where(vecpps gt 0,nidx)
   endelse

   ; Create an index vector with just the good points.
   if nidx eq 0 then begin
      if maxima then begin
         idx = where(vec eq max(vec))
      endif else begin
         idx = where(vec eq min(vec))
      endelse
   endif else begin
      idx = z+1
   endelse
   ; Sort the extrema (actually, the absolute value)
   sidx = reverse(sort(abs(vec[idx])))
   ; Scan down the list of extrema, start with the brightest and take out
   ;   all extrema within width of the position.  Any that are too close should
   ;   be removed from further consideration.
   nsidx=nidx
   if width ge 1 then begin
      i=0
      while i lt nsidx-1 do begin
         z=where(abs(idx[sidx[i+1:nsidx-1]]-idx[sidx[i]]) le width,count)
         if count ne 0 then begin
            idx[sidx[z+i+1]] = -1
            z=where(idx[sidx] ge 0, nsidx)
            if nsidx eq 0 then return,-1 ; this should never happen
            sidx=sidx[z]
         endif
         i = i + 1
      endwhile
   endif
   count = nsidx
   if keyword_set(point_order) then begin
      idx=idx[sidx]
      idx=idx[sort(idx)]
      return,idx
   endif else begin
      return,idx[sidx]
   endelse

end

pro haar_transform_verynew, data, range, dilat, haar_trf, haar_var 
  del_range = range(1) - range(0)
  translat = range
  n_translat = n_elements(translat)
  haar_trf = fltarr(n_translat)
  ran = [reverse(min(range)-(findgen(dilat/del_range)+1.)*del_range),range,$
         max(range)+(findgen(dilat/del_range)+1.)*del_range]
  dat = [replicate(data(0), dilat/del_range), data, $
         replicate(data(n_elements(data)-1), dilat/del_range)]
  for j=0, n_translat-1 do begin
    up_ind = where(ran gt translat(j) AND ran le translat(j) + dilat)
    down_ind = where(ran lt translat(j) AND ran ge translat(j) - dilat)
    haar_trf(j) = (total(-dat(up_ind)) + total(dat(down_ind))) / dilat
  ENDFOR
  haar_var = total(haar_trf^2)
END

FUNCTION findzi,altitude,sig,dilation=dilation
;altitude is an array of heights
;sig is the signal you want to transform (temperature, etc..)
;dilation is the smoothing value, default is 100
  IF n_elements(dilation) NE 1 THEN dilation = 100
  num_recs = n_elements(z)
  haar_transform_verynew, sig, altitude, dilation, haar_trf, haar_var
  alt_ind = local_maxmin(abs(haar_trf),0,/maxima)
  if alt_ind(0) ne -1 then BEGIN
;    k = n_elements(alt_ind) > 2
    za = altitude[alt_ind]
    hta = abs(haar_trf(alt_ind))
    vals = (sort(za))           ;[0:k-2]
    za = za[vals]
    hta = hta[vals]
    gradient_alt = za(min(where(hta eq max(hta))))
  ENDIF ELSE gradient_alt = !values.f_nan
  return,gradient_alt
END