;+
; NAME:
;       INSIDE
;
; PURPOSE:
;       Check if a given point is outside or inside a convex quadrilater.
;
; CALLING SEQUENCE:
;
;	Result = INSIDE(x,y, XV, YV)
;
; INPUTS:
;       XV, YV:	X and Y coordinates of the vertices of poligon
;
;       x,y:	X and Y coordinates of the point to check
;
; OUTPUTS:
;       Result:	0 if outside, 1 if inside or on border.
;
; PROCEDURE:
;	Computes the t1 and t2 parameters of the point in a Gouraud parametrization
;	of the surface of the poligon and check if they are or not between 0 and 1.
;
; MODIFICATION HISTORY:
;          G. Li Causi, Rome Astronomical Observatory: 17-01-2000
;          G. Li Causi, 03-10-2008: fixed a bug on a zero over zero case 
;-

FUNCTION INSIDE, x, y, xv, yv

ON_ERROR, 2             ;Return to caller if error

x=double(x)
y=double(y)
xv=double(xv)
yv=double(yv)

ax = xv(0)-xv(1)+xv(2)-xv(3)
bx = xv(3)-xv(0)
cx = xv(1)-xv(0)
dx = xv(0)

ay = yv(0)-yv(1)+yv(2)-yv(3)
by = yv(3)-yv(0)
cy = yv(1)-yv(0)
dy = yv(0)

p = ax * cy - ay * cx
q = ax * (dy - y) + ay * (x - dx) + bx * cy - by * cx
r = bx * (dy - y) + by * (x - dx)
rad = q*q - 4*p*r

IF p EQ 0 THEN BEGIN

	IF q EQ 0 THEN RETURN,0	ELSE t2_a = - r / q

ENDIF ELSE BEGIN

	IF rad LT 0 THEN RETURN,0 ELSE t2_a = (sqrt(rad) - q) / (2*p)

ENDELSE

kx = ax * t2_a + bx
kx1 = x - cx * t2_a - dx

pr = 1d-6                   ;numerical precision

IF kx EQ 0 THEN BEGIN
  IF ABS(kx1) LT pr THEN BEGIN    ;zero over zero case
    ky = ay * t2_a + by
    ky1 = y - cy * t2_a - dy
    IF ky EQ 0 THEN RETURN,0 ELSE t1_a = ky1 / ky
  ENDIF ELSE RETURN, 0    ;infinite
ENDIF ELSE t1_a = kx1 / kx


IF t1_a LT (-pr) OR t1_a GT (1+pr) OR t2_a LT (-pr) OR t2_a GT (1+pr) THEN RETURN,0 ELSE RETURN,1


END