pro azam_sphtri0, olat, b, c, x, y 
;-----------------------------------------------------------------------------
;
;	subroutine:  azam_sphtri0
;
;	purpose:  Transform a unit vector between two spherical reference
;		  points.  (Same effect as sphtri but does not calculate
;		  parallactic angle).
;
;	reference:  Fundamental spherical trigonometry equations from Smart:
;		    Text-Book on Spherical Astronomy, eq. A, B, C, and D.
;
;	authors:  kcjones@sunspot.noao.edu (Phil Wiborg),   paul@ncar
;			modified 13Jan2012 by BWL for processing vectors
;
;-----------------------------------------------------------------------------
;
;	usage:	call sphtri0_sbsp( olat,  az0, el0,  az1, el1 )
;
;		Transform a unit vector between two spherical
;		reference points.
;
;		Definition:
;
;			sct0	~ The great circle sector that
;				  connects the two reference points
;				  on the sphere
;
;		Input Arguments (radians):
;
;			olat	- !pi/2.-abs(sct0)
;			az0	- Azimuth about first reference
;				- point, from sector sct0
;			el0	- Elevation at first reference
;				  point, positive away from surface
;
;		Output Arguments (radians):
;
;			az1	- Azimuth about second reference
;				  point, from sector sct0
;			el1	- Elevation at second reference
;				  point, positive away from surface
;
;		Note: If az0 is positive CCW, the az1 is positive CCW
;		      If az0 is positive  CW, the az1 is positive  CW
;
;	The routine can be used to convert between coordinate
;	systems as follows:
;
;	(1) Let: az0 = hour angle      (2) Let: az0 = azimuth
;	         el0 = declination              el0 = elevation
;	        olat = earth latitude          olat = earth latitude
;
;	   Then: az1 = azimuth            Then: az1 = hour angle
;	         el1 = elevation                el1 = declination
;
;ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

	sina = sin(olat)
	cosa = cos(olat)

	sinb = sin(b)
	cosb = cos(b)
	sinc = sin(c)
	cosc = cos(c)

	tmp1 = cosc*cosb
	tmp2 = sinc*cosa
	sinx = -cosc*sinb
	cosx = tmp2-tmp1*sina

;  old single value computation
;    siny = max([-1.,min([1.,sinc*sina+tmp1*cosa])])
;  revisions for array calculations
	siny = sinc*sina+tmp1*cosa
	whr = where(siny gt 1.,countc)
	if countc gt 0 then siny(whr) = 1.
	whr = where(siny lt -1.,countc)
	if countc gt 0 then siny(whr) = -1.

	x    = atan(sinx,cosx)
	y    = asin(siny)

	return
	end