pro dist_ellipse,im,n,xc,yc,ratio,pos_ang, DOUBLE = double
;+
; NAME:
;       DIST_ELLIPSE
; PURPOSE:
;       Create a mask array useful for elliptical aperture photometry
; EXPLANATION:
;       Form an array in which the value of each element is equal to the
;       semi-major axis of the ellipse of specified center, axial ratio, and 
;       position  angle, which passes through that element.  Useful for 
;       elliptical aperture photometry.
;
; CALLING SEQUENCE:
;       DIST_ELLIPSE, IM, N, XC, YC, RATIO, [ POS_ANG] , /DOUBLE
;
; INPUTS:
;       N = either  a scalar specifying the size of the N x N square output
;               array, or a 2 element vector specifying the size of the
;               M x N rectangular output array.
;       XC,YC - Scalars giving the position of the ellipse center.   This does
;               not necessarily have to be within the image
;       RATIO - Scalar giving the ratio of the major to minor axis.   This 
;               should be greater than 1 for position angle to have its 
;               standard meaning.
;
; OPTIONAL INPUTS:
;       POS_ANG - Position angle of the major axis in degrees, measured counter-clockwise
;               from the Y axis.  For an image in standard orientation 
;               (North up, East left) this is the astronomical position angle.
;               Default is 0 degrees.
;
; OPTIONAL INPUT KEYWORD:
;       /DOUBLE - If this keyword is set and nonzero, the output array will
;               be of type DOUBLE rather than floating point.
;
; OUTPUT:
;       IM - REAL*4 elliptical mask array, of size M x N.  THe value of each 
;               pixel is equal to the semi-major axis of the ellipse of center
;                XC,YC, axial ratio RATIO, and position angle POS_ANG, which 
;               passes through the pixel.
;
; EXAMPLE:
;       Total the flux in a elliptical aperture with a major axis of 3', an
;       axial ratio of 2.3, and a position angle of 25 degrees centered on 
;       a specified RA and DEC.   The image array, IM is 200 x 200, and has 
;       an associated FITS header H.
;
;       ADXY, H, ra, dec, x, y       ;Get X and Y corresponding to RA and Dec
;       GETROT, H, rot, cdelt        ;CDELT gives plate scale degrees/pixel
;       cdelt = abs( cdelt)*3600.    ;CDELT now in arc seconds/pixel
;       DIST_ELLIPSE, ell, 200, x, y, 2.3, 25  ;Create a elliptical image mask
;       ell = ell*cdelt[0]           ;Distances now given in arcseconds
;       good = where( ell lt 180 )   ;Within 3 arc minutes
;       print,total( im(good) )      ;Total pixel values within 3'
;
; RESTRICTIONS:
;       The speed of DIST_ELLIPSE decreases and the the demands on virtual
;       increase as the square of the output dimensions.   Users should
;       dimension the output array as small as possible, and re-use the
;       array rather than re-calling DIST_ELLIPSE
;
; REVISION HISTORY:
;       Written    W. Landsman             April, 1991
;       Somewhat faster algorithm          August, 1992
;       Allow rectangular output array     June, 1994
;       Added /DOUBLE keyword   W. Landsman   July 2000
;       Make POS_ANG optional, as documented  W. Landsman Aug 2015
;-
 On_error,2                             ;Return to caller

 if N_params() LT 5 then begin
    print,'Syntax - DIST_ELLIPSE, im, n, xc, yc, ratio, [pos_ang], /DOUBLE'
    print,'   im - output elliptical mask image array'
    print,'   n -  size of output image mask, scalar or 2 element vector'
    print,'   xc,yc - coordinates of ellipse center, scalars'
    print,'   ratio - ratio of major to minor axis of ellipse, scalar'
    print,'   pos_ang - position angle, counterclockwise from up'
    return
 endif
                                          ;Check some parameters
 if N_elements(ratio) NE 1 then message, $
     'ERROR - Axial ratio (fifth parameter) must be a scalar value'

 if N_elements(pos_ang) GT 1 then message, $
     'ERROR - Position angle (sixth parameter) must be a scalar value'

 if N_elements(pos_ang) EQ 0 then pos_ang = 0
 ang = pos_ang /!RADEG                      ;Convert to radians
 cosang = cos(ang)
 sinang = sin(ang)

 if N_elements(N) EQ 2 then begin
        nx = n[0]
        ny = n[1] 
 endif else if N_elements(N) EQ 1 then begin
        ny = n
        nx = n                    ;Make a row
 endif else message, $
        'ERROR - Output size parameter N must contain 1 or 2 elements'
        
 if keyword_set(double) then begin
    x = dindgen(nx) - xc
    y = dindgen(ny) - yc
    im = dblarr(nx, ny, /NOZERO)
 endif else begin
    x = findgen( nx ) - xc
    y = findgen( ny ) - yc
    im = fltarr( nx, ny, /NOZERO )
 endelse
                         ;Rotate pixels to match ellipse orientation
 xcosang = x*cosang
 xsinang = x*sinang

 for i = 0,ny-1 do begin
   xtemp =  xcosang + y[i]*sinang
   ytemp = -xsinang + y[i]*cosang
   im[0,i] = sqrt( (xtemp*ratio)^2 + ytemp^2 )
 endfor

 return
 end