C+
C NAME:
C	SphereWeight
C PURPOSE:
C	Assign a weight to a point on the surface of the sphere, based on the
C	angular distance of the point relative to a reference point ('origin')
C CATEGORY:
C	Auxilliary
C CALLING SEQUENCE:
	function SphereWeight(DXL,YL0,YL1,DL,NTYPE,CUTOFF,WHALF)
C INPUTS:
C	DXL	real	difference in longitude between data point and origin
C			(degrees; longitude [0,360])
C	YL0	real	latitude of origin (degrees; latitude [-90,90])
C	YL1	real	latitude of data point (degrees)
C	DL	real	half-width of weight function (degrees)
C	NTYPE	integer	defines type of weight function
C			= 1	block function (weight 1 inside, 0 outside DL)
C			= 2	(1/x)*exp(-x**2); x = Ang.dist/DL
C			= 3	exp(-x**2); x = Ang.dist/DL
C			(other functions can be added if necessary)
C	CUTOFF	real	if angular distance > CUTOFF*DL, the weight is set to
C			zero (may be necessary to prevent floating underflow)
C OUTPUTS:
C	SphereWeight	real	weighting factor
C	WHALF		real	weigthing factor at distance DL
C CALLS:
C	Say
C INCLUDE:
	include		'math.h'
C MODIFICATION HISTORY:
C    >	The weight function is defined on the surface of a sphere with unit
C	radius, i.e. is a function of two variables THETA,PHI in any spherical
C	coordinate system defined on the sphere.
C    >	The weight function will generally be a rapidly decreasing function of
C	the angular distance between the origin and a data point
C-
	    real	DXL
	    real	YL0
	    real	YL1
	    real	DL
	    integer	NTYPE
	    real	CUTOFF
	    real	WHALF

	logical		LTYPE
	save		LTYPE
!$omp threadprivate(LTYPE)
	data		LTYPE /.TRUE./

	!-------
	! Angular distance between origin and data points (in deg) using cosine rule

	AngDis = sind(YL0)*sind(YL1)+cosd(YL0)*cosd(YL1)*cosd(DXL)
	AngDis = acosd(AngDis)

	SphereWeight = 0.

	if (NTYPE .eq. 1) then

	    !-------
	    ! Block function: points closer than DL get the same non-zero weight; points
	    !	farther away then DL get zero weight)

	    if (LTYPE) call Say('SphereWeight','I','Fnc','Block function')
	    WHALF = 1.
	    D = AngDis/DL
	    if (D .le. CUTOFF .and. D .le. 1) SphereWeight = 1.

	else if (NTYPE .eq. 2) then

	    !-------
	    ! Weight function drops off as (1/x)*exp(-x**2)

	    if (LTYPE) call Say('SphereWeight','I','Fnc','exp(-angle**2)/angle')
	    WHALF = 1./MATH__E
	    D = AngDis/DL
	    if (D .le. CUTOFF) SphereWeight = exp(-D*D)/max(1.E-3,D)

	else if (NTYPE .eq. 3) then

	    !-------
	    ! Add your own weight function here

	    if (LTYPE) call Say('SphereWeight','I','Fnc','exp(-angle**2)')
	    WHALF = 1./MATH__E			! Weight at half-width distance
	    D = AngDis/DL
	    if (D .le. CUTOFF) SphereWeight = exp(-D*D)

	else

	    if (LTYPE) call Say('SphereWeight','E','Fnc','No weighting function available')

	endif

	LTYPE = .FALSE.

	return
	end