Pro Tvcirc, radius, xc, yc
;+
; NAME:
;     TVCIRC
; PURPOSE:
;     Draw circle(s) of specified radius at specified position(s)
; EXPLANATION: 
;     If a position is not specified, and device has a cursor, then a circle
;     is drawn at the current cursor position.
;
; CALLING SEQUENCE:
;     TVCIRCLE, rad, x, y         
;
; INPUTS:
;     RAD - radius of circle(s) to be drawn, positive numeric scalar
;
; OPTIONAL INPUT:
;      X - x position for circle center, vector or scalar
;      Y - y position for circle center, vector or scalar
;               If X and Y are not specified, and the device has a cursor, 
;               then program will draw a circle at the current cursor position
; OUTPUTS:
;       None
;
   On_Error, 2   ; Return to caller

  if N_params() LT 3 then begin
        if (!D.WINDOW EQ -1) then message, $
                'ERROR - Cursor not available for device ' + !D.NAME
        cursor, xc, yc, /DEVICE, /NOWAIT
        if (xc LT 0) or (yc LT 0) then begin
        message,'Position cursor in window ' + strtrim(!D.WINDOW,2) + $
                ' -- then hit mouse button',/INF
        cursor, xc, yc, /DEVICE, /WAIT
        message,'Circle is centered at (' + strtrim(xc,2) + ',' + $
                strtrim(yc,2) + ')',/INF
  endif

  endif 

    N_circle = min( [ N_elements(xc), N_elements(yc) ] )

    if keyword_set( DATA ) then begin 
                coord = abs(convert_coord(radius,0,/data,/to_dev) - $
                        convert_coord(0,0,/data,/to_dev)) 
                irad =  round( coord[0] )
    endif else irad = round(radius)

   x = 0
   y = irad 
   d = 3 - 2 * irad


   ; Find the x and y coordinates for one eighth of a circle.
   ; The maximum number of these coordinates is the radius of the circle.

   xHalfQuad = Make_Array( irad + 1, /Int, /NoZero )
   yHalfQuad = xHalfQuad

   path = 0

   WHILE x lt y $
   DO BEGIN

      xHalfQuad[path] = x
      yHalfQuad[path] = y

      path = path + 1

      IF d lt 0 $
      THEN d = d + (4*x) + 6 $
      ELSE BEGIN

           d = d + (4*(x-y)) + 10
           y = y - 1

           END

      x = x + 1

      END

   IF x eq y $
   THEN BEGIN ; Fill in last point

        xHalfQuad[path] = x
        yHalfQuad[path] = y

        path = path + 1

        END ; Filling in last point

   ; Shrink the arrays to their correct size

   xHalfQuad = xHalfQuad[ 0:path-1 ]
   yHalfQuad = yHalfQuad[ 0:path-1 ]

   ; Convert the eighth circle into a quadrant

   xQuad = [ xHalfQuad, Rotate(yHalfQuad, 5) ]
   yQuad = [ yHalfQuad, Rotate(xHalfQuad, 5) ]

   ; Prepare for converting the quadrants into a full circle

   xQuadRev = Rotate( xQuad[0:2*path-2], 5 )
   yQuadRev = Rotate( yQuad[0:2*path-2], 5 )

   ; Create full-circle coordinates

   x = [ xQuad, xQuadRev, -xQuad[1:*], -xQuadRev ]
   y = [ yQuad, -yQuadRev, -yQuad[1:*], yQuadRev ]

   ; Plot the coordinates about the given center

   if keyword_set(data) then begin        ;Convert to device coordinates
        coord = convert_coord( xc, yc, /DATA, /TO_DEVICE)
        xcen = round(coord[0,*]) & ycen = round(coord[1,*])
   endif else begin
        xcen = round(xc) & ycen = round(yc)
   endelse

   for i = 0l, N_circle-1 do begin
      if keyword_set(fill) then begin
            polyfill, x+xcen[i],  y + ycen[i], 255, /DEV, $
                _Extra = _extra
      endif else begin
            PlotS, x + xcen[i], y+ ycen[i], 255, /DEV, $
                _Extra = _extra
      endelse

   endfor

   Return
   End; TVcircle