FUNCTION NEWBILINEAR,P,IX,JY,FLAG
; $Id: bilinear.pro,v 1.4 1994/08/18 21:52:59 chadwick Exp $
;+
; NAME:
;	NEWBILINEAR
; PURPOSE:
;	Allows user to determine bilinear interpolate at a set of reference
;	points.
; CALLING SEQUENCE:
;	Z=NEWBILINEAR(P,X,Y,FLAG)
; INPUT:                 
;	P, FLTARR(N0,M0), the data array, it must be two dimensional.
;	X and Y contain the "virtual subscripts" of P to look up values
;	  for the output.  
;	X, two choices:
;	     1)	X=FLTARR(N), subscripts to look up in P.  Same
;		  set of subscripts is used for all rows in the output array.
;	     2)	X=FLTARR(N,M), also contains "x-axis" subscripts, except
;	       	  specified at all points in the output array.  
;	    in either case X must satisfy,
;		    0 <= MIN(X) < N0  and 0 < MAX(X) <= N0
;	   where N0 is the total number of subscripts in the first dimension
;	   of P.
;	Y, two choices:
;	     1) Y=FLTARR(M), "y-axis" subscript to look up in P.  Same mask
;		  is used in all columns in the output.
;	     2) Y=FLTARR(N,M), at all points in output.
;	    in either case Y must satisfy,
;		    0 <= MIN(Y) < M0  and 0 < MAX(Y) <= M0
;	    where M0 is the total number of subscripts in the second dimension
;	    of P.
;  It is better to use two dimensional arrays for X and Y when calling
;  BILINEAR because the algorithm is somewhat faster.  If X and Y are 
;  one-dimensional then they are converted to two-dimensional arrays on
;  return from the function.  They can then be reused on subsequent calls
;  and thereby save time.  The 2-D array P is unchanged on return.
;        FLAG, point value for which the point should have no weight in the interpolation
; OUTPUT:
;	Z, FLTARR(N,M), the interpolated array.  
; SIDE EFFECTS:
;	can be time consuming.
; RESTRICTIONS:
;	none.
; EXAMPLE:
;	suppose P=FLTARR(3,3) and X=[.1,.2], Y=[.6,2.1]  then Z(0,0) will
;	be returned as though it where equal to P(.1,.6) interpolated from
;	the nearest neighbors at P(0,0),P(1,0),P(1,1) and P(0,1).
; PROCEDURE:
;	invokes bilinear interpolation algorithm to evaluate each element
;	in Z at virtual coordinates contained in X and Y with the data
;	in P.                                                          
; REVISION HISTORY:
;	Nov. 1985	written L.Kramer (U. of Maryland/U. Res. Found.)
;       Nov. 1996       Bill Snyder(UCAR) added the FLAG parameter
;-
	ON_ERROR,2              ;Return to caller if an error occurs	
	IF((N_ELEMENTS(IX) EQ 0) AND (N_ELEMENTS(JY) EQ 0)) THEN BEGIN
	  I=FIX(IX) & J=FIX(JY) & IP=I+1 & JP=J+1
	  DX=IX-FLOAT(I) & DY=JY-FLOAT(J)
	  DX1=(1.-DX) & DY1=(1.-DY) 
	  RETURN,( P(I,J)*DX1*DY1 + P(I,JP)*DX1*DY $
	  	+ P(IP,J)*DX*DY1 + P(IP,JP)*DX*DY)
	ENDIF

	A=SIZE(IX)  & B=SIZE(JY)
	NX=A(1)
	IF(B(0) EQ 1) THEN BEGIN
	  NY=B(1)
	ENDIF ELSE BEGIN
	  NY=B(2)
	ENDELSE                                                   
        IF(A(0) EQ 1) THEN BEGIN
	  TEMP=IX
	  IX=FLTARR(NX,NY)
	  FOR I=0L,NY-1 DO IX(0,I)=TEMP
	ENDIF
	IF(B(0) EQ 1) THEN BEGIN
	  TEMP=JY
	  JY=FLTARR(NY,NX)
	  FOR I=0L,NX-1 DO JY(0,I)=TEMP
	  JY=TRANSPOSE(JY)
	ENDIF
	I=FIX(IX) & J=FIX(JY)
	IP=I+1   &  JP=J+1
	DX=IX-I & DY=JY-J
	DX1=1.-DX & DY1=1.-DY
	Z=FLTARR(N_ELEMENTS(I(*,0)),N_ELEMENTS(J(0,*)))
	NUMX=N_ELEMENTS(I)
	PZ=FLTARR(N_ELEMENTS(P(*,0))+1,N_ELEMENTS(P(0,*))+1)
	PZ(0,0)=P(0:*,0:*)
	TMP=[0.0, 0.0, 0.0, 0.0]
	FOR N=0L,NUMX-1 DO BEGIN
;	  TMP(0) = (PZ(I(N), J(N)) NE FLAG)*PZ(I(N), J(N))*DX1(N)*DY1(N)	
;	  TMP(1) = (PZ(I(N), JP(N)) NE FLAG)*PZ(I(N),JP(N))*DX1(N)*DY(N)	
;	  TMP(2) = (PZ(IP(N), J(N)) NE FLAG)*PZ(IP(N),J(N)) *DX(N) *DY1(N)	
;	  TMP(3) = (PZ(IP(N), JP(N)) NE FLAG)*PZ(IP(N),JP(N))*DX(N) *DY(N)
          IF(PZ(I(N), J(N)) EQ FLAG) THEN TMP(0) = 0.0 ELSE TMP(0)=PZ(I(N), J(N))*DX1(N)*DY1(N)
          IF(PZ(I(N), JP(N)) EQ FLAG) THEN TMP(1) = 0.0 ELSE TMP(1)=PZ(I(N), JP(N))*DX1(N)*DY(N)
	  IF(PZ(IP(N), J(N)) EQ FLAG) THEN TMP(2) = 0.0 ELSE TMP(2)=PZ(IP(N), J(N))*DX(N)*DY1(N)
	  IF(PZ(IP(N), JP(N)) EQ FLAG) THEN TMP(3) = 0.0 ELSE TMP(3)=PZ(IP(N), JP(N))*DX(N)*DY(N)
	  Z(N)=TMP(0)+TMP(1)+TMP(2)+TMP(3)
	ENDFOR
;       flag_locs is the location in the array of 0 values.
;       because it is the result of the where function it 
;       is always an array even if only a -1 is returned.
	flag_locs=where(z eq 0)
	if (flag_locs(0) gt -1) then begin
	   Z(flag_locs) = FLAG
	endif

	RETURN,Z
END