FUNCTION Distance2Line, line, pointx, pointy, grid=grid, close_area=close_area
    @compile_opt.pro        ; On error, return to caller

close_area = keyword_set(close_area)
grid = keyword_set(grid)

r1 = line                       ; Don't change input arrays

IF n_params() EQ 2 THEN $
    pp = pointx         $

ELSE IF NOT grid THEN BEGIN
    nx = n_elements(pointx) < n_elements(pointy)
    pp = [reform(pointx[0:nx-1],1,nx),reform(pointy[0:nx-1],1,nx)]
    IF nx EQ 1 THEN pp = reform(pp)

ENDIF ELSE BEGIN
    nx = n_elements(pointx)
    ny = n_elements(pointy)
    pp = replicate(pointx[0],2,nx,ny)
    pp[0,*,*] = SuperArray(pointx,/trail,ny)
    pp[1,*,*] = SuperArray(pointy,/lead ,nx)

ENDELSE

sz = size(pp)

nl = n_elements(r1)/2           ; # lines
np = n_elements(pp)/2           ; # points to be tested

r1 = reform(r1,2,nl,/overwrite)
IF np GT 1 THEN pp = reform(pp,2,np,/overwrite)

r2 = shift(r1,0,-1)

IF NOT close_area THEN BEGIN
    nl = nl-1
    r1 = r1[*,0:nl-1]           ; array[2,nl] (array[2] if nl=1)
    r2 = r2[*,0:nl-1]
ENDIF

r21 = r2-r1
r21 = sqrt( total(r21*r21,1) )  ; array[nl] (scalar if nl=1)

IF nl GT 1 THEN $
    pp  = SuperArray(pp , nl, after=1)  ; 2 x nl x np array

IF np GT 1 THEN BEGIN
    r1  = SuperArray(r1 , np, /trail)   ; 2 x nl x np array
    r2  = SuperArray(r2 , np, /trail)   ; 2 x nl x np array
    r21 = SuperArray(r21, np, /trail)   ;     nl x np array
ENDIF

r1 = r1-pp
r2 = r2-pp

pp  = (SubArray(r2,elem=0)*SubArray(r1,elem=1)-SubArray(r2,elem=1)*SubArray(r1,elem=0))/r21

; Remove leading dimension of 2 from 'point'

IF np EQ 1 THEN             $
    sz = [0,sz[sz[0]+1],1]  $
ELSE                        $
    sz = [sz[0]-1,sz[2:sz[0]],sz[sz[0]+1],sz[sz[0]+2]/2]

; Add leading dimension of nl elements

IF nl GT 1 THEN BEGIN
    IF np EQ 1 THEN                 $
        sz = [1,nl,sz[sz[0]+1],nl]  $
    ELSE                            $
        sz = [sz[0]+1,nl,sz[1:sz[0]],sz[sz[0]+1],nl*sz[sz[0]+2]]
ENDIF

SyncDims, pp, size=sz

RETURN, pp  &  END