Polar averaging and errors

Function polarAvgErr(inWave,opStr,degrad)
    Wave inWave //single wave with list of angles
   
    //opStr = "avg", "sdev", or "sem"
    //degrad = "deg" or "rad" depending on the input data type
    String opStr,degrad
    Variable i,size,avgX,avgY,avg,sdev,sem,errX,errY,value,multiplier
   
    size = DimSize(inWave,0)
    Make/FREE/N=(size) xTotal,yTotal
       
    //input is degrees or radians
    strswitch(degrad)
        case "deg":
            xTotal = cos(inWave * pi/180)
            yTotal = sin(inWave * pi/180)
            multiplier = 180/pi
            break
        case "rad":
            xTotal = cos(inWave)
            yTotal = sin(inWave)
            multiplier = 1
            break
    endswitch

    //polar average
    avgX = mean(xTotal)
    avgY = mean(yTotal)
    avg = atan2(avgY,avgX)*multiplier
   
    If(avg < 0)
        avg += 360
    EndIf
   
    //polar standard deviation
    For(i=0;i<size;i+=1)       
        errX += (xTotal[i] - avgX)^2/(size-1)
        errY += (yTotal[i] - avgY)^2/(size-1)
    EndFor
   
    errX = sqrt(errX)      
    errY = sqrt(errY)
    sdev = atan2(errY,errX)*multiplier
           
    //polar standard error of the mean
    sem = sdev/sqrt(size)
   
    //select return parameter
    strswitch(opStr)
        case "avg":
            value = avg
            break
        case "sdev":
            value = sdev
            break
        case "sem":
            value = sem
            break
    endswitch
   
    return value
End