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