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