pro intersect,lis1,lis2,y1y2,y1n2,n1y2,val=val,$
	ny1y2=ny1y2,ny1n2=ny1n2,nn1y2=nn1y2,verbose=verbose,$
	_extra=e
;+
;procedure	intersect
;	procedure to obtain common elements of two numeric arrays
;	sorted in ascending order and with no duplicates
;
;syntax
;	intersect,lis1,lis2,y1y2,y1n2,n1y2,/val,$
;	ny1y2=ny1y2,ny1n2=ny1n2,nn1y2=nn1y2,verbose=verbose
;
;parameters
;	lis1	[INPUT; required] input list #1, array
;	lis2	[INPUT; required] input list #2, array
;	y1y2	[OUTPUT; required] elements in LIS1 also in LIS2
;	y1n2	[OUTPUT] elements in LIS1 that are not in LIS2
;	n1y2	[OUTPUT] elements not in LIS1 that are in LIS2
;
;keywords
;	val	[INPUT] if set, returns the actual values in x1x2 rather
;		than the array position numbers
;	ny1y2	[OUTPUT] number of elements in Y1Y2
;	ny1n2	[OUTPUT] number of elements in Y1N2
;	nn1y2	[OUTPUT] number of elements in N1Y2
;	verbose	[INPUT] controls chatter
;	_extra	[JUNK] here only to prevent crashing the program
;
;description
;	Y1Y2 contains indices of elements of LIS1 that are common to both LIS1 and LIS2.
;	Y1N2 contains indices of elements of LIS1 that are not in LIS2.
;	N1Y2 contains indices of elements of LIS2 that are not in LIS1.
;	if Y1Y2/Y1N2/N1Y2 is empty, it is set equal to -1L
;	if keyword VAL is set, actual array elements are returned instead of indices,
;	and if Y1Y2/Y1N2/N1Y2 is empty, it is set equal to !values.F_NAN
;
;history
;	algorithm by mala viswanath
;	idl version (but not optimized for it) by vinay kashyap (8/27/1992)
;	added keyword VAL (3/26/93)
;	added keyword VERBOSE (4/14/02)
;	added keywords NY1Y2,NY1N2,NN1Y2, cleaned up edge case outputs, check for sorted arrays (VK; 7/18/15)
;-

;	usage
ok='ok' & np=n_params() & n1=n_elements(lis1) & n2=n_elements(lis2)
if np lt 2 then ok='Insufficient input parameters' else $
 if n1 eq 0 then ok='LIS1 is undefined' else $
  if n2 eq 0 then ok='LIS2 is undefined' else $
   if n1 eq 1 then ok='LIS1 has only one element, use something else' else $
    if n2 eq 1 then ok='LIS2 has only one element, use something else' else $
     if not arg_present(y1y2) then ok='Y1Y2 absent, why bother running this?'
if ok ne 'ok' then begin
  print,'Usage: intersect,lis1,lis2,y1y2,y1n2,n1y2,/val,$'
  print,'       ny1y2=ny1y2,ny1n2=ny1n2,nn1y2=nn1y2,verbose=verbose'
  print,'  returns common and non-common elements of 2 numeric lists sorted in ascending order'
  if np ne 0 then message,ok,/info
  return
endif

;	keywords
getval = 0 & if keyword_set(val) then getval = 1
vv = 0 & if keyword_set(verbose) then vv=long(verbose(0)) > 1

y1y2 = -1L & y1n2 = y1y2 & n1y2 = y1y2
if getval then y1y2 = y1y2*!values.F_NAN
if getval then y1n2 = y1n2*!values.F_NAN
if getval then n1y2 = n1y2*!values.F_NAN

i = 0 & j = 0 & k = 0 & i1 = 0 & i2 = 0

if vv gt 0 then begin
  print,'# elements in LIS1: '+strtrim(n1,2)
  print,'# elements in LIS2: '+strtrim(n2,2)
endif

z1=lis1 & z2=lis2
;	check that the arrays are sorted, and if not, quit with error
dtmp1=z1[1:*]-z1 & dtmp2=z2[1:*]-z2
oy1=where(dtmp1 lt 0,moy1) & oy2=where(dtmp2 lt 0,moy2)
if moy1 gt 0 then begin & message,'LIS1 must be sorted in ascending order',/informational & return & endif
if moy2 gt 0 then begin & message,'LIS2 must be sorted in ascending order',/informational & return & endif
;	(some day, implement local sorting to avoid above problem)
;ii1=lindgen(n1) & ii2=lindgen(n2)
;os1=sort(lis1) & z1=lis1[os1] & jj1=ii1[os1] & ks1=sort(jj1) & k1=(ii1[jj1])[ks1]
;os2=sort(lis2) & z2=lis2[os2] & jj2=ii2[os2] & ks2=sort(jj2) & k2=(ii2[jj2])[ks2]

;first, get y1y2
i = 0 & i1 = 0 & i2 = 0
while i1 lt n1 and i2 lt n2 do begin
  if vv gt 5 and i eq 1000L*long(i/1000) then print,form="($,a)",'.'
  x = z1(i1) - z2(i2)
  case 1 of
    x lt 0: i1=i1+1L
    x gt 0: i2=i2+1L
    else: begin
      if i eq 0 and getval eq 0 then y1y2 = [ i1 ]
      if i eq 0 and getval eq 1 then y1y2 = [ z1(i1) ]
      if i gt 0 and getval eq 0 then y1y2 = [y1y2,i1]
      if i gt 0 and getval eq 1 then y1y2 = [y1y2,z1(i1)]
      i1=i1+1L & i2=i2+1L & i=i+1L
    end
  endcase
endwhile
ny1y2=n_elements(y1y2)
if ny1y2 eq 1 then begin
  if getval and finite(y1y2[0]) eq 0 then ny1y2=0L
  if not getval and y1y2[0] lt 0 then ny1y2=0L
endif
if vv gt 2 then message,'Y1Y2: done',/info

if np eq 3 then return
;now, get y1n2
i=0 & i1 = 0 & i2 = 0 & flag = 1
while flag do begin
  if vv gt 5 and i eq 100L*long(i/100) then print,form="($,a)",strtrim(i,2)
  x = z1(i1) - z2(i2)
  case 1 of
    x eq 0: begin
      i1=i1+1L & i2=i2+1L
    end
    x lt 0: begin
      if i eq 0 and getval eq 0 then y1n2 = [ i1 ]
      if i eq 0 and getval eq 1 then y1n2 = [ z1(i1) ]
      if i gt 0 and getval eq 0 then y1n2 = [y1n2,i1]
      if i gt 0 and getval eq 1 then y1n2 = [y1n2,z1(i1)]
      i1=i1+1L & i=i+1L
    end
    else: i2=i2+1L
  endcase
  if i2 eq n2 and i1 lt n1 then begin
    flag = 0
    if i eq 0 and getval eq 0 then y1n2 = indgen(n1)
    if i eq 0 and getval eq 1 then y1n2 = z1
    if i gt 0 and getval eq 0 then y1n2 = [y1n2,indgen(n1-i1)+i1]
    if i gt 0 and getval eq 1 then y1n2 = [y1n2,z1(i1:*)]
  endif
  if i1 eq n1 then flag = 0
endwhile
ny1n2=n_elements(y1n2)
if ny1n2 eq 1 then begin
  if getval and finite(y1n2[0]) eq 0 then ny1n2=0L
  if not getval and y1n2[0] lt 0 then ny1n2=0L
endif
if vv gt 2 then message,'Y1N2: done',/info

if np eq 4 then return
;finally, get n1y2
i=0 & i1 = 0 & i2 = 0 & flag = 1
while flag do begin
  if vv gt 5 and i eq 100L*long(i/100) then print,form="($,a)",strtrim(i,2)
  x = z2(i2) - z1(i1)
  case 1 of
    x eq 0: begin
      i2=i2+1L & i1=i1+1L
    end
    x lt 0: begin
      if i eq 0 and getval eq 0 then n1y2 = [ i2 ]
      if i eq 0 and getval eq 1 then n1y2 = [ z2(i2) ]
      if i gt 0 and getval eq 0 then n1y2 = [n1y2,i2]
      if i gt 0 and getval eq 1 then n1y2 = [n1y2,z2(i2)]
      i2=i2+1L & i=i+1L
    end
    else: i1=i1+1L
  endcase
  if i1 eq n1 and i2 lt n2 then begin
    flag = 0
    if i eq 0 and getval eq 0 then n1y2 = indgen(n1)
    if i eq 0 and getval eq 1 then n1y2 = z2
    if i gt 0 and getval eq 0 then n1y2 = [n1y2,indgen(n2-i2)+i2]
    if i gt 0 and getval eq 1 then n1y2 = [n1y2,z2(i2:*)]
  endif
  if i2 eq n2 then flag = 0
endwhile
nn1y2=n_elements(n1y2)
if nn1y2 eq 1 then begin
  if getval and finite(n1y2[0]) eq 0 then nn1y2=0L
  if not getval and n1y2[0] lt 0 then nn1y2=0L
endif
if vv gt 2 then message,'N1Y2: done',/info

return
end