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