Wave arithmetic for mismatched XY wave pairs

// Constants to use for the "operation" input
Constant XYArith_Add=0
Constant XYArith_Subtract=1
Constant XYArith_Multiply=2
Constant XYArith_Divide=3
 
Function UnmatchedXYPairArithmetic(axwave, aywave, bxwave, bywave, operation, outwavename)
 
    WAVE axwave, aywave, bxwave, bywave
    Variable operation
    String outwavename
 
    Sort axwave, axwave, aywave
    Sort bxwave, bxwave, bywave
    
    WaveStats/Q axwave
    Variable amin = V_min
    Variable aminloc = V_minloc
    Variable amax = V_max
    Variable amaxloc = V_maxloc
    WaveStats/Q bxwave
    Variable bmin = V_min
    Variable bminloc = V_minloc
    Variable bmax = V_max
    Variable bmaxloc = V_maxloc
    
    if ((amin > bmax) || (amax < bmin))
        DoAlert 0, "X waves don't overlap"
    endif
    Variable beginx = amin < bmin ? bmin : amin
    Variable endx = amax < bmax ? amax : bmax
    Variable beginp = binarysearch(axwave, beginx)+1
    Variable endp = binarysearch(axwave, endx)
    Duplicate/FREE/R=[beginp, endp] axwave,srcax
    Duplicate/FREE/R=[beginp, endp] aywave,srcay
    beginp = binarysearch(bxwave, beginx)+1
    endp = binarysearch(bxwave, endx)
    Duplicate/FREE/R=[beginp, endp] bxwave,srcbx
    Duplicate/FREE/R=[beginp, endp] bywave,srcby
    
    KillWaves/Z $(outwavename+"_X")
    Concatenate/NP/O {srcax, srcbx}, $(outwavename+"_X")/WAVE=xw
    Duplicate/O xw, $outwavename/WAVE=yw
 
    switch(operation)
        case 0:
            yw = aywave[binarysearchinterp(axwave, xw[p])] + bywave[binarysearchinterp(bxwave, xw[p])]
            break;
        case 1:
            yw = aywave[binarysearchinterp(axwave, xw[p])] - bywave[binarysearchinterp(bxwave, xw[p])]
            break;
        case 2:
            yw = aywave[binarysearchinterp(axwave, xw[p])] * bywave[binarysearchinterp(bxwave, xw[p])]
            break;
        case 3:
            yw = aywave[binarysearchinterp(axwave, xw[p])] / bywave[binarysearchinterp(bxwave, xw[p])]
            break;
    endswitch
end
 
// IsMonotonic() returns true if the wave has dwave(x)/dx >= 0 for all x.
Static Function IsMonotonic(WAVE wv)
    
    Differentiate/METH=1/EP=1 wv/D=difwave
    WAVE difwave
    Variable wmin = WaveMin(difwave)
    Variable wmax = WaveMax(difwave)
    
    KillWaves difwave
    return (wmin < 0 && wmax <= 0) || (wmin >= 0 && wmax > 0)
End