;+
; NAME:
;       OVERLAP
;
; PURPOSE:
;       Computes the area of the intersection of two convex quadrilaters.
;
; CATEGORY:
;
;		Mathematics
;
; CALLING SEQUENCE:
;
;		Result = OVERLAP( XV1, YV1, XV2, YV2 [, frac1] [, frac2] [, DRAW=draw])
;
; INPUTS:
;       XV1, YV1:  X and Y coordinates of the vertices of poligon 1
;
;       XV2, YV2:  X and Y coordinates of the vertices of poligon 2
;
; OUTPUTS:
;       Result:	The area of the intersection of the two quadrilaters.
;
; KEYWORD PARAMETERS:
;
;		DRAW:	If set, show the position of the poligons and the found intersection points
;
; OPTIONAL OUTPUTS:
;       frac1:	The area of the intersection as fraction of the area of poligon 1
;
;       frac2:	The area of the intersection as fraction of the area of poligon 2
;
; PROCEDURE:
;	Computes the intersection of each edge of the two poligons in parametrized form
;	and excludes the points outside the intersection area. Then sort the vertices
;	and computes the area using P_AREA.
;
;
; CALLS:
;	Calls the function INSIDE for checking if a point is inside a quadrilater and
;	the function P_AREA to compute the areas.
;
; MODIFICATION HISTORY:
;       G. Li Causi, Rome Astronomical Observatory: 17-01-2000
;-


FUNCTION OVERLAP, XV1, YV1, XV2, YV2, frac1, frac2, DRAW=draw

ON_ERROR, 2             ;Return to caller if error

XV1 = double(XV1)
XV2 = double(XV2)
YV1 = double(YV1)
YV2 = double(YV2)

sx1 = size(XV1)
sy1 = size(YV1)
sx2 = size(XV2)
sy2 = size(YV2)

IF sx1(0) NE 1 OR sy1(0) NE 1 OR sx2(0) NE 1 OR sy2(0) NE 1 THEN Message,'Array of vertices must be 1-Dimensional'

IF sx1(1) NE sy1(1) OR sx2(1) NE sy2(1) THEN Message,'X and Y arrays of vertices coordinates must be same size'

area_1 = P_AREA(XV1, YV1)
area_2 = P_AREA(XV2, YV2)

IF area_1 EQ 0 OR area_2 EQ 0 THEN BEGIN
	frac1 = 0.
	frac2 = 0.
	RETURN,0.
ENDIF

IF keyword_set(draw) THEN BEGIN
	plot,XV1,YV1, xrange=[min([XV1,XV2])-.1,max([XV1,XV2])+.1], yrange=[min([YV1,YV2])-.1,max([YV1,YV2])+.1]
	oplot,XV2,YV2
	oplot, [XV1[sx1(1)-1], XV1[0]], [YV1[sx1(1)-1], YV1[0]]
	oplot, [XV2[sx2(1)-1], XV2[0]], [YV2[sx2(1)-1], YV2[0]]
ENDIF


;Intersezioni dei lati
xi = fltarr(16)
yi = fltarr(16)
v = 0
FOR v1=0,3 DO BEGIN				;vertice quad.1
	FOR v2=0,3 DO BEGIN				;vertice quad.2

		v11 = (v1+1) MOD 4
		v22 = (v2+1) MOD 4

		ax = XV1(v11) - XV1(v1)
		bx = XV2(v2) - XV2(v22)
		cx = XV2(v2) - XV1(v1)
		ay = YV1(v11) - YV1(v1)
		by = YV2(v2) - YV2(v22)
		cy = YV2(v2) - YV1(v1)

		det = ax*by - ay*bx

		IF det NE 0 THEN BEGIN
			t1 = (cx*by - cy*bx) / det
			t2 = (ax*cy - ay*cx) / det
			IF (t1<t2) GE 0 AND (t1>t2) LE 1 THEN BEGIN
				xi(v) = ax*t1 + XV1(v1)
				yi(v) = ay*t1 + YV1(v1)
				v = v + 1
			ENDIF
		ENDIF

	ENDFOR
ENDFOR


;Vertici dei quadrilateri
FOR v1=0,3 DO BEGIN				;vertice quad.1
	IF INSIDE(XV1(v1), YV1(v1), XV2,YV2) THEN BEGIN
		xi(v) = XV1(v1)
		yi(v) = YV1(v1)
		v = v + 1
	ENDIF
ENDFOR
FOR v2=0,3 DO BEGIN				;vertice quad.2
	IF INSIDE(XV2(v2), YV2(v2), XV1,YV1) THEN BEGIN
		xi(v) = XV2(v2)
		yi(v) = YV2(v2)
		v = v + 1
	ENDIF
ENDFOR

IF v LE 2 THEN BEGIN				;i due poligoni non si sovrappongono
	frac1 = 0.
	frac2 = 0.
	RETURN, 0.
ENDIF

xi = xi(0 : v - 1)
yi = yi(0 : v - 1)


;Ordinamento dei vertici
xbar = TOTAL(xi) / v
ybar = TOTAL(yi) / v

fi = ATAN(yi-ybar,xi-xbar)
ind = sort(fi)
xi = xi(ind)
yi = yi(ind)
ind = uniq(fi(ind))
xi = xi(ind)
yi = yi(ind)

IF n_elements(xi) LE 2 THEN BEGIN				;i due poligoni non si sovrappongono
	frac1 = 0.
	frac2 = 0.
	RETURN, 0.
ENDIF

IF keyword_set(draw) THEN oplot, xi, yi, psym=2

;Risultato
area_i = P_AREA(xi,yi)
frac1 = area_i / area_1
frac2 = area_i / area_2

RETURN, area_i


END