Polar densityplot

Function fDemo()

    Make/O/N=(2500, 360) polarMatrix
    SetScale x, 0, 1, "radius", polarMatrix
    SetScale/P y, 0, 1, "degrees", polarMatrix
   
    polarMatrix= p*q
   
    Make/O/N=(2500,2500) matrix
    Variable radiusMax= DimOffset(polarMatrix,0)+DimDelta(polarMatrix,0)*(DimSize(polarMatrix,0)-1)
    SetScale/I x, -radiusMax, radiusMax, "", matrix
    SetScale/I y, -radiusMax, radiusMax, "", matrix
    MultiThread matrix= FindZInPolarMatrix(polarMatrix,x,y)
End

ThreadSafe Function FindZInPolarMatrix(polarMatrix,x,y)
    Wave polarMatrix
    Variable x,y
   
    Variable/C rTheta= r2polar(cmplx(x,y))  // cmplx(radius,angleinRadians)
   
    Variable radius= real(rTheta)
    Variable degrees= imag(rTheta)*180/pi
    if( degrees < 0 )
        degrees += 360  // assuming all angles in polarMatrix are >=0
    endif
    if( degrees > 359 )
        degrees= 359    // assuming all angles in polarMatrix are <=359
    endif
   
    Variable radiusMax= DimOffset(polarMatrix,0)+DimDelta(polarMatrix,0)*(DimSize(polarMatrix,0)-1)
    Variable radiusMin= DimOffset(polarMatrix,0)
    Variable zVal= NaN
    if( (radius >= radiusMin) && (radius <= radiusMax) )
        zVal = polarMatrix(radius)(degrees)
    endif
    return zVal
End