;+
; NAME:
;  vt3d_locavg
; PURPOSE: (one line)
;  Return spatial average over location
; DESCRIPTION:
;  Return spatial average over location
; CATEGORY:
;  Util
; CALLING SEQUENCE:
;  avg = vt3d_locavg(val, weight, indx=indx, default=default)
; INPUTS:
;  val = a scalar, array[n_loc], or array[n_loc,n_time]
;  weight = array float[n_loc] 
;  - usually angular area of surface elements (ster)
; OPTIONAL INPUT PARAMETER
;  indx = array long[n_indx], 0 <= indx < nloc.  
;    If not passed, all indices are used.
;    If indx = -1, use "default"
;  default = returned value if indx = -1 (e.g., do not choose any elements)
; OUTPUTS:
;  average of val or val[indx], weighted by area (weight)
; PROCEDURE
;  At its core this is just total(val * weight)/ total(weight)
;  This is very similar to wmean, except:
;  1. it accepts weights instead of errors
;  2. it can accpet the indx keyword
;  3. the names are taylored for weighting a map
;  4. If val is 2-D, then average over the first dimension only
;  Young, L. A. 2016, Volatile transport on inhomogeneous surfaces: II. Numerical calculations (VT3D)
;   Resubmitted to Icarus.
;   Eq. 5.3-3
; MODIFICATION HISTORY:
;  Written 2011 Aug 3, by Leslie Young, SwRI
;  2016 Apr 21, moved into layoung IDL library
;-

function vt3d_locavg, val, weight, indx=indx, default=default

n_loc = n_elements(weight)

; get dimensions of val
n_dim_val = size(val, /n_dim)
dim_val = size(val, /dim)

; val is a scalar
if n_dim_val eq 0 then begin
   return, val
endif

; val is an array, [nloc]
if n_dim_val eq 1 then begin
   ; all the values
   if not keyword_set(indx) then begin
      return, total(val * weight)/ total(weight)
   endif 
   ; none of the values
   if indexcount(indx) eq 0 then begin
      if keyword_set(default) then return, default else return, 0.
   endif
   ; some of the values
   return, total(val[indx] * weight[indx])/ total(weight[indx])
endif

; val is an matrix [nloc, nelse]
if n_dim_val eq 2 then begin
   nxnyn, val, nloc2, nelse
   ; all the values
   if not keyword_set(indx) then begin
      return, total(( val * (weight#onev(nelse)) ),1)/total(weight)
   endif 
   ; none of the values
   if indexcount(indx) eq 0 then begin
      if not keyword_set(default) then default else return = 0.
      return, replicate(default, nelse)
   endif
   ; some of the values
   return, total(( val[indx,*] * (weight[indx]#onev(nelse)) ),1)/total(weight[indx])
end

return, -1

end