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