# Returns polynomial coefficients of eight Wiggins' basis function at level #i of
# type 'itype':
#	i = 0 for the surface, etc.
#	itype=0 for 'derivative' function, and 1 for 'value' function.

function [ above below ] = bfun( i, itype )

global z;
global dz;
global nlayers;

  if ( itype == 0 )		# derivative

 	if ( i == 1 )
		above = [ 1 -1 0 0 ]; 	# * z(i);
  	else
		above = [ 1 -1 0 0 ];   # * dz(i);
	endif

  	if ( i < nlayers )
		below = [ 1 -2 1 0 ]; #  * dz(i+1);
  	else
		below = zeros(1,4);
  	endif	

	return;
  endif

 	# otherwise value

  if ( i == 0 )
		above = zeros(1,4);
		below = [ 2 -3 0 1 ]; # * z(1);
		return;
  endif

  above = [ -2 3 0 0 ];

  if ( i < nlayers )
  	below = [ 2 -3 0 1 ];
  else
	below = zeros(1,4);
  endif

endfunction