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