function curvesect,y1,y2,xpt, verbose=verbose, _extra=e
;+
;function	curvesect
;	finds and returns the point(s) of intersection of two curves
;
;syntax
;	xy=curvesect(y1,y2,xpt,verbose=verbose)
;
;parameters
;	y1	[INPUT; required] points of curve Y1(X)
;	y2	[INPUT; required] points of curve Y2(X)
;		* Y1 and Y2 must be defined on the same grid
;	xpt	[INPUT] x-coordinates of Y1 and Y2
;		* if not supplied, uses the array index
;
;keywords
;	verbose	[INPUT] controls chatter
;	_extra	[JUNK] here only to prevent crashing the program
;
;restrictions
;	only returns the intersections within the domain of XPT
;	inputs must be sorted in X
;
;description
;	a very simple algorithm that checks whether two curves
;	have any crossings, and if they do, finds the point of
;	intersection using the line segments that connect across
;	the crossings.  Don't push it.
;
;example
;	.run curvesect
;
;history
;	Vinay Kashyap (2012dec)
;-

;	usage
ok='ok' & np=n_params() & n1=n_elements(y1) & n2=n_elements(y2)
if np lt 2 then ok='Insufficient parameters' else $
 if n1 eq 0 then ok='Y1 is not defined' else $
  if n2 eq 0 then ok='Y2 is not defined' else $
   if n1 ne n2 then ok='Y1 and Y2 must be on same grid' else $
    if n1 lt 2 then ok='curves must have at least 2 points'
if ok ne 'ok' then begin
  print,'Usage: xy=curvesect(y1,y2,xpt,verbose=verbose)'
  print,'  finds and returns points of intersection between Y1 and Y2'
  if np ne 0 then message,ok,/informational
  return,-1L
endif

;	initialize
vv=0L & if keyword_set(verbose) then vv=long(verbose[0])>1L
xx=findgen(n1) & nx=n_elements(xpt)
if nx eq n1 then xx=xpt else if vv gt 0 then message,$
	'XPT incompatible with Y1 and Y2; using array indices',/informational

;	are there any intersections at all?
i1=y1 gt y2
ti1=total(i1)
if ti1 eq 0 or ti1 eq n1 then begin	;(one curve lies entirely above the other
  message,'No intersections in this range',/informational
  return,[!values.F_NAN,!values.F_NAN]
endif					;ti1=0 or ti1=N1)

;	how many intersections?
di1=i1[1:*]-i1
od1=where(di1 ne 0,nsec)
;if mod1 ne nsec then message,'BUG!'

;for each intersection, find the adjacent points and find the intersect of the segment
;	y-y10 = m1 * (x-x0) ==> y = m1*x + (y10-m1*x0)
;	similarly, y = m2*x + (y20-m2*x0)
;	so (m1-m2)*x_int = (y20-m2*x0)-(y10-m1*x0) ==> x_int = ((y20-y10)-x0*(m2-m1))/(m1-m2)
;	and y_int = m1*x_int + (y10-m1*x0)
xy=fltarr(2,nsec)
for i=0L,nsec-1L do begin
  k0=od1[i] & k1=k0+1
  dx=xx[k1]-xx[k0] & dy1=y1[k1]-y1[k0] & dy2=y2[k1]-y2[k0]
  m1=dy1/dx & m2=dy2/dx
  if m1 eq m2 then begin
    xy[*,i]=[xx[k0],y1[k0]]
    continue	;skip to next point
  endif
  xint=((y2[k0]-y1[k0])-xx[k0]*(m2-m1))/(m1-m2)
  yint=0.5*((m1*xint+(y1[k0]-m1*xx[k0])) + (m2*xint+(y2[k0]-m2*xx[k0])))
  xy[*,i]=[xint,yint]
endfor

if vv gt 10000 then stop,'HALTing; type .CON to continue'

return,xy
end

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;example usage

if not keyword_set(verbose) then verbose=1
xpt=findgen(20) & y1=xpt^2-20*x+100. & y2=1*xpt+40
plot,xpt,y1,psym=-1,xtitle='X',ytitle='Y',title='[PINTofALE] curvesect.pro',line=1 & oplot,xpt,y2,psym=-1,line=1

xy=curvesect(y1,y2,xpt,verbose=verbose)

oplot,xy[0,*],xy[1,*],psym=4,symsize=2

end