Create Density Matrix From XY Data Points

// NOTE: This is experimental - not thoroughly tested.
//
//	CreateDensityMatrix(xw, yw, numXElements, numYElements, matName)
//	xw and yw contain the location of stars in a galaxy.
//	We want to create a contour plot of the density of stars.
//	To create a contour plot we first need to create a matrix whose X and Y
//	values cover the range of data in xw and yw and whose Z values count
//	the number of stars in a given box (pixel) in the matrix.
//
//	Example:
//		Make/O/N=10000 xData=gnoise(1), yData=gnoise(1)
//		CreateDensityMatrix(xData, yData, 100, 100, "densityMat")
//		NewImage densityMat
//		AppendMatrixContour/T densityMat
Function CreateDensityMatrix(xw, yw, numXElements, numYElements, matName)
	Wave xw, yw
	Variable numXElements			// Number of X elements in output matrix
	Variable numYElements			// Number of Y elements in output matrix
	String matName					// Name to use for output matrix
	
	// Make matrix with specified number of X and Y elements.
	// NOTE: This overwrites an existing wave with the specified name.
	Make/O/N=(numXElements,numYElements) $matName
	Wave mat = $matName			// Create wave reference for the matrix
	mat = 0
	
	WaveStats/Q xw				// Find the X range
	Variable xMin = V_min
	Variable xMax = V_max
	Variable xSpan = xMax - xMin
	
	WaveStats/Q yw
	Variable yMin = V_min
	Variable yMax = V_max
	Variable ySpan = yMax - yMin

	// Set the X and Y scaling of the matrix
	SetScale x xMin, xMax, mat
	SetScale y yMin, yMax, mat
	
	// Create temporary wave containing the X element number for a given point in the input x wave
	Variable numXPoints = numpnts(xw)
	Make/O/N=(numXPoints) xElementNumbersTemp
	xElementNumbersTemp = round(numXElements*(xw-xMin)/xSpan)
	
	// Create temporary wave containing the Y element number for a given point in the input x wave
	Variable numYPoints = numpnts(yw)
	Make/O/N=(numYPoints) yElementNumbersTemp
	yElementNumbersTemp = round(numYElements*(yw-yMin)/ySpan)
	
	Variable xElement, yElement
	for (xElement=0; xElement<numXElements; xElement+=1)
		Variable xEl = xElementNumbersTemp[xElement]
		for (yElement=0; yElement<numXElements; yElement+=1)
			Variable yEl = yElementNumbersTemp[yElement]
			mat[xEl][yEl] += 1
		endfor
	endfor
	
	KillWaves/Z xElementNumbersTemp, yElementNumbersTemp
End