2D interpolation

Function bilinearFromZMat(minx,dx,nx,miny,dy,ny,xWave,yWave,zMatrix)
    Variable minx,dx,nx,miny,dy,ny
    Wave xWave,yWave,zMatrix
    // first check that x&y waves have one more point than zMatrix:
    if(DimSize(xWave,0)-DimSize(zMatrix,0)!=1)
        doAlert 0,"Input x-wave should have 1 more point than the first dimension of zMatrix"
        return nan
    endif
    if(DimSize(yWave,0)-DimSize(zMatrix,1)!=1)
        doAlert 0,"Input y-wave should have 1 more point than the second dimension of zMatrix"
        return nan
    endif
   
    // Create the output wave:
    Make/O/N=(nx,ny)/Y=(WaveType(zMatrix)) M_InterpolatedValues
    SetScale/P x (minx),(dx),"", M_InterpolatedValues
    SetScale/P y (miny),(dy),"", M_InterpolatedValues
   
    M_InterpolatedValues=findZValue(x,y,xWave,yWave,zMatrix)
End


Function findZValue(inX,inY,xWave,yWave,zMatrix)
    Variable inX,inY
    Wave xWave,yWave,zMatrix
   
    findLevel/Q xWave,inX
    if(V_flag!=0 || GetRTError(1))
        return NaN
    endif
    Variable ax=V_Levelx
    Variable aa=trunc(ax)
    Variable alpha=ax-aa
    findLevel/Q yWave,inY
    if(V_flag!=0 || GetRTError(1))
        return NaN
    endif
    Variable ay=V_Levelx
    Variable bb=trunc(ay)
    Variable beta=ay-bb

    Variable outVal=(1-alpha)*(1-beta)*zMatrix[aa][bb]
            outVal+=(1-alpha)*beta*zMatrix[aa+1][bb]
            outVal+=alpha*(1-beta)*zMatrix[aa][bb+1]
            outVal+=alpha*beta*zMatrix[aa+1][bb+1]
    return outVal
End