FUNCTION flat_centerofmass, box, centroid   , $
        polar       = polar         , $
        shift_origin= shift_origin  , $
        degrees     = degrees

    @compile_opt.pro            ; On error, return to caller


InitVar, polar        , /key
InitVar, shift_origin , /key
InitVar, get_centroid , /key
add2centroid = IsType(centroid, /defined)

IF polar THEN BEGIN                             ; Polar coordinates

    CASE add2centroid OF
    0: BEGIN

        rpm = ToRadians(degrees=degrees)

        ; box specifies the phase angle between [-!pi,!pi]. The phase angle
        ; covered runs from box[0,0] to box[0,1]. The phase angle for the centroid
        ; is the mean of the two, but some care must be taken when the box
        ; straddles phase angle !pi.

        phi = box[0,*]                          ; Phase angle limits
        phi += (2*!pi/rpm)*[0B, phi[1] LT phi[0]]; In case the box straddles !pi
        phi_centroid = 0.5*total(phi)           ; Phase angle of centroid

        ; The radial location of the centroid follow from an equal area argument:
        ; the area between the lower radius and the centroid must be the same as
        ; between the centroid and the outer radius.

        rad = box[1,*]                          ; Radial limits
        rad_centroid = sqrt(0.5*total(rad*rad)) ; Radius of centroid

        CASE shift_origin OF
        0: result = [AngleRange(phi_centroid, /pi, degrees=degrees), rad_centroid]
        1: result = [phi-phi_centroid, rad-rad_centroid]
        ENDCASE

    END
    1: BEGIN

        ; Call dr1 = box[1] and dr2 = box[3]
        ; Dropping the box at radius r the radial boundaries are:
        ;   r1 = r+dr1, r2 = r+dr2
        ; The centroid of this box is located at s determined by:
        ;   s^2 = (r1^2+r2^2)/2 = r^2+r*(dr1+dr2)+(dr1^2+dr2^2)/2
        ; The left hand side is the value stored in centroid[1]
        ; Solving this for r (taking the positive root) gives
        ;   r = -(dr1+dr2)/2+sqrt(s^2-(dr1-dr2)^2/4)
        ; Dropping the box at r will put the centroid at the required value s.

        rad = sqrt(centroid[1]^2-0.25*(box[1]-box[3])^2) - 0.5*(box[1]+box[3])

        result = [centroid[0],rad]#[1,1]+box
        result[0,*] = AngleRange(result[0,*], /pi, degrees=degrees)

    END
    ENDCASE

ENDIF ELSE BEGIN                                ; Rectangular coordinates

    CASE add2centroid OF
    0: BEGIN
        xy_centroid = 0.5*total(box,2)

        CASE shift_origin OF
        0: result = xy_centroid
        1: result = box-xy_centroid#[1,1]
        ENDCASE
    END
    1: result = centroid#[1,1]+box
    ENDCASE

ENDELSE

RETURN, result  &  END