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.;       Postition 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)/2IF dim LE 2 THEN BEGIN    nz=1    IF dim EQ 1 THEN ny=1ENDIFnxny = 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)    ENDIFENDELSE; 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)*nxnyENDCASE;-------------------------------; 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])ENDIFreturn,fieldEND  ; End of function ngp.