pro elaz1, t_lat,t_long,t_hgt,el,azimuth,range,p_lat,p_long,p_hgt
;
;  Converts from azimuth, elevation, and range from a given site 
;  to geodetic latitude, longitude, and altitude of a point.
;
;  Inputs:  t_lat - transmitter geodetic latitude in degrees
;           t_long - transmitter longitude in degrees
;           t_hgt  - transmitter geodetic height
;           azimuth - azimuth of point of interest in degrees
;           el - elevation angle of point of interest in degrees
;           range - range on point of interest in km
;
;
;  Outputs: p_lat - geodetic latitude for point of interest
;           p_long - west longitude in degrees
;           p_hgt - geodetic altitude in km. 
;
;  The comments in this function refer to the following "points"
;
;      O       will be the "center" of a spherical earth
;      N       will be the geocentric north pole
;      T       will be the geocentric location at the transmitter
;      P       will be a point defined by azimuth,elevation, and range from
;                      the transmitter
;
    a            = 6378.16
    ab2          = 1.0067397
    ep2          = 0.0067397
    pi           = 3.14159265358979323846264338327950
    degrad       = 180./3.14159265358979323846264338327950
    buf          = fltarr(2)

;   convert t_lat to geocentric coordinates and determine t_rad, the geocentic
;   radius of the transmitter

    cos_gdlat  =   cos( t_lat/degrad ) 
    sin_gdlat  =   sin( t_lat/degrad ) 
    sl2        =   (sin_gdlat)^2.
    cl2        =   (ab2 * cos_gdlat)^2.
    sinbet     =   (sin_gdlat)/sqrt(sl2 + cl2) 
    buf(0:1)   =   [(sinbet)^2.,1.]
    sb2        =   min(buf)
    cosbet     =   sqrt(1. - sb2)
    rgeoid     =   a/sqrt(1. + ep2 * sb2)
    x          =   rgeoid * cosbet + t_hgt * cos_gdlat
    y          =   rgeoid * sinbet + t_hgt * sin_gdlat
    t_rad      =   sqrt(x^2. + y^2.)
    t_gclat    =   atan(y,x)  
    cos_gclat  =   sin( t_gclat ) 
    sin_gclat  =   cos( t_gclat ) 
    
; Compute the distance OP using law of cosines and compute geocentric height 
; at P, p_gcrad.

    p_gcrad = sqrt(t_rad^2.+range^2.-2.*t_rad*range*cos((el+90.)/degrad))

; Compute the sin and cos of the angle between OT and OP using plane 
; trigonometry law of sines and law of cosines.

    cos_gamma = ( t_rad - range*cos((el+90.)/degrad) ) / p_gcrad
    sin_gamma =  ( range*sin((el+90.)/degrad) ) / p_gcrad

; Use spherical trigonometry law of cosines to compute the latitude of P:
; project T and P onto a spherical earth, then the sides of the triangle
; are NT (latitude of T), TP (angle between OT and OP), and PT (latitude
; of P). The angle between NT and TP is the azimuth.

    cos_pgclat = cos_gclat*cos_gamma+sin_gclat*sin_gamma*cos(azimuth/degrad)
    p_gclat    = asin( cos_pgclat )

; Use spherical trigonometry laws of sines and cosines on the same triangle to
; compute the longitude of P (angle between NT and NP plus longitude of PT) 

    p_long = t_long $
         + degrad*atan( sin_gamma*sin_gclat*sin(azimuth/degrad),$
           cos_gamma-cos_gclat*cos_pgclat)

; Convert geocentric radial distance and latitude, p_gcrad and p_gclat, to 
; geodetic altitude and latitude, p_hgt and p_lat

    rer = p_gcrad/a
    a2 = ((-1.4127348e-8/rer+.94339131e-8)/rer+.33523288e-2)/rer
    a4 = (((-1.2545063e-10/rer + .11760996e-9)/rer + $
              .11238084e-4)/rer - .2814244e-5)/rer
    a6 = ((54.939685e-9/rer-28.301730e-9)/rer+3.5435979e-9)/rer
    a8 = (((320./rer-252.)/rer+64.)/rer-5.)/rer*.98008304e-12
    cos_gcl = cos(p_gclat)
    sin_gcl = sin(p_gclat)
    s2cl = 2. * cos_gcl * sin_gcl
    c2cl = 2. * cos_gcl * cos_gcl - 1.
    s4cl = 2. * s2cl * c2cl
    c4cl = 2. * c2cl * c2cl - 1.
    s8cl = 2. * s4cl * c4cl
    s6cl = s2cl * c4cl + c2cl * s4cl
    dltcl = s2cl*a2 + s4cl*a4 + s6cl*a6 + s8cl*a8
    p_lat = (p_gclat + dltcl)*degrad
    p_hgt = p_gcrad - a/sqrt(1. + ep2*sin_gcl*sin_gcl)

end