Interpolate a path with equal steps

// Interpolate to create a path with equal spacing between successive points in the XY plane.
// the path described by xwave and ywave as inputs needs to be quite smooth for this to work.
// play with step size (step), number of nodes for interpolation (numNodes) and linear vs. cubic spline
// interpolation (flagT = 1, 2) to optimise the path output as xout and yout.
// step size must be rather small with respect to spacing of path points in xwave, ywave
   
function StepInterpolate(xwave, ywave, step, flagT, numNodes)
    wave xwave, ywave
    variable step, flagT, numNodes
       
    make /o/n=1 xout, yout
    xout=xwave[0]
    yout=ywave[0]
   
    make /o/n=(numNodes) w_nodesX, w_nodesY
    make /o/n=200 w_interp
   
    variable inPnt=0, outPnt=0, minPnt
    variable nextX, nextY

    variable i, inc=(wavemax(xwave)-wavemin(xwave))/1e5
   
    do     
        minPnt=ceil(inPnt-numNodes/2)
        minPnt=max(0, Minpnt)
        minPnt=min(minPnt, numpnts(xwave)-numNodes)
           
        w_nodesX=xwave[minPnt+p]
        w_nodesY=ywave[minPnt+p]
       
        if ( (wavemax(w_nodesY)-wavemin(w_nodesY)) < (wavemax(w_nodesX)-wavemin(w_nodesX)) )
        // heading in a 'horizontal' direction 
           
            Interpolate2/T=(flagT)/N=200/E=2/Y=w_interp w_nodesX, w_nodesY
           
            nextX=xout[outPnt]
            i=0
            do
                i+=1
                if (i>1e5)
                    return 0
                endif
                nextX+=inc*sign(w_nodesX[numpnts(w_nodesX)-1]-w_nodesX[0])
            while (sqrt( (w_interp(nextX)-yout[outPnt])^2 + (nextX-xout[outPnt])^2 ) < step )
            nextY=w_interp(nextX)
            if (nextX<min(xwave[inPnt], xwave[inpnt+1]) || nextX>max(xwave[inPnt], xwave[inpnt+1]) )
                do
                    inPnt+=1
                    if(inPnt>(numpnts(ywave)-1))
                        return 1
                    endif
                while ( (nextX<min(xwave[inPnt], xwave[inpnt+1])) || (nextX>max(xwave[inPnt], xwave[inpnt+1])) )
            endif      
           
        else // heading in a "more vertical" direction
           
            Interpolate2/T=(flagT)/N=200/E=2/Y=w_interp w_nodesY, w_nodesX
           
            nextY=Yout[outPnt]
            i=0
            do
                i+=1
                if (i>1e5)
                    return 0
                endif
                nextY+=inc*sign(w_nodesY[numpnts(w_nodesY)-1]-w_nodesY[0])
            while (sqrt( (nextY-yout[outPnt])^2 + (w_interp(nextY)-xout[outPnt])^2 ) < step )
            nextX=w_interp(nextY)
           
            if (nextY<min(ywave[inPnt], ywave[inpnt+1]) || nextY>max(ywave[inPnt], ywave[inpnt+1]) )
                do
                    inPnt+=1
                    if(inPnt>(numpnts(ywave)-1) )
                        return 2
                    endif                          
                while (nextY<min(ywave[inPnt], ywave[inpnt+1]) || nextY>max(ywave[inPnt], ywave[inpnt+1]) )
            endif              
        endif
               
        outPnt+=1
        insertpoints outPnt, 1, xout, yout
        xout[outpnt]=nextX
        yout[outpnt]=nextY
   
    while(inPnt<numpnts(xwave))
    return 3
end