FUNCTION cic,value,posx,nx,posy,ny,posz,nz, \$ AVERAGE=average,WRAPAROUND=wraparound,ISOLATED=isolated, \$ NO_MESSAGE=no_message ;+ ; NAME: ; CIC ; ; PURPOSE: ; Interpolate an irregularly sampled field using Cloud in Cell method ; ; EXPLANATION: ; This function interpolates an irregularly sampled field to a ; regular grid using Cloud In Cell (nearest grid point gets ; weight 1-dngp, point on other side gets weight dngp, where ; dngp is the distance to the nearest grid point in units of the ; cell size). ; ; CATEGORY: ; Mathematical functions, Interpolation ; ; CALLING SEQUENCE: ; Result = CIC, VALUE, POSX, NX[, POSY, NY, POSZ, NZ, ; AVERAGE = average, WRAPAROUND = wraparound, ; ISOLATED = isolated, NO_MESSAGE = 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 weighted average of all the ; samples allocated to it. If this keyword is not ; set, the value at each grid point will be the ; weighted 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 you want the first grid point ; to contain samples of both sides of the volume ; (see below). ; ISOLATED: Set this keyword if the data is isolated, i.e. not ; periodic. In that case total `mass' is not conserved. ; This keyword cannot be used in combination with the ; keyword WRAPAROUND. ; NO_MESSAGE: Suppress informational messages. ; ; Example of default allocation of nearest grid points: 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 CIC interpolation is being performed of x ; samples to y grid points, unless NO_MESSAGE is set. ; ; RESTRICTIONS: ; Field data is assumed to be periodic with the sampled volume ; the basic cell, unless ISOLATED is set. ; All input arrays must have the same dimensions. ; Position coordinates should be in `index units' of the ; desired grid: POSX=[0,NX>, etc. ; Keywords ISOLATED and WRAPAROUND cannot both be set. ; ; PROCEDURE: ; Nearest grid point is determined for each sample. ; CIC weights are computed for each sample. ; Samples are interpolated to the grid. ; Grid point values are computed (sum or average of samples). ; NOTES: ; Use tsc.pro for a higher-order interpolation scheme, ngp.pro for a lower ; order interpolation scheme. A standard reference for these ; interpolation methods is: R.W. Hockney and J.W. Eastwood, Computer ; Simulations Using Particles (New York: McGraw-Hill, 1981). ; EXAMPLE: ; nx=20 ; ny=10 ; posx=randomu(s,1000) ; posy=randomu(s,1000) ; value=posx^2+posy^2 ; field=cic(value,posx*nx,nx,posy*ny,ny,/average) ; surface,field,/lego ; ; MODIFICATION HISTORY: ; Written by Joop Schaye, Feb 1999. ; Avoid integer overflow for large dimensions P.Riley/W.Landsman Dec. 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: CIC, VALUE, POSX, NX[, POSY, NY, POSZ, NZ,' + \$ ' AVERAGE = average, PERIODIC = periodic]' 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 keyword_set(isolated) AND keyword_set(wraparound) THEN \$ message,'Keywords ISOLATED and WRAPAROUND cannot both be set!' 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 CIC...' ;----------------------- ; Calculate CIC weights. ;----------------------- ; Compute weights per axis, in order to reduce memory (everything ; needs to be in memory if we compute all nearest grid points first). ;************* ; X-direction. ;************* ; Coordinates of nearest grid point (ngp). IF keyword_set(wraparound) THEN ngx=fix(posx+0.5) \$ ELSE ngx=fix(posx)+0.5 ; Distance from sample to ngp. dngx=ngx-posx ; Index of ngp. IF keyword_set(wraparound) THEN kx1=temporary(ngx) \$ ELSE kx1=temporary(ngx)-0.5 ; Weight of ngp. wx1=1.0-abs(dngx) ; Other side. left=where(dngx LT 0.0,nrleft) ; samples with ngp to the left. ; The following is only correct if x(ngp)>posx (ngp to the right). kx2=kx1-1 ; Correct points where x(ngp)posy (ngp to the right). ky2=ky1-1 ; Correct points where y(ngp)posz (ngp to the right). kz2=kz1-1 ; Correct points where z(ngp) --> cube length different from EDFW paper). index=kx1+ky1*nx+kz1*nxny cicweight=wx1*wy1*wz1 IF keyword_set(average) THEN BEGIN FOR j=0l,nrsamples-1l DO BEGIN field[index[j]]=field[index[j]]+cicweight[j]*value[j] totcicweight[index[j]]=totcicweight[index[j]]+cicweight[j] ENDFOR ENDIF ELSE FOR j=0l,nrsamples-1l DO \$ field[index[j]]=field[index[j]]+cicweight[j]*value[j] index=kx2+ky1*nx+kz1*nxny cicweight=wx2*wy1*wz1 IF keyword_set(average) THEN BEGIN FOR j=0l,nrsamples-1l DO BEGIN field[index[j]]=field[index[j]]+cicweight[j]*value[j] totcicweight[index[j]]=totcicweight[index[j]]+cicweight[j] ENDFOR ENDIF ELSE FOR j=0l,nrsamples-1l DO \$ field[index[j]]=field[index[j]]+cicweight[j]*value[j] IF dim GE 2 THEN BEGIN index=kx1+ky2*nx+kz1*nxny cicweight=wx1*wy2*wz1 IF keyword_set(average) THEN BEGIN FOR j=0l,nrsamples-1l DO BEGIN field[index[j]]=field[index[j]]+cicweight[j]*value[j] totcicweight[index[j]]=totcicweight[index[j]]+cicweight[j] ENDFOR ENDIF ELSE FOR j=0l,nrsamples-1l DO \$ field[index[j]]=field[index[j]]+cicweight[j]*value[j] index=kx2+ky2*nx+kz1*nxny cicweight=wx2*wy2*wz1 IF keyword_set(average) THEN BEGIN FOR j=0l,nrsamples-1l DO BEGIN field[index[j]]=field[index[j]]+cicweight[j]*value[j] totcicweight[index[j]]=totcicweight[index[j]]+cicweight[j] ENDFOR ENDIF ELSE FOR j=0l,nrsamples-1l DO \$ field[index[j]]=field[index[j]]+cicweight[j]*value[j] IF dim EQ 3 THEN BEGIN index=kx1+ky1*nx+kz2*nxny cicweight=wx1*wy1*wz2 IF keyword_set(average) THEN BEGIN FOR j=0l,nrsamples-1l DO BEGIN field[index[j]]=field[index[j]]+cicweight[j]*value[j] totcicweight[index[j]]=totcicweight[index[j]]+cicweight[j] ENDFOR ENDIF ELSE FOR j=0l,nrsamples-1l DO \$ field[index[j]]=field[index[j]]+cicweight[j]*value[j] index=kx2+ky1*nx+kz2*nxny cicweight=wx2*wy1*wz2 IF keyword_set(average) THEN BEGIN FOR j=0l,nrsamples-1l DO BEGIN field[index[j]]=field[index[j]]+cicweight[j]*value[j] totcicweight[index[j]]=totcicweight[index[j]]+cicweight[j] ENDFOR ENDIF ELSE FOR j=0l,nrsamples-1l DO \$ field[index[j]]=field[index[j]]+cicweight[j]*value[j] index=kx1+ky2*nx+kz2*nxny cicweight=wx1*wy2*wz2 IF keyword_set(average) THEN BEGIN FOR j=0l,nrsamples-1l DO BEGIN field[index[j]]=field[index[j]]+cicweight[j]*value[j] totcicweight[index[j]]=totcicweight[index[j]]+cicweight[j] ENDFOR ENDIF ELSE FOR j=0l,nrsamples-1l DO \$ field[index[j]]=field[index[j]]+cicweight[j]*value[j] index=kx2+ky2*nx+kz2*nxny cicweight=wx2*wy2*wz2 IF keyword_set(average) THEN BEGIN FOR j=0l,nrsamples-1l DO BEGIN field[index[j]]=field[index[j]]+cicweight[j]*value[j] totcicweight[index[j]]=totcicweight[index[j]]+cicweight[j] ENDFOR ENDIF ELSE FOR j=0l,nrsamples-1l DO \$ field[index[j]]=field[index[j]]+cicweight[j]*value[j] ENDIF ENDIF ; Free memory (no need to free any more local arrays, will not lower ; maximum memory usage). index=0 ;-------------------------- ; Compute weighted average. ;-------------------------- IF keyword_set(average) THEN BEGIN good=where(totcicweight NE 0,nrgood) field[good]=temporary(field[good])/temporary(totcicweight[good]) ENDIF return,field END ; End of function cic.