FUNCTION ngp,value,posx,nx,posy,ny,posz,nz, $
             AVERAGE=average,WRAPAROUND=wraparound,NO_MESSAGE=no_message
;+
; NAME:
;       NGP
;
; PURPOSE:
;       Interpolate an irregularly sampled field using Nearest Grid Point
;
; EXPLANATION:
;       This function interpolates irregularly gridded points to a
;       regular grid using Nearest Grid Point.
;
; CATEGORY:
;       Mathematical functions, Interpolation
;
; CALLING SEQUENCE:
;       Result = NGP, VALUE, POSX, NX[, POSY, NY, POSZ, NZ, 
;                     /AVERAGE, /WRAPAROUND, /NO_MESSAGE]
;
; INPUTS:
;       VALUE: Array of sample weights (field values). For e.g. a
;              temperature field this would be the temperature and the
;              keyword AVERAGE should be set. For e.g. a density field
;              this could be either the particle mass (AVERAGE should
;              not be set) or the density (AVERAGE should be set).
;       POSX:  Array of X coordinates of field samples, unit indices: [0,NX>.
;       NX:    Desired number of grid points in X-direction.
;       
; OPTIONAL INPUTS:
;      POSY: Array of Y coordinates of field samples, unit indices: [0,NY>.
;      NY:   Desired number of grid points in Y-direction.
;      POSZ: Array of Z coordinates of field samples, unit indices: [0,NZ>.
;      NZ:   Desired number of grid points in Z-direction.
;
; KEYWORD PARAMETERS:
;       AVERAGE:    Set this keyword if the nodes contain field samples
;                   (e.g. a temperature field). The value at each grid
;                   point will then be the average of all the samples
;                   allocated to it. If this keyword is not set, the
;                   value at each grid point will be the sum of all the
;                   nodes allocated to it (e.g. for a density field from
;                   a distribution of particles). (D=0). 
;       WRAPAROUND: Set this keyword if the data is periodic and if you
;                   want the first grid point to contain samples of both
;                   sides of the volume (see below). (D=0).
;       NO_MESSAGE: Suppress informational messages.
;
; Example of default NGP allocation: n0=4, *=gridpoint.
;
;     0   1   2   3     Index of gridpoints
;     *   *   *   *     Grid points
;   |---|---|---|---|   Range allocated to gridpoints ([0.0,1.0> --> 0, etc.)
;   0   1   2   3   4   posx
;
; Example of NGP allocation for WRAPAROUND: n0=4, *=gridpoint.
;
;   0   1   2   3         Index of gridpoints
;   *   *   *   *         Grid points
; |---|---|---|---|--     Range allocated to gridpoints ([0.5,1.5> --> 1, etc.)
;   0   1   2   3   4=0   posx
;
;
; OUTPUTS:
;       Prints that a NGP interpolation is being performed of x
;       samples to y grid points, unless NO_MESSAGE is set. 
;
; RESTRICTIONS:
;       All input arrays must have the same dimensions.
;       Position coordinates should be in `index units' of the
;       desired grid: POSX=[0,NX>, etc.
;
; PROCEDURE:
;       Nearest grid point is determined for each sample.
;       Samples are allocated to nearest grid points.
;       Grid point values are computed (sum or average of samples).
;
; EXAMPLE:
;       nx = 20
;       ny = 10
;       posx = randomu(s,1000)
;       posy = randomu(s,1000)
;       value = posx^2+posy^2
;       field = ngp(value,posx*nx,nx,posy*ny,ny,/average)
;       surface,field,/lego
;
; NOTES:
;       Use tsc.pro or cic.pro for a higher order interpolation schemes.    A 
;       standard reference for these interpolation methods is:   R.W. Hockney 
;       and J.W. Eastwood, Computer Simulations Using Particles (New York: 
;       McGraw-Hill, 1981).
; MODIFICATION HISTORY:
;       Written by Joop Schaye, Feb 1999.
;       Check for LONG overflow  P. Riley/W. Landsman   December 1999
;-

nrsamples=n_elements(value)
nparams=n_params()
dim=(nparams-1)/2

IF dim LE 2 THEN BEGIN
    nz=1
    IF dim EQ 1 THEN ny=1
ENDIF
nxny = long(nx)*long(ny)


;---------------------
; Some error handling.
;---------------------

on_error,2  ; Return to caller if an error occurs.

IF NOT (nparams EQ 3 OR nparams EQ 5 OR nparams EQ 7) THEN BEGIN
    message,'Incorrect number of arguments!',/continue
    message,'Syntax: NGP, VALUE, POSX, NX[, POSY, NY, POSZ, NZ,' + $
      ' /AVERAGE, /WRAPAROUND, /NO_MESSAGE]'
ENDIF 

IF (nrsamples NE n_elements(posx)) OR $
  (dim GE 2 AND nrsamples NE n_elements(posy)) OR $
  (dim EQ 3 AND nrsamples NE n_elements(posz)) THEN $
  message,'Input arrays must have the same dimensions!'

IF NOT keyword_set(no_message) THEN $
  print,'Interpolating ' + strtrim(string(nrsamples,format='(i10)'),1) $
  + ' samples to ' + strtrim(string(nxny*nz,format='(i10)'),1) + $
  ' grid points using NGP...'


;-----------------------------
; Compute nearest grid points.
;-----------------------------

IF keyword_set(wraparound) THEN BEGIN
    ; Coordinates of nearest grid point (ngp).
    ngx=fix(posx+0.5)
    ; Periodic boundary conditions.
    bad=where(ngx EQ nx,count)
    IF count NE 0 THEN ngx[bad]=0
    IF dim GE 2 THEN BEGIN 
        ngy=fix(posy+0.5)
        bad=where(ngy EQ ny,count)
        IF count NE 0 THEN ngy[bad]=0
        IF dim EQ 3 THEN BEGIN
            ngz=fix(posz+0.5)
            bad=where(ngz EQ nz,count)
            IF count NE 0 THEN ngz[bad]=0
        ENDIF
    ENDIF 
    bad=0  ; Free memory.
ENDIF ELSE BEGIN
    ; Coordinates of nearest grid point (ngp).
    ngx=fix(posx)
    IF dim GE 2 THEN BEGIN  
        ngy=fix(posy)
        IF dim EQ 3 THEN ngz=fix(posz)
    ENDIF
ENDELSE

; Indices of grid points to which samples are assigned.
CASE dim OF
    1: index=temporary(ngx)
    2: index=temporary(ngx)+temporary(ngy)*nx
    3: index=temporary(ngx)+temporary(ngy)*nx+temporary(ngz)*nxny
ENDCASE


;-------------------------------
; Interpolate samples to grid.
;-------------------------------

field=fltarr(nx,ny,nz)

FOR i=0l,nrsamples-1l DO field[index[i]]=field[index[i]]+value[i]


;--------------------------
; Compute weighted average.
;--------------------------

IF keyword_set(average) THEN BEGIN
    ; Number of samples per grid point.
    frequency=histogram(temporary(index),min=0,max=nxny*nz-1l)
    
    ; Normalize.
    good=where(frequency NE 0,nrgood)
    field[good]=temporary(field[good])/temporary(frequency[good])
ENDIF

return,field

END  ; End of function ngp.